clear;

load('ssac_moscap.ldev');
switchtolayout;
error_count = 0;

#analytic

eps_ox = 3.9;
eps_si = 11.7;
d_ox = 1e-6; # cm
NA = 1e16;  #  /cm^3
ni = 1.5e10; #  /cm^3
kT_q = 0.0259;  # eV
VFB = 0.675;
Winv = 0.002e-4; # cm
Wx = getnamed('CHARGE simulation region','x span');
Vdc = linspace(-2.5,2.5,21);


Cox = eps_ox*eps0*1e-2/d_ox; # F/cm^2
phi_F = kT_q*log(NA/ni);
Wm = 2*sqrt(eps_si*eps0*1e-2*phi_F/e/NA);  # cm
Qd = -e*NA*Wm;
VT = -Qd/Cox + 2*phi_F;

Cd = eps_si*eps0*1e-2/Wm;
Cmin = Cox*Cd/(Cox+Cd);

EOT = d_ox + (eps_ox/eps_si)*Winv;
Cinv = eps_ox*eps0*1e-2/EOT; # F/cm^2


# SSAC

switchtolayout;
select("CHARGE::boundary conditions::gate");
set('bc mode','steady state');
set('sweep type','range');
set('range start',-2.5);
set('range stop',2.5);
set('range num points',21);
set('apply ac small signal','all');
setnamed("CHARGE","solver mode","ssac");
setnamed("CHARGE","fermi statistics",1);
setnamed("CHARGE::mesh","max edge length",1e-8);

run("CHARGE");

ac_gate = getresult("CHARGE","ac_gate");
Yd = pinch(imag(ac_gate.dI))/getnamed("CHARGE","perturbation amplitude");
w = 2*pi*ac_gate.f;
Va = linspace(-2.5,2.5,21);
w2 = meshgridy(Va,w);
Cac = Yd/w2/(Wx*1e2);  # F/cm^2

matlabload('moscap_static_Cap');

## plot results

plot(Vdc,Cac(1:21,1)*1e9);
holdon;
plot(Vdc,Cac(1:21,10)*1e9);
plot(Vdc,C_static*1e9);
plot(Vdc,matrix(21,1)+Cox*1e9);
plot(Vdc,matrix(21,1)+Cmin*1e9);

legend('C_ac_low_f','C_ac_high_f','C_static','Cox','Cmin');

C_ac_LF = Cac(1:21,1);  # F/cm^2
C_ac_HF = Cac(1:21,10);  # F/cm^2