### Calculates Cj from small signal simulation ###
####### and compares with analytic value #########

clear;

load('ac_pn_diode.ldev');
switchtolayout;
run('CHARGE');

# get simulation results

contact = 'emitter';
ac_data = getresult('CHARGE','ac_'+contact);
dI = pinch(getdata('CHARGE','ac_'+contact,'dI'));
dV = pinch(getdata('CHARGE','ac_'+contact,'dV_'+contact));
f = pinch(getdata('CHARGE','ac_'+contact,'f'));
Vdc = pinch(getdata('CHARGE','ac_'+contact,'V_'+contact));

# get the length of voltage and frequency sweeps

kern = size(dI);
nV = kern(1);
nf = kern(2);

f_matrix = matrix(nV,nf);
for (i=1:nV) {
f_matrix(i,1:nf) = f;  # frequency matrix
}

Ycomp = dI/dV;  # diode admittance

G = real(Ycomp);  # conductance
C = imag(Ycomp)/(2*pi*f_matrix);  # capaciatnce

# Compare to analytic junction capacitance

epsr = 11.7;
ND = 1e21;
NA = 1e23;
ni = 1.05e16;
Vbi = 0.025*log(NA*ND/ni^2);
A = 1e-6 * 1e-2;   # 0.2 um * 1 cm

Cj = A * sqrt((e*eps0*epsr*ND)/(2*(Vbi-Vdc)));  # unit is F

plot(Vdc,Cj,C(1:11,1),'V_emitter (V)','Capaciatnce (F)');  # plot junction capaciatnce
legend('C_analytic','C_simulation');

Cj_f = matrix(nf,1) + Cj(11);
plot(f,Cj_f,C(11,1:nf),'f (Hz)','Capaciatnce (F)','','log10x');  # plot junction capaciatnce with freq.
legend('C_analytic','C_simulation');