function [gn,gp] = DLD_Boundaries(N, L_bulk, W_bulk, t, PL, Pa, R_cut)
% Calculates the negative and positive boundary gaps for a DLD array
% Inputs are 
% N: inverse of row shift fraction. Number of rows before repeat
% L_Bulk : length of the array unit cell (in flow direction) in microns
% W_Bulk : width of array unit cells (lateral to flow) in microns
% t : array depth (in microns)
% PL: post diameter in lateral direction (direction of primary gap) (in microns)
% Pa: post diameter in axial direction (direction of flow) (in microns)
% R_cut : how much to cut the axial radius of the two pillars at the Nth positive gap

GL=W_bulk-PL;
Ga=L_bulk-Pa;
n=[1:N];
mu=1;

%%%%%%%%%%%%% fitting parameters
a=1.702;
b=2.682;
c=1.833;
%%%%%%%%%%%%%

[Ra_bulk, RL_bulk]=ArrayResistance(W_bulk,L_bulk,Pa,PL,t,mu);

%negative boundaries
Wn=W_bulk*ones(1,N);
Rn=Ra_bulk./n*N;
for i=1:N
    func =@(W) (10^6)^(-3)*12*mu*L_bulk/W/t^3 * (1+ a*(0.6+tan(pi()/2*PL/W))^b   *(t/W)^c) -Rn(i);
     Wn(i)=fzero(func,[1.05*PL 2*W_bulk]);
end
gn=Wn-PL;

%positive boundaries
Wp=W_bulk*ones(1,N);
Rp=Ra_bulk./(2-n/N);

%calculate widths for all but the Nth gap.
for i=1:N-1
    func =@(W) (10^6)^(-3)*12*mu*L_bulk/W/t^3 * (1+ a*(0.6+tan(pi()/2*PL/W))^b   *(t/W)^c) -Rp(i);
     Wp(i)=fzero(func,[1.05*PL 2*W_bulk]);
end
gp=Wp-PL;

%calculate resistance of Nth lateral gap after cut is made
[RNP_axial, RNP_lateral] = ArrayResistance(W_bulk,L_bulk,Pa-R_cut,PL,t,mu);
%calculate required resistance of Nth axial gap
Rp(N) = Ra_bulk +RL_bulk/N-RNP_lateral;
%calculate gap that will give this resistance
func =@(W) (10^6)^(-3)*12*mu*L_bulk/W/t^3 * (1+ a*(0.6+tan(pi()/2*PL/W))^b   *(t/W)^c) -Rp(N);
     Wp(N)=fzero(func,[1.05*PL 2*W_bulk]);
gp(N)=Wp(N)-PL;
gp=gp'/GL;
gn=gn'/GL;

    
    function [Ra, RL] = ArrayResistance(W,L, Pa, PL, t, mu)
%   ArrayResistance calculates the resistance of the unit cell of a DLD array
%   Resistance is return in units of Ns/m^5
%   It uses a function that is fit to modelled data points and covers the
%   range of 
%       0.5 < t/W < 4.5
%       0.3 < P/W <0.9
%   Input dimensions are all in microns
%   W = unit width (lateral direction), L = unit length (axial direction), 
%   PL = Post width in lateral direction, Pa = Post width
%   in axial (flow) direction, t = array depth, mu = fluid viscosity.
GL=W-PL;
Ga=L-Pa;
% model fitting parameters
a=1.702;
b=2.682;
c=1.833;
Ra=(10^6)^(-3)*12*mu*L/W/t^3 * (1+ a*(0.6+tan(pi()/2*PL/W))^b   *(t/W)^c);
RL=(10^6)^(-3)*12*mu*L/W/t^3 * (1+ a*(0.6+tan(pi()/2*Pa/W))^b   *(t/W)^c);
end

end


