function [gres,Cres,chip,chim,reg,epsLp,epsLm,insigp,insigm,coa,cor] = ...
                            ElectroWetting(eps,alp,bst,coe,kap,xii,Ks,met)

% Implements the 'Statistical Mechanics of Electrowetting' of 
% Louge and Wang for a periodic square array of square pillars.
%
% Distances are made dimensionless with H = overall thickness of plate.
% Capacitances are dim(-) with [Ks usp0 (A0+As)]/H where
%                           Ks= dielectric cst of plate solid
%                           ups0 = 8.854fF/mm diel perm vacuum
%                           A0 = opening area of cavity
%                           As = flat solid area of unit cell
%
% Michel Louge May 15, 2019 MYL3@cornell.edu
%
%	input variables
%           -   eps     =   epsilon = surface area of cavity openings 
%           -	alp	    =	alpha = ratio of cavity opening to wall areas
%           -   bst     =   ratio pillar width to overall plate thickness H
%           -   coe     =   cosine of the Young contact angle
%           -   kap     =   dim(-) Laplace pressure parameter
%           -   xii     =   electrowetting param U^2 Ks ups0 /(2 H gam)
%                               where   U = voltage
%                                       gam = liquid-gas interfacial energ
%           -   Ks      =   dielectric constant of solid
%           -   met     =   if not 0, impose metastable regimes III and V
%
%	output variables
%           -   gres    =   rowvec geom parameters [eps alp lam bst dd hh]
%           -   Cres  	=   matrix of caps for cases 1-4
%           -   chip  	=   advancing transition probability
%           -   chim  	=   receding transition probability
%           -   reg     =   regime; I, II, III, IV, V, or VI are denoted by
%                                   1,  2,  3,  4,  5, and 6, respectively;
%                                       electrowetting-only regimes 
%                                       IIE, IIIE, and VE are 7, 8 and 9.
%           -   epsLp   =   epsilon*(advancing latent energy)
%           -   epsLm   =   epsilon*(receding latent energy)
%           -   insigp  =   integral of sigma bar d(chi) in advance
%           -   insigm  =   integral of sigma bar d(chi) in recession
%           -   coa     =   advancing contact angle cosine
%           -   cor     =   receding contact angle cosine
%
%   global variables
%           -   none  	=   none
%
% Stable example:
% [gres,Cres,chip,chim,reg,epsLp,epsLm,insigp,insigm,coa,cor] = ...
%                            ElectroWetting(.5,2,.3,cos(pi/3),0,.1,3.2,0);
%

%
% 	other internal variables
%
%           -   lam     =   ratio cross-section links to nearest neighbors
%                               and cavity opening
%           -   hh      =   ratio cavity height and 
%                               plate overall thickness H
%           -   dd      =   ratio gap between pillars and H
%           -   p1,..,p5=   five parameters best-fitted to simulations
%           -   nn      =   number of unlike near neighbors
%           -   Cd      =   capacitance from cavity bottom to ground
%           -   Cp      =   pillar capacitance below the base
%           -   CqI     =   inverse of pillar capacitance above the base
%           -   Cw      =   pillar to cavity coupling capacitance
%           -   Cl      =   pillar to next neighbor cell coupling capacit
%           -   C20     =   case 2 dim(-) cap with all like near-neighbors
%           -   Zab     =   intermediate impedance in case 2 
%           -   Zbc     =   intermediate impedance in case 2 
%           -   ZcaI    =   inverse of intermediate impedance in case 2 
%           -   Za      =   intermediate impedance in case 2 
%           -   Zb      =   intermediate impedance in case 2 
%           -   Zc      =   intermediate impedance in case 2 
%           -   Zl      =   intermediate impedance in case 2 
%           -   Zr      =   intermediate impedance in case 2 
%           -   Z2      =   intermediate impedance in case 2 
%           -   C24     =   case 2 dim(-) cap with nn unlike near-neighbors
%           -   C10     =   case 1 dim(-) cap with all like near-neighbors
%           -   ZfI     =   inverse of intermediate impedance in case 1
%           -   Z1I     =   inverse of intermediate impedance in case 1
%           -   C14     =   case 1 dim(-) cap with nn unlike near-neighbors
%           -   C30     =   case 3 dim(-) cap with all like near-neigh (=0)
%           -   Cm      =   cap frm cav base to crss-sec unlike neigh case3
%           -   CmCd    =   ratio Cm/Cd
%           -   C34     =   case 3 dim(-) cap with nn unlike near-neighbors
%           -   CdCc    =   ratio Cd/Cc
%           -   C40     =   case 4 dim(-) cap with all like near-neighbors
%           -   CkCc    =   ratio Ck/Cc
%           -   C44     =   case 4 dim(-) cap with nn unlike near-neighbors
%           -   deltac  =   (C20-C10)/C40
%           -   C14eps  =   expanded expression to o(eps) used in chip
%           -   C40eps  =   expanded expression to o(eps) used in chip
%           -   C10eps  =   expanded expression to o(eps) used in chim
%           -   C44eps  =   expanded expression to o(eps) used in chim
%           -   CmmC    =   Cbar(sig=-1;sigbar=-1;chi)
%           -   CmpC    =   Cbar(sig=-1;sigbar=+1;chi)
%           -   CpmC    =   Cbar(sig=+1;sigbar=-1;chi)
%           -   CppC    =   Cbar(sig=+1;sigbar=+1;chi)
%

%%
% Warn and stop if input parameters are out of bounds
%
if eps == 1 || eps > 1 || eps < 0
    error(['epsilon = ',num2str(eps),' must be >= 0 and < 1'])
end
%
if alp < 1 || alp == 1
    error(['alpha = ',num2str(alp),' must be > 1'])
end
%
    if eps ~= 0
if bst < 0 || bst == 0 || bst > 4*(1-eps)/eps/(alp-1)
    error(['bst =',num2str(bst),' > ',...
                num2str(4*(1-eps)/eps/(alp-1)),' is out of bounds'])
% in general, substrate thickness below pillars must be > 0; then, h/H > 1,
% which implies bst > 4*(1-eps)/eps/(alp-1) according to Eq. (13)
end
    else
if bst < 0 || bst == 0
    error(['bst =',num2str(bst),' out of bounds'])
end
    end
%
if coe < -1 || coe > +1
    error(['cos(thetae) =',num2str(coe),' out of bounds'])
end
%
if xii < 0
    error(['xi =',num2str(xii),' out of bounds'])
end
%
%%
% In the geometry of periodically-spaced square pillars, transform the
% input variables epsilon, alpha and b to other dimensions of texture:
%
lam = (alp-1)*(1/sqrt(1-eps)-1)     ;
hh  = (alp-1)*bst*eps/4/(1-eps)     ;
dd  = bst*lam/(alp-1)               ;
%
gres= [eps alp lam bst dd hh] ; % build rowvec of geometrical parameters

%%
% Calculation of lumped-parameter capacitances
%
% best-fit lumped-parameters specifications:
%
p1  =   0.929   ;
p2  =   0.197   ;
p3  =   0.044   ;
p4  =   0.065   ;
p5  =   0.100   ;

nn  =   4       ; % theory only needs four unlike near-neighbors

Cd = ((bst+dd)^2-bst^2)/(1-hh)/(bst+dd)^2   ;
Cp = bst^2/(1-hh)/(bst+dd)^2                ;
CqI= hh*(bst+dd)^2/bst^2                    ; % calculate inverse of Cq
                                              % to avoid blowup with h=0
Cw = 4*hh/p2/(bst+dd)^2                     ;
Cl = hh/p3/(bst+dd)^2                       ;
%
C20 = Cd+1/(CqI/(1/p1+CqI*Cw)+(1-p1)*CqI+1/Cp)  ;
%
Zab = CqI*p1/nn             ;
Zbc = CqI/(1/p1+Cw*CqI)     ;
ZcaI= nn*Cl                 ; % calc inverse of Zca to avoid h=0 blowup
%
Za  = Zab/(1+Zab*ZcaI+Zbc*ZcaI)             ;
Zb  = Zab*Zbc*ZcaI/(1+Zab*ZcaI+Zbc*ZcaI)  	;
Zc  = Zbc/(1+Zab*ZcaI+Zbc*ZcaI)             ;
%
Zl  = 1/Cp+(1-p1)*CqI   ;
Zr  = Zl/nn             ;
%
Z2  = Zl+Zc+Zb*(1+(Zl+Zc)/(Zr+Za))  ;
C24 = Cd+1/Z2                       ;
%
C10 = Cd+Cw/(1+(1-p1)*Cw*CqI+Cw/Cp) ;
%
Zg  = 1/Cp+(1-p1)*CqI               ;
ZfI = nn*Cl/(1+(1-p1)*Cl*CqI+Cl/Cp) ;
Z1I = Cw/(1+Zg*ZfI+Zg*Cw)           ;
C14 = Cd+Z1I                        ;
%
C30 = 0     ; % always zero
%
Cm  = hh/p4/Ks/(bst+dd)^2           ;
CmCd= (alp-1)*(1-hh)/4/p4/Ks/bst    ; % ratio Cm/Cd is finite
C34 = Cm/(CmCd+1/nn)                ;
%
CdCc= Ks*hh/(1-hh)                  ;
C40 = Cd/(1+CdCc)+Cp/(1+Cp*CqI)     ;
%
CkCc= (alp-1)^2*eps/16/(1-eps)/p5               ;
C44 = Cp/(1+Cp*CqI)+Cd/(1+CdCc/(1+nn*CkCc))     ;
%
% build matrix of capacitances for homogeneous (first row) and
% four unlike near-neighbors (second row):
%
Cres= [C10 C20 C30 C40;C14 C24 C34 C44];
%
deltac  =   (C20-C10)/C40 ; 
%
%%
% Calculation of transition probabilities
%
    if eps > 0.000001 % for sufficiently large eps, no singularity
%
chip=(1-alp*coe+lam+kap-xii*(C14-C30)/eps)/(2-xii*(C40-C24+C14-C30)/eps);
chim=(1-alp*coe-lam+kap-xii*(C10-C34)/eps)/(2-xii*(C44-C20+C10-C34)/eps);
%
    else % if eps < 10^(-6), expand (C14-C30)/eps, (C10-C34)/eps,
         % (C40-C24+C14-C30)/eps and (C44-C20+C10-C34)/eps to remove the
         % singularity in expressions of chip and chim:
%
C14eps = 1+(alp-1)/bst/p2+...
         eps/4*(bst*(alp-1)-4*(p2+p3)*(alp-1)^2/bst^2/p2^2/p3)          ;
C40eps = 1+(alp-1)/bst/p2-...
         eps/4*(bst*(alp-1)*(Ks-1)+4*(p2+p3)*(alp-1)^2/bst^2/p2^2/p3)   ;
C10eps = (alp-1)/bst/p2+bst*Ks*p4/(alp-1+bst*Ks*p4)+...
         eps/4*(bst*(alp-1)-4*(alp-1)^2/bst^2/p2^2-...
                                    bst*(alp-1)^3/(alp-1+bst*Ks*p4)^2)  ;
C44eps = (alp-1)/bst/p2+bst*Ks*p4/(alp-1+bst*Ks*p4)-...
         eps/4*(bst*(alp-1)*(Ks-1)+4*(alp-1)^2/bst^2/p2^2+...
                                    bst*(alp-1)^3/(alp-1+bst*Ks*p4)^2)  ;
%
chip=(1-alp*coe+lam+kap-xii*C14eps)/(2-xii*C40eps);
chim=(1-alp*coe-lam+kap-xii*C10eps)/(2-xii*C44eps);
%
    end
%
%%
% Regime identification
%
% Regimes are mutually exclusive.
%
reg = 1*... % regime I
        ((chip>0)+(chip==0))*((chip<1)+(chip==1))*... % 0<=chip<=1
        ((chim>0)+(chim==0))*((chim<1)+(chim==1))+... % 0<=chim<=1
      2*... % regime II
	    ((chip>0)+(chip==0))*((chip<1)+(chip==1))*... % 0<=chip<=1
         (chim<0)                                +... % chim < 0
      3*... % regime III
         (chip>1)                                *... % chip > 1
        ((chim>0)+(chim==0))*((chim<1)+(chim==1))+... % 0<=chim<=1
      4*... % regime IV
         (chip<0)                                *... % chip < 0
         (chim<0)                                +... % chim < 0
      5*... % regime V
         (chip>1)                                *... % chip > 1
         (chim<0)                                +... % chim < 0
      6*... % regime VI
         (chip>1)                                *... % chip > 1
         (chim>1)                                +... % chim > 1
      7*... % regime IIE (electrowetting only)
	    ((chim>0)+(chim==0))*((chim<1)+(chim==1))*... % 0<=chim<=1
         (chip<0)                                +... % chip < 0
      8*... % regime IIIE (electrowetting only)
         (chim>1)                                *... % chim > 1
        ((chip>0)+(chip==0))*((chip<1)+(chip==1))+... ;% 0<=chip<=1
      9*... % regime VE (electrowetting only)
         (chim>1)                                *... % chim > 1
         (chip<0)                                   ; % chip < 0      
%
%%
% Latent energy calculation
%
epsLp = 0 ; % advancing latent energy is zero by default
epsLm = 0 ; % receding latent energy is zero by default
%
if reg == 1 || reg == 2 || reg == 8 % advance regimes I, II and IIIE
%
    CmmC  = chip*C20+(1-chip)*C10 ; % Cbar(sig=-1;sigbar=-1;chi)
    CmpC  = chip*C24+(1-chip)*C14 ; % Cbar(sig=-1;sigbar=+1;chi)
    CpmC  = chip*C44+(1-chip)*C34 ; % Cbar(sig=+1;sigbar=-1;chi)
    CppC  = chip*C40+(1-chip)*C30 ; % Cbar(sig=+1;sigbar=+1;chi)
%
    epsLp = -eps*(2*chip-1+alp*coe-kap) -...
            xii/2*(CmmC-CpmC-CppC+CmpC)         ; % eps*latent energy
end
%
if reg == 1 || reg == 7 % recession regime I or IIE
%
    CmmC  = chim*C20+(1-chim)*C10 ; % Cbar(sig=-1;sigbar=-1;chi)
    CmpC  = chim*C24+(1-chim)*C14 ; % Cbar(sig=-1;sigbar=+1;chi)
    CpmC  = chim*C44+(1-chim)*C34 ; % Cbar(sig=+1;sigbar=-1;chi)
    CppC  = chim*C40+(1-chim)*C30 ; % Cbar(sig=+1;sigbar=+1;chi)
%
    epsLm = +eps*(2*chim-1+alp*coe-kap) +...
            xii/2*(CmmC-CpmC-CppC+CmpC)         ; % eps*latent energy
end
%
if reg == 3 % recession regime III
%
    if met == 0 % if regime III is not metastable
%
    CmmC  = chim*C20+(1-chim)*C10 ; % Cbar(sig=-1;sigbar=-1;chi)
    CmpC  = chim*C24+(1-chim)*C14 ; % Cbar(sig=-1;sigbar=+1;chi)
    CpmC  = chim*C44+(1-chim)*C34 ; % Cbar(sig=+1;sigbar=-1;chi)
    CppC  = chim*C40+(1-chim)*C30 ; % Cbar(sig=+1;sigbar=+1;chi)
%
    epsLm = +eps*(2*chim-1+alp*coe-kap) +...
            xii/2*(CmmC-CpmC-CppC+CmpC)         ; % eps*latent energy
    else
    epsLm = 0 ; % if regime III is metastable, then latent energy vanishes
    end
%
end
%
%%
% Integral sum of sigma bar d(chi)
%
insigp = 0 ; % initialize advance integral
insigm = 0 ; % initialize recession integral
%
if reg == 1 % regime I
    insigp = 2*chip-1 ; % advance
    insigm = 2*chim-1 ; % recession
end
%
if reg == 2 % regime II
    insigp = 2*chip-1 ; % advance
    insigm =       -1 ; % recession
end
%
if reg == 3 % regime III
%
    if met == 0 % if regime III is not metastable
%
    insigp =       +1 ; % advance
    insigm = 2*chim-1 ; % recession
%
    else % if reg III metastable, advance and recession have same integral
%
    insigp =       +1 ; % advance
    insigm =       +1 ; % recession
%
    end
%
end
%
if reg == 4 % regime IV
    insigp =       -1 ; % advance
    insigm =       -1 ; % recession
end
%
if reg == 5 % regime V
%
    if met == 0 % if regime V is not metastable
%
    insigp =       +1 ; % advance
    insigm =       -1 ; % recession
%
    else % if regime V metastable, advance and recession have same integral
%
    insigp =       +1 ; % advance
    insigm =       +1 ; % recession
%
    end
%
end
%
if reg == 6 % regime VI
    insigp =       +1 ; % advance
    insigm =       +1 ; % recession
end
%
if reg == 7 % regime IIE (electrowetting only)
    insigp =       -1 ; % advance
    insigm = 2*chim-1 ; % recession
end
%
if reg == 8 % regime IIIE (electrowetting only)
    insigp = 2*chip-1 ; % advance
    insigm =       +1 ; % recession
end
%
if reg == 9 % regime VE (electrowetting only)
    insigp =       -1 ; % advance
    insigm =       +1 ; % recession
end
%
%%
% Find contact angles
%
coa = (1-eps)*coe-eps*insigp+... % advance
      xii*C40*(1+deltac)/2*(1+(1-deltac)/(1+deltac)*insigp)+epsLp ;
%
cor = (1-eps)*coe-eps*insigm+... % recession
      xii*C40*(1+deltac)/2*(1+(1-deltac)/(1+deltac)*insigm)-epsLm ;
%
%%
