use_coarse_mesh = true;

#
# Set-up a voltage sweep with points V and V + dV
#

vmin = 0;
vmax = 1;
N = 6;
dV = 0.025;

vdv = matrix(2*N,1);
vdv(1:2:(2*N)) = linspace(vmin,vmax,N);
vdv(2:2:(2*N)) = vdv(1:2:(2*N)) + sign(vmax-vmin)*dV;


switchtolayout;
select('CHARGE::boundary conditions::cathode');
set('bc mode','steady state');
set('sweep type','value');
set('value table',vdv);

#
# Use coarse mesh to test calculation, or fine mesh for accurate CV calculation
#
if (use_coarse_mesh) {
  mesh_constraint = 5e-2*(1e-6);  # um
} else {
  mesh_constraint = 0.5e-2*(1e-6);  # um
}
setnamed('CHARGE::wg mesh','max edge length',mesh_constraint);   
setnamed('CHARGE::slab mesh','max edge length',mesh_constraint);

#
# Run the simulation
#
run;

#
# Analyze the results
#
rho_total = getresult('CHARGE::monitor','total_charge');
unit_len_um = getnamed('CHARGE simulation region','y span')*1e6;

n_total = pinch(rho_total.n);
p_total = pinch(rho_total.p);
Va = pinch(rho_total.getparameter('V_cathode'));
nva = length(Va);

dQn = n_total(2:2:nva) - n_total(1:2:nva);
dQp = p_total(2:2:nva) - p_total(1:2:nva);

#
# Calculate capacitance from both electron and hole densities
#
Cn = -e*dQn/dV;
Cp = -e*dQp/dV;

V = Va(1:2:nva);

#
# Plot result. If performing mesh converence tests, enable comparison plot 
# between Cn and Cp
#
#plot(V,Cn*1e15/unit_len_um,Cp*1e15/unit_len_um,"Voltage (V)","Capacitance (fF/um)");
plot(V,0.5*(Cn+Cp)*1e15/unit_len_um,"Voltage (V)","Capacitance (fF/um)");