%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Filename:      'Microchannel_cell_constants_for_coplanar_electrodes.m'
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Author:         Pontus Linderholm
%       Date:           01 June 2005  
%       Affiliation:    LMIS4, EPFL, Lausanne, Switzerland    
%         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Description:    Calculates the cell constant between two coplanar
%       electrodes in a microchannel. The cell constant is the proportionality 
%       constant between measured resistance and the resistivity of the solution.
%                   
%                       R = rho * cell_constant / channel_width    [ohm m * {} / m]
%       
%       where rho is the resistivity of the solution, cell_constant is the 2D-cell 
%       constant and channel_width is the width of the channel
%       (in a direction perpendicular to the flow direction).
%
%       NOTE that the algorithm assumes that the microchannel has a
%       rectangular cross-section (surrounded by insulating material), or
%       that it is infinitely wide. 
%
%       Geometry:
%
%
%                        ELECTRODE 1      GAP      ELECTRODE 2                          
%          \     \                                              \     \     \
%          xxxxxxxxxxxxx_______________xxxxxxxxxx_______________xxxxxxxxxxxxx   _           
%                                                                                |
%           ->                                                          ->       |             
%      Flow -->                                                          -->     |  CHANNEL
%           ->                                                          ->       |             
%                                                                               _|
%          xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%          /     /     /     /     /     /     /     /     /     /     /
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       INPUT PARAMETERS
%            electode_width             the width of the electrode along the flow direction    
%            gap_width                  the distance between the electrodes    
%            channel_height             the distance from the electrodes to the opposite wall              
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       OUTPUT PARAMETER    
%            cell_constant              the 2D cell constant, which must be
%                                       divided by the width of the channel to obtain the 
%                                       3D cell constant. (a wider channel yields of course 
%                                       a lower resistance)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       REFERENCES
%                   (1)  LINDERHOLM, P. and RENAUD, P.: 'Comment on "AC frequency 
%                          characteristics of coplanar impedance sensors as design parameters" 
%                          by Jongin Hong, et al.', Lab Chip, 2005. 5(12):1416-1417
%                               - presents the mathematical derivation used here
%                   (2)  HILBERG, W.: 'From Approximations to Exact Relations for 
%                          Characteristic Impedances', IEEE Trans. Microw. Theory Tech., 
%                          1969. MT17(5): pp. 259-261
%                               
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cell_constant =Microchannel_cell_constants_for_coplanar_electrodes (electrode_width,gap_width,channel_height);

digits (100);   % for the variable precision, this value should be between 20-1000, 
                % the higher it is the slower the calculation becomes
                
%variable precision arithmetic value of Pi              
pi_long = vpa ( 3.141592653589793238462643383279502884197169399375105820974944592307816406286); 


%+++++++++++++++++++++++++++++++++++++++++++
%       Physical geometry 
%        (Z-plane)
%+++++++++++++++++++++++++++++++++++++++++++
za = -channel_height/2 + i * (gap_width/2 + electrode_width);
zb = -channel_height/2 + i * gap_width/2;
zc = -channel_height/2;
zd = channel_height/2;

%+++++++++++++++++++++++++++++++++++++++++++
%       SINE-TRANSFORMATION
%        (U-plane)
%+++++++++++++++++++++++++++++++++++++++++++
p=vpa(real(sin (pi*za/channel_height)));
q=vpa(real(sin (pi*zb/channel_height)));
r=vpa(real(sin (pi*zc/channel_height)));
s=vpa(real(sin (pi*zd/channel_height)));

%+++++++++++++++++++++++++++++++++++++++++++
%       BILINEAR TRANSFORMATION
%       the bilinear coordinates are not needed, we just need the value of
%       the modulus k (which is equal to 1/u4)
%+++++++++++++++++++++++++++++++++++++++++++

k=vpa((p*(q+r-2*s)+2*sqrt( (q-p)*(r-p))*sqrt(q-s)*sqrt(r-s) +r*s+q*(s-2*r)));
k1=vpa((q-r)*(p-s));

k=k/k1;
%+++++++++++++++++++++++++++++++++++++++++++
%       HILBERG APPROXIMATIONS
%       The elliptic functions are calculated using the
%       Hilberg-approximations, because the ellipke function 
%       is not defined in vpa (Variable Precision Arithmetics).
%       For k-values close to 1 or zero, the vpa yields higher resolution.
%+++++++++++++++++++++++++++++++++++++++++++

k1=vpa(sqrt(1-k*k));
bigk=double(k)>sqrt(1/2);

if (~bigk) % The formula used depends on the value of the modulus k.
    KK_Hilberg = 2 * pi_long / log (2*(sqrt(1+k1)+(4*k1)^0.25)/(sqrt(1+k1)-(4*k1)^0.25));
end;

if (bigk)
    KK_Hilberg =  log (2*(sqrt(1+k)+(4*k)^0.25)/(sqrt(1+k)-(4*k)^0.25))/(2 * pi_long);
end;
    
%+++++++++++++++++++++++++++++++++++++++++++
%       FINAL ANSWER
%+++++++++++++++++++++++++++++++++++++++++++

% The resistance, and therefore the cell_constant, is twice that found for the half-plane 
%   mapping, since only half the geometry was mapped in the sin-mapping

cell_constant = double(4*KK_Hilberg)