######################################################
#
# Run the R and C calculations for the bandwidth of the
# nanobeam PC electro-optic modulator
#
# Inputs 
# - calculate_R: also perform the calculation of the resistance
# - Nhole: number of holes in the grating excluding taper
# - pc_length: length (m) of the 1D PC
# - taper_endpoint: the position (m) of the end of the taper
#     region in the 'y' direction (includes first six holes
#     from the middle).
#
# (c) 2015 Lumerical Solutions, Inc.
######################################################

calculate_R = 1;
Nhole = 24;
pc_length = 14.8e-6;
taper_endpoint = 2.251e-6;

#
# Set up the {V,V+dV} voltage sweep
#
N = 6;     dV = 0.025;
vmin = 0;  vmax = 0.5;
vdv = matrix(2*N,1);
Va = linspace(vmin,vmax,N);
vdv(1:2:(2*N)) = Va;
vdv(2:2:(2*N)) = Va + dV;

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


#
# Run the sweep for the taper region
#
setnamed('CHARGE simulation region', 'dimension', '3D');
setnamed('CHARGE simulation region', 'y min', 0);
setnamed('CHARGE simulation region', 'y max', taper_endpoint);
run('CHARGE');

#
# Get the net charge and calculate dQ/dV for the taper/cavity
#
net_charge = getresult('CHARGE::cv_monitor','net_charge');
Q = pinch(net_charge.Q);
dQ = Q(2:2:(2*N)) - Q(1:2:(2*N));
Ccav = 2*dQ/dV;

#
# Move the simulation region and run for a single hole
#
switchtolayout;
setnamed('CHARGE simulation region','y min',taper_endpoint+3*0.423e-6);
setnamed('CHARGE simulation region','y max',taper_endpoint+4*0.423e-6);
run('CHARGE');

#
# Get the net charge for a single hole, and calculate dQ/dV
#
net_charge = getresult('::model::CHARGE::cv_monitor','net_charge');
Q = pinch(net_charge.Q);
dQ = Q(2:2:(2*N)) - Q(1:2:(2*N));
Choles = Nhole*dQ/dV;

#
# If enabled, calculate Rslab by moving the simulation region
# to the test cross-section, and using the test contact
#
if (calculate_R) {
  switchtolayout;
  setnamed('CHARGE simulation region', 'dimension',2); # y-normal
  setnamed('CHARGE simulation region', 'y', 8.5e-6);
  run;

  test_contact = getresult('CHARGE','anode_test');
  I_test = pinch(test_contact.I);
  R_test = 0.1/I_test(2);          # ohm-cm
  Rslab = R_test/(pc_length*100);  # 0.04 ohm-cm
}

# 
# Calculate the 3dB bandwidth
#
f3dB = 1/(2*pi*2*Rslab*(Choles+Ccav));

?"Bandwidth calculation results:";
?"  R = " + num2str(R_test) + " ohm-cm";
?"  C = " + num2str(mean(Choles+Ccav)*1e15) + " fF";
?"  f3dB = " + num2str(mean(f3dB)/1e9) + " GHz";

#plot(Va,f3dB/1e9,"Voltage (V)", "Bandwidth (GHz)");