function [xx,yy,Volt,Cap,Vguess1] = ...
  SmartTrayFiniteDifference(L,W,nn,aspRat,aa,keff,Ds,Vguess,flG,flP,flF)

% This prgm solves the electrostatic problem of the Smart Tray using finite 
% differences, assuming a mesh of individual square cells. It then inverts
% the matrix directly to find the solution, or it uses the iterative
% Gauss-Seidel technique for large systems. If the latter, a guess of the
% solution can be input from a previous run of the program. Convergence
% of the Gauss-Seider technique is either achieved on the entire voltage
% field for best accuracy, or on the capacitance recorded by the ST sensor.
%f
% Michel Louge Nov 26, 2023
%
% for fastest handling of complex variables, see
% https://www.mathworks.com/matlabcentral/answers/...
%         372428-why-is-assigning-integers-to-a-complex-matrix-slow
%
%	input variables
%           -	L   	=	depth of powder (mm)
%           -   W       =   distance between outer singularities (mm)
%           -   nn      =   # of grid cells in powder height (-)
%           -   aspRat  =   Smart Tray aspect ratios = [a1 a2 a3 asm asp] 
%                           relative to W:
%                           a1  = width of left & right grounds (-)
%                           a2  = width of left and right guards (-)
%                           a3  = width left guard + inner ground (-)
%                           asm = left outer singul to left sensor pt (-)
%                           asp = left outer singul to right sensor pt (-)
%           -   aa      = air buffer above powder surface/powder depth (-)
%           -   keff    = input func or interpolant of complex diel cst (-)
%                           = f[downward dist rowvec from powder surf (mm)]                       
%           -   Ds      = sensor depth in dir perpendicular to (x,y) [mm]
%           -   Vguess  = initial colvec size(ii*jj,1) voltage solution
%                         guess for Gauss-Seidel technique; also serves as 
%                         a flag if its dimension is (1,1)
%                          if == scalar 0, direct matrix inversion
%                          if == scalar 1, Gauss-Seidel iteration
%           -   flG     =   if flG == 0, converge Gauss-Seidel with norm
%                           if flG == 1, converge GS with capacitance
%           -   flP     =   if flP == 1, print iter state (interact mode)
%                                        otherwise automatic GaussSeid mode
%           -   flF     =   if flF == 1, plot figure, otherwise no plot
%
%	output variables
%           -   Volt    =   solution voltge(x,y) or (i,j) array size(ii,jj)
%           -   xx      =   (jj,ii) meshgrid array of x-coord (mm) for plot
%           -   yy      =   (jj,ii) meshgrid array of y-coord (mm) for plot
%           -   Cap     =   complex capacitance (fF)
%           -   Vguess1 =   new dim(-) voltage guess, colvec size(ii*jj,1)
%
%   global variables
%           -   none    =   none
%
%   function called
%           -   none
%
% Stable examples:
%
% With matrix inversion (Ke=1 for air alone):
% [xx,yy,Volt,Cap,Vguess] =SmartTrayFiniteDifference(10.0,72.0,20,...
% [0.5312 0.3623 0.6377 0.7076 0.7492],2.0,@(x)ones(size(x)),31.5,0,1,1);
%
% Using Gauss-Seidel iteration converging to a stable capacitance with Ke=2
% [xx,yy,Volt,Cap,Vguess]=SmartTrayFiniteDifference(10.0,72.0,20,...
% [0.5312 0.3623 0.6377 0.7076 0.7492],2.0,@(x)2*ones(size(x)),31.5,1,1,1);
%
% Example with humid MCC at several RH:
%
% cc=0.1156;dd=0.7833; % coefficients of loss tangent vs RH at 35C
% aa=1.431;bb=3.761;ee=1.294; % coefficients of K'e dependence on RH
%
% RH=(0:0.1:0.8); % array of RH considered (no pure air yet) size 9
% KeReal=[1 ee+aa.*RH.^bb]; % arr real parts with pure air as 1st point
% losstan=[0 cc.*RH.^2+dd.*RH.^8]; % arr of loss tangs with air as 1st pt
% KeImag=KeReal.*losstan+1e-9; % array of imaginary parts 
% (ensures always complex with small imag part at RH=0 and in air)
% Cap=NaN(size(KeReal)); % init capacitance result arr size 10 (with air)
% for ii=1:size(KeReal,2) % loop: air, followed by MCC @ RH=0 and higher RH
% ke=KeReal(ii)+1i*KeImag(ii) % complex dielectric constant
% [~,~,~,Cap(ii),~] =SmartTrayFiniteDifference(11.0,72.0,28,...
%[0.5312 0.3623 0.6377 0.7076 0.7492],1.2,@(x)ke*ones(size(x)),31.5,0,0,1);
% end
% 
% KeRealApp = real(Cap(2:end))/real(Cap(1)) ; % array size 9 
% apparent real part relative to air = first element in Cap array
% KeImagApp = imag(Cap(2:end))/real(Cap(1)) ;
% losstanApp = KeImagApp./KeRealApp ;
%


%
% 	other internal variables
%
%           -   atheo   =   theoretical distance (mm)
%           -   a1      =   aspect ratio, see above
%           -   a2      =   aspect ratio, see above
%           -   a3      =   aspect ratio, see above
%           -   asm     =   aspect ratio, see above
%           -   asp     =   aspect ratio, see above
%           -   cc      =   colvec size ii*jj
%           -   CapDiff =   capacitance error
%           -   Cap_old =   previous capacitance guess
%           -   Ctheo   =   theoretical capacitance (fF)
%           -   dx      =   square grid cell width (mm)
%           -   eps0    =   pernittivity of free space (fF/mm)
%           -   flC     =   if == 1, complex variables throughout
%           -   fDM     =   finite diff matrix, size ii*jj x ii*jj (-)
%           -   g1theo  =   theoretical distance (mm)
%           -   g2theo  =   theoretical distance (mm)
%           -   gradX   =   x-gradient of dim(-) voltage magnitude (mm^-1)
%           -   gradY   =   y-gradient of dim(-) voltage magnitude (mm^-1)
%           -   i       =   index along downward distance from domain top
%           -	ii  	=	# grid cells in x-dir from top surf to base
%           -   iter    =   iteration number
%           -   j       =   index along rightward hor dist from domain left
%           -   jj      =   # grid cells in lateral y-dir
%           -   [j1 j2 j3 j4 j5 j6] base node geometry, where
%                                  1:j1     = left ground
%                               j1+1:j2     = left guard
%                               j2+1:j3     = inner ground
%                               j3+1:j4     = inner guard
%                               j4+1:j5     = sensor
%                               j5+1:j6     = right guard
%                               j6+1:jj     = right ground
%           -   ke      =   complex diel cst rowvec ke(i=1:ii)
%           -   kk      =   index of air-powder interface, 1 < kk < ii
%           -   l       =   equation index (or temperature index)
%           -   LGrad   =   gradient at left sensor-guard interf (mm^-1)
%           -   nmax    =   max # of iterations
%           -   normVal =   relative norm of error
%           -   RGrad   =   gradient at right sensor-guard interf (mm^-1)
%           -   sens    =   points in sensor line (mm)
%           -   s1theo  =   theoretical distance (mm)
%           -   s2theo  =   theoretical distance (mm)
%           -   theta1  =   polar angle on left ground-guard singularity
%           -   theta2  =   polar angle on left inner ground-guard singul
%           -   theta3  =   polar angle on right inner ground-guard singul
%           -   theta4  =   polar angle on right ground-guard singularity
%           -   tol     =   error tolerance between successive iterations
%           -   tolCap  =   capacitance tolerance between successive iter
%           -   V       =   voltage colvec solution, size(ii*jj,1)
%           -   Vmod    =   model relative voltage (-)
%           -   V_old   =   previous guess of voltage col vec
%           -   Vs      =   integrated voltage gradient over sensor to be
%                           ratioed with Vs in air to find the effective
%                           dielectric constant.
%           -   VsG     =   voltage difference over sensor nodes
%           -   x       =   dist rowvec x(i=1:ii) from domain top (mm)
%           -   xmod    =   model x-distance (mm)
%           -   ymod    =   model y-distance (mm)
%

a1   =   aspRat(1)   ;   % unpack aspect ratios
a2   =   aspRat(2)   ;
a3   =   aspRat(3)   ;
asm  =   aspRat(4)   ;
asp  =   aspRat(5)   ;

% theoretical capacitance (fF)

eps0    =   8.854           ;       % permittivity of free space (fF/mm)
atheo   =   W/2             ;    	% a in ref
g1theo  =   (a2-a3)*W/2     ;       % g1 < 0 in ref
g2theo  =   (a3-a2)*W/2     ;       % g2 > 0 in ref
s1theo  =   (asm-(a3+a2)/2)*W ;     % s1 in ref
s2theo  =   (asp-(a3+a2)/2)*W ;     % s2 in ref
Ctheo   =   eps0*Ds/pi*log( (atheo+s2theo)/...
                            (atheo+s1theo)*...
                            (s1theo-g1theo)/...
                            (s2theo-g1theo)*...
                            (s2theo-g2theo)/...
                            (s1theo-g2theo)*...
                            (atheo-s1theo)/...
                            (atheo-s2theo))     ;
                        
                        if L == 0 || L < 0
error('Powder height is zero. If air alone, just make ke=1.')
                        end

dx  =   L/nn                    ; % grid cell resolution (mm)

j1  =   1+floor(a1*W/dx)        ; % index left ground-guard interface
j2  =   j1+floor(a2*W/dx)       ; % index left guard-inner ground interface
j3  =   j2+floor((a3-a2)*W/dx)	; % index right inner ground-guard interf
j4  =   j3+floor((asm-a3)*W/dx)	; % index left guard-sensor interface
j5  =   j4+floor((asp-asm)*W/dx); % index right sensor-guard interface
j6  =   j5+floor((1-asp)*W/dx)	; % index right guard-ground interface
jj  =   j6+floor(a1*W/dx)       ; % index of right domain boundary

warning(['number of sensor discrete points = ' num2str(j5-j4-1)])

kk  =   1+floor(nn*aa)      	; % index of air-powder interface
ii  =   kk+nn                   ; % index of Smart Tray PCB base

% air indices:              1:kk-1
% powder indices:           kk+1:ii-1

% define the series of complex dielectric constants in the depth:

x   =   dx*(1:ii)                            ; % rowvec vert dist (mm)

ke  =   [ones(1,kk-1) keff(x(kk:ii)-x(kk))]  ; % rowvec complex diel cst(-)

flC =   isreal(ke)                ; % if flC==0 if ke has any complex value

fDM =   zeros(ii*jj,ii*jj)  ; % initialize finite difference matrix
cc  =   zeros(ii*jj,1)      ; % initialize eqn RHS i*j colvec

                % equation for left corner in air (l=1, first equation):
                if flC == 0 % complex variables
fDM(1,1)        =  3+ke(2)/ke(1)+1i*1e-9 	; % add a dummy complex value
                                              % to the first element
                                              % to speed up inspection
                                              % whether the array is
                                              % complex, see size above
                else        % real pb
fDM(1,1)        =  3+ke(2)/ke(1)            ;
                end
                
fDM(1,1+jj)     = -(1+ke(2)/ke(1))  ;
%
fDM(1,2)        = -2                ;
%
cc(1,1)         =  0                ;

for l=2:jj-1   	% eqn for upper zero-flux boundary cdt:
    fDM(l,l)    = 3+ke(2)/ke(1) 	;
    fDM(l,l+jj) = -(1+ke(2)/ke(1))  ;
%    
    fDM(l,l+1)  = -1                ;
    fDM(l,l-1)  = -1                ;
    cc(l,1)     =  0                ;
end

                % equation for right corner in air (l=jj):
fDM(jj,jj)    	=  3+ke(2)/ke(1) 	;
fDM(jj,jj+jj)   = -(1+ke(2)/ke(1))  ;
%
%
fDM(jj,jj-1)    = -2                ;
cc(jj,1)        =  0                ;

for i=(2:ii-1)  % eqn for left side domain:
    l=(i-1)*jj + 1                              ;
    fDM(l,l)    = 3+(ke(i-1)+ke(i+1))/ke(i)/2   ;
    fDM(l,l+jj) = -(1+ke(i+1)/ke(i))/2          ;
    fDM(l,l-jj) = -(1+ke(i-1)/ke(i))/2          ;
    fDM(l,l+1)  = -2                            ;
%
    cc(l,1)     =  0                            ;
end

for i=(2:ii-1) 	% eqn for right side domain:
    l=(i-1)*jj + jj                             ;
    fDM(l,l)    = 3+(ke(i-1)+ke(i+1))/ke(i)/2   ;
    fDM(l,l+jj) = -(1+ke(i+1)/ke(i))/2          ;
    fDM(l,l-jj) = -(1+ke(i-1)/ke(i))/2          ;
%
    fDM(l,l-1)  = -2                            ;
    cc(l,1)     =  0                            ;
end

for i=2:ii-1  	% eqn for entire domain (air+powder):
        for j=2:jj-1
        l = (i-1)*jj+j                          ; % equation index
    fDM(l,l)    =  3+(ke(i-1)+ke(i+1))/ke(i)/2	;
    fDM(l,l+jj) = -(1+ke(i+1)/ke(i))/2          ;
    fDM(l,l-jj) = -(1+ke(i-1)/ke(i))/2          ;
    fDM(l,l+1)  = -1                            ;
    fDM(l,l-1)  = -1                            ;
    cc(l,1)     =  0                            ;    
        end
end

for j=1:j1      % specify voltages on PCB left ground
% left ground indices:      1:j1
        l = (ii-1)*jj+j     ; % equation index
    fDM(l,l)    =  +1       ;
    cc(l,1)     =   0       ;    
end

for j=j1+1:j2 	% specify voltages on PCB left guard
% left guard indices:       j1+1:j2
        l = (ii-1)*jj+j     ; % equation index
    fDM(l,l)    =  +1       ;
    cc(l,1)     =  +1      	; % assume guard voltage = 1 (arbitrary unit)  
end

for j=j2+1:j3 	% specify voltages on PCB inner ground
% inner ground indices:     j2+1:j3
    l = (ii-1)*jj+j         ; % equation index
    fDM(l,l)    =  +1       ;
    cc(l,1)     =   0      	;  
end

for j=j3+1:j4 	% specify voltages on PCB inner guard
% inner guard indices:      j3+1:j4
    l = (ii-1)*jj+j         ; % equation index
    fDM(l,l)    =  +1       ;
    cc(l,1)     =  +1      	; % assume guard voltage = 1 (arbitrary unit)  
end

for j=j4+1:j5 	% specify voltages on PCB sensor
% sensor indices:           j4+1:j5
    l = (ii-1)*jj+j         ; % equation index
    fDM(l,l)    =  +1       ;
    cc(l,1)     =  +1      	; % assume sensor voltage = 1 (arbitrary unit)  
end

for j=j5+1:j6 	% specify voltages on PCB right guard
% right guard indices:      j5+1:j6
    l = (ii-1)*jj+j         ; % equation index
    fDM(l,l)    =  +1       ;
    cc(l,1)     =  +1      	; % assume guard voltage = 1 (arbitrary unit)  
end

for j=j6+1:jj 	% specify voltages on PCB right ground
% right ground indices:     j6+1:jj
    l = (ii-1)*jj+j         ; % equation index
    fDM(l,l)    =  +1       ;
    cc(l,1)     =   0      	;  
end

[xx,yy] = meshgrid(1:jj,1:ii)	;   % mesh grid for plots
xx      = xx*dx ; % convert to mm
yy      = yy*dx ; % convert to mm

% define sensor line

sens      = zeros(2,2)  ;
sens(1,1) = dx*(j4+1/2) ;
sens(1,2) = dx*ii       ;
sens(2,1) = dx*(j5+1/2) ;
sens(2,2) = dx*ii       ;

    if sum(size(Vguess))==2 && Vguess == 0 % seek direct matrix inversion

%
% Solve voltage field
% https://www.mathworks.com/help/matlab/ref/inv.html
%

V       = fDM\cc        ; % solution with standard Matlab matrix inversion

Vguess1 =   V   ; % record the solution as new guess for next run

% see also https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method
% Gauss-Seidel in Matlab: https://www.mathworks.com/matlabcentral/...
%                                 answers/518877-gauss-seidel-method-help
% p.1026 Incropera & Dewitt 7th-ed, Appendix D

% convert V to an (x,y) array for plotting and output:

Volt    = zeros(ii,jj)  ; % initialize voltage array

for l=1:ii*jj
    i = floor((l-1)/jj)+1 	; % find i,j indices
    j = l - jj*(i-1)        ;
    Volt(i,j) = V(l)        ; % fill in voltage array
end

% plot the voltage, its contours and field contours

[gradX,gradY]   =   gradient(abs(Volt))         ; % gradient field (mm^-1)
gradEx          =   sqrt(gradX.^2+gradY.^2)     ; % gradient magn field

RGrad = gradEx(size(gradEx,1),j5)    ; % grad @ right sensor-guard interf
LGrad = gradEx(size(gradEx,1),j4+1)  ; % grad @ left sensor-guard interface

if flF == 1
figure
surf(xx,yy,abs(Volt))    % surface plot; if complex, return magnitude

figure
contour(xx,yy,abs(Volt),'LineWidth',2,'Color','k') % voltage contour plot
hold on
contour(xx,yy,gradEx,[LGrad RGrad],'LineWidth',2,'Color','k') 
% contours of E-field from sensor
hold on
%plot(sens(1:2,1),sens(1:2,2),'-r','LineWidth',4);hold on  % sensor line
plot([xx(1,1) xx(end,end)],...
             [yy(kk,end) yy(kk,end)],':k','LineWidth',2)  % interface
set(gca,'ydir','reverse') % reverses the plot perpendicular to Smart Tray
axis equal
end

% calculate integrated voltage gradient over the sensor:

Vs  =   ( Volt(ii,j4)-Volt(ii-1,j4) )/2 + ...     % contributions of
        ( Volt(ii,j5)-Volt(ii-1,j5) )/2         ; % sensor extremities 
                                                    
for j = j4+1:j5-1

Vs = Vs + Volt(ii,j) - Volt(ii-1,j) ;

end

Cap    =   eps0*Vs*Ds*ke(ii)	; % complex capacitance (fF)

    end

% calculate the solution with Gauss-Seidel
% https://www.mathworks.com/matlabcentral/answers/...
%         518877-gauss-seidel-method-help

if sum(size(Vguess))==2 && Vguess > 0  % Gauss-Seidel technique
                                       % with exact solution as guess

% Calculate the exact sol for a guess with uniform ke and infinite domain:

Vmod    = zeros(ii,jj)      ; % initialize model relative voltage
Vguess  = zeros(ii*jj,1)    ; % init guess of solution colvec from model

for i = 1:ii
    for j = 1:jj
        
        xmod = (i-1/2)*dx ; % model x-distance
        ymod = (j-1/2)*dx ; % model y-distance
        
        theta1 = atan((ii*dx-xmod)/(ymod-a1*W)) + ...
                     ((ii*dx-xmod)/(a1*W-ymod)>0)*pi            ;
        theta2 = atan((ii*dx-xmod)/(ymod-(a1+a2)*W)) + ...
                     ((ii*dx-xmod)/((a1+a2)*W-ymod)>0)*pi       ;
        theta3 = atan((ii*dx-xmod)/(ymod-(a1+a3)*W)) + ...
                     ((ii*dx-xmod)/((a1+a3)*W-ymod)>0)*pi       ;
        theta4 = atan((ii*dx-xmod)/(ymod-(a1+1)*W)) + ...
                     ((ii*dx-xmod)/((a1+1)*W-ymod)>0)*pi        ;
                 
        Vmod(i,j) 	= (-theta1+theta2-theta3+theta4)/pi     ; % model
        
        Vguess(jj*(i-1)+j) = Vmod(i,j)  ; % solution guess
        
    end
end

end

if sum(size(Vguess))>2  % implements Gauss-Seidel technique

tol         =   1e-5            ; % Tolerance for method
tolCap      =   Ctheo*0.001     ; % capacitance tolerance in air (fF)
normVal     =   2*tol           ; % set error to number > tolerance
CapDiff     =   2*tolCap*abs(ke(ii)) ; % set error in capac diff to large #
nmax        =   1000            ; % max number of iterations

iter=0          ;

V   = Vguess    ;
Cap = 0         ;

while (flG==0)*normVal+(flG==1)*CapDiff > ...
      (flG==0)*tol+(flG==1)*tolCap*abs(ke(ii)) && iter<nmax
    
    V_old   =   V       ; % save earlier voltage matrix
    Cap_old =   Cap     ; % save earlier capacitance value
        
    for i=1:ii*jj
        
        guess=0;
        
        for j=1:i-1
            guess=guess+fDM(i,j)*V(j);
        end
        
        for j=i+1:ii*jj
            guess=guess+fDM(i,j)*V_old(j);
        end
        
        V(i)=(1/fDM(i,i))*(cc(i)-guess);
        
    end
    iter=iter+1 ;
    
% calculate the sensor capacitance

VsG =   ( V((ii-1)*jj+j4)-V((ii-2)*jj+j4) )/2 + ... % contributions of
        ( V((ii-1)*jj+j5)-V((ii-2)*jj+j5) )/2     ; % sensor extremities 
                                                    
for j = j4+1:j5-1

VsG     =   VsG + V((ii-1)*jj+j) - V((ii-2)*jj+j) ;

end

Cap         =   eps0*VsG*Ds*ke(ii)	; % complex capacitance (fF)

if flP == 1 % interactive mode, print iteration status, if flP = 1:
fprintf(['capacitance = ' num2str(Cap) ' fF with iteration ' ...
                                    num2str(iter) char(10)])     ;
end

CapDiff     =   abs(Cap_old-Cap)    ; % absol change of cap from last iter

normVal     =   norm(V_old-V)/ii/jj ; % norm of change from last iteration
    
end

Vguess1 =   V   ; % record the solution as new guess for next run

Volt   = zeros(ii,jj)  ; % initialize voltage array

for l=1:ii*jj
    i = floor((l-1)/jj)+1 	; % find i,j indices
    j = l - jj*(i-1)        ;
    Volt(i,j) = V(l)        ; % fill in voltage array
end

[gradX,gradY]   =   gradient(abs(Volt))         ; % gradient field (mm^-1)
gradEx          =   sqrt(gradX.^2+gradY.^2)     ; % gradient magn field

RGrad = gradEx(size(gradEx,1),j5)    ; % grad @ right sensor-guard interf
LGrad = gradEx(size(gradEx,1),j4+1)  ; % grad @ left sensor-guard interface

if flF==1
figure
surf(xx,yy,abs(Volt))    % surface plot; if complex, return magnitude

figure
contour(xx,yy,abs(Volt),'LineWidth',2,'Color','k') % voltage contour plot
hold on
contour(xx,yy,gradEx,[LGrad RGrad],'LineWidth',2,'Color','k') 
% contours of E-field from sensor
hold on
plot(sens(1:2,1),sens(1:2,2),'-r','LineWidth',2);hold on  % sensor line
plot([xx(1,1) xx(end,end)],...
             [yy(kk,end) yy(kk,end)],':k','LineWidth',2)  % interface
set(gca,'ydir','reverse') % reverses the plot perpendicular to Smart Tray
axis equal
end

end
