%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Filename:       'ElectricFieldInMicroChannel.m'
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Author:         Pontus Linderholm, LMIS4
%       Date:           06 March 2006  
%       Affiliation:    LMIS4, EPFL, Lausanne, Switzerland    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Description:    Calculates the electric field created by two coplanar
%       electrodes in a microchannel. 
%
%       NOTE that the algorithm assumes that the microchannel has a
%       rectangular cross-section (surrounded by insulating material). 
%
%       Geometry (upside down):
%
%
%                        ELECTRODE 1      GAP      ELECTRODE 2                          
%          \     \                                              \     \     \
%          xxxxxxxxxxxxx_______________xxxxxxxxxx_______________xxxxxxxxxxxxx   _           
%                                                                                |
%                                      (x,y) -->  Ex-component                   |             
%                                        |                                       |   MICROCHANNEL            
%                                        v                                       |             
%                                    Ey-component                                |
%                                                                               _|
%          xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%          /     /     /     /     /     /     /     /     /     /     /
%
%          NOTE: The origin is between the electrodes.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   
%   Theory: The function uses two conformal mappings to bring
%           a 2D geometry with two electrodes in a microchannel 
%           to a rectangle with the electrodes defining the sidewalls.
%           The electric field is calculated between one of the electrodes
%           and the symmetry plane along the y-axis.
%           The physical electrode is defined by z3-z4 and the symmetry
%           electrode z1-z2.
%           First a hyperbolic transformation is used to
%           "remove" the channel top wall (it is mapped onto
%           the real axis of the U-plane). The two electrodes are now
%           asymmetric. A bilinear transform is used to stretch the
%           electrodes to equal lengths. Lastly, the S-C transform is used
%           to transform the V-plane into a rectangle. This transformation
%           is not carried out in this function. The electric field in the
%           rectangle is known (V/length). 
%           Using a property of conformal mappings, the electric field in
%           the V- and U and finally Z-plane can be calculated from the
%           derivatives of the conformal transformations. 
%
%           If all this seems incomprehensible, I recommend you just try
%           to run the code to see what it gives.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       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
%            x                          the x-poisition
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       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) DEMIERRE et al.
%                               
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [Ex, Ey]= Electric_Field_in_Microchannel ( elec, gap, height, x, y, voltage);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              TRANSFORMING THE CORNER POINTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+++++++++++++++++++++++++++++++++++++++++++
%       Physical geometry 
%       Definition of the corner points
%        (Z-plane)
%+++++++++++++++++++++++++++++++++++++++++++
z1= i*height;
z2= 0;
z3= gap/2;
z4= gap/2 + elec;


%+++++++++++++++++++++++++++++++++++++++++++
%       SINH 
%       Corner points in U-plane
%+++++++++++++++++++++++++++++++++++++++++++
u1= sinh( pi/height * ( z1-i*height/2 ));
u2= sinh( pi/height * ( z2-i*height/2 ));
u3= sinh( pi/height * ( z3-i*height/2 ));
u4= sinh( pi/height * ( z4-i*height/2 ));


%+++++++++++++++++++++++++++++++++++++++++++
%       BILINEAR TRANSFORMATION I
%           V=(U+B)/(C*U+D)
%       Finding B,C,D values for correct transformation
%+++++++++++++++++++++++++++++++++++++++++++
tmp1= -u2*u2 +u1*(u2+u3-u4-u4) +2*sqrt((u1-u2)*(u1-u3)*(u2-u4)*(u3- u4)) +u2*u4 +u3*(u4-u3);
tmp2= (u2-u3)*(-u1+u2+u3-u4);

C= tmp1/tmp2;
C= real(C); % est purement réel, enlever tout résidu imaginaire

D= (u3*(1-C)-u2*(1+C))/2;
D= i*imag(D); % est purement imaginaire, enlever tout résidu réel

B= -u2*(1+C)-D; % est égalment purement imaginaire

%+++++++++++++++++++++++++++++++++++++++++++
%       BILINEAR TRANSFORMATION II
%       The interior electrode edges are at +/-1,
%        and the outer edges are given by 
%        v1 = -v4
%+++++++++++++++++++++++++++++++++++++++++++

v4= (u4+B)/(C*u4+D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              DEFINING THE MODULUS FOR THE S-C TRANSFORMATION (TO A W-PLANE NOT SHOWN)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k= 1/v4;
k= real(k); % to make sure that k is completely real


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              DEFINING THE TRANSFORMATIONS FOR THE POINT (x,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z=x+i*y;
u= sinh( pi/height*( z - i*height/2 ));
v= (u+B)./(C*u+D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              DEFINING THE DERIVATIVES OF THE TRANSFORMATIONS AT THE POINT (x,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dudz= pi/height * cosh(pi * z / height - i*pi/2);
dvdu= (D-B*C)./(D+C*u)./(D+C*u);
dwdv= 1./sqrt(1-k*k.*v.*v)./sqrt(1- v.*v);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              ELECTRIC FIELD CALCULATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ew= voltage / (4 * ellipke (k^2)); % In the W-plane (after S-C mapping) the electric field is 
                                   % given by the  


Ez= Ew.*conj(dwdv).*conj(dvdu).*conj(dudz);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       FINAL ANSWER!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ex= real(Ez);
Ey= imag(Ez);