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

% Implements the 'Statistical Mechanics of Electrowetting' of 
% Louge and Wang for a periodic square array of cylindrical pillars on a
% textured surface lined with a dielectric layer of uniform thickness H
% insolating the conductive textured solid from its exposed surface, as
% done in Krupenkin, T.N., Taylor, J.A., Schneider, T.M. and Yang, S., 2004
% From rolling ball to complete wetting: the dynamic tuning of liquids on 
% nanostructured surfaces. Langmuir, 20(10), pp.3824-3827. 
% https://doi.org/10.1021/la036093q
%
% Distances are made dimensionless with H
% 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 January 18, 2024 MYL3@cornell.edu
%
%	input variables
%           -   eps     =   epsilon = surface area of cavity openings 
%           -	alp	    =	alpha = ratio of cavity opening to wall areas
%           -   lam     =   ratio cross-section links to nearest neighbors
%                               and cavity opening
%           -   bst     =   ratio pillar diameter to dielectric 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 insulator
%           -   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] = ...
%                 ElectroWettingKrupenkin(.5,3,2,2.1,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 insulator thickness H
%           -   dd      =   ratio gap between pillars and H
%           -   p4,p5   =   two parameters best-fitted to simulations
%           -   nn      =   number of unlike near neighbors
%           -   Cd      =   capacitance from cavity surface to ground
%           -   Cp      =   pillar capacitance below the base
%           -   C20     =   case 2 dim(-) cap with all like near-neighbors
%           -   C24     =   case 2 dim(-) cap with nn unlike near-neighbors
%           -   C10     =   case 1 dim(-) cap with all like near-neighbors
%           -   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
%           -   Ck      =   cap frm cav base to crss-sec unlike neigh case4
%           -   C34     =   case 3 dim(-) cap with nn unlike near-neighbors
%           -   C40     =   case 4 dim(-) cap with all like near-neighbors
%           -   C44     =   case 4 dim(-) cap with nn unlike near-neighbors
%           -   deltac  =   (C20-C10)/C40 (here always equal to 1)
%           -   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 < 1-pi/4
    error(['In this geometry, epsilon = ',...
            num2str(eps),...
            ' must be >= pi/4 and < 1'])
end
%
if alp < 1 || alp == 1
    error(['alpha = ',num2str(alp),' must be > 1'])
end
%
if bst < 2
    error(['bst =',num2str(bst),...
        ' is smaller than twice the dielectric substrate thickness '])
% in general, substrate thickness must be > 2 pillar diameter
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 a square array of cylindrical pillars, transform the
% input variables epsilon, alpha and b to other dimensions of texture:
%
%lam = 4/pi*(alp-1)*(sqrt(pi/4/(1-eps))-1)   ;
hh  = (alp-1)*bst*eps/4/(1-eps)             ;
dd  = bst*lam/(alp-1)*pi/4                  ;
%
gres= [eps alp lam bst dd hh] ; % build rowvec of geometrical parameters

%%
% Calculation of lumped-parameter capacitances
%
% best-fit lumped-parameters specifications:
%
p4  =   0.065   ;
p5  =   0.100   ;

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

Cd = eps*alp    ;
Cp = 1-eps    	;
%
C20 = Cd+Cp     ;
%
C24 = C20     	;
%
C10 = Cd        ;
%
C14 = Cd        ;
%
C30 = 0         ; % always zero
%
Cm  = (alp-1)*eps/pi/bst/p4/Ks  ;
C34 = Cd*Cm/(Cm+Cd/nn)         	;
%
C40 = Cp        ;
%
Ck  = (alp-1)*eps/pi/bst/p5/Ks   	;
C44 = Cd*Ck/(Ck+Cd/nn)           	;
%
% 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 ; % deltac should be 1 by def
%
%%
% Calculation of transition probabilities
%
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);
%
%%
% 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 ;
%
%%
