#########################################################
#
# Simulate the SiGe LPD response under dark and illuminated
# conditions.
#
#(c) Lumerical Solutions 
#########################################################

clear;

project_name = 'lpd_sige';
sweep_name = 'compare_illum';
L = 15.3e-6;     # active photodetector length
compare_all = 0; # 1=compare all illumination cases, 0=dark and point/DBR only

#
# General setup
#
switchtolayout;
setnamed('CHARGE','norm length',L);
save;
orig_project = currentfilename;

slash = '\';  # or '/' for Linux
sweep_path = pwd+slash+project_name+'_'+sweep_name;
if(!fileexists(sweep_path)) { 
  system('mkdir "'+sweep_path+'"'); 
}

#
# Setup configurations for illumination (select OGR)
# -- first entry: -1=dark, 0=no DBR, 1=DBR
# -- second entry: -1=dark, 0=no contacts, 1=continous contacts, 2=point contacts
#
if (compare_all) {
  combos = [-1,-1; 0,0; 0,1; 1,1; 1,2];
} else { # just compare dark and best illum
  combos = [-1,-1; 1,2];
}
ncombos = size(combos); ncombos = ncombos(1);

#
# Load the job queue (simulate forward and reverse bias separately
#
clearjobs;
for (icombo = 1:ncombos) {
  ?'Building job for DBR '+num2str(combos(icombo,1))+', contact '+num2str(combos(icombo,2));

  if (combos(icombo,1) < 0) {
    file_tag = '_dark';
  } else {
    file_tag = '_'+num2str(combos(icombo,1))+num2str(combos(icombo,2));
  }

  switchtolayout;
  setnamed('::model','dbr',combos(icombo,1));
  setnamed('::model','contact',combos(icombo,2));

  #
  #reverse bias
  #
  select('CHARGE::boundary conditions::anode');
  set('range start',0);
  set('range stop',-3);
  set('range num points',16);

  fname = sweep_path+slash+project_name+file_tag+'_rvs.ldev';
  save(fname);
  addjob(fname);

  #
  # forward bias
  #
  select('CHARGE::boundary conditions::anode');
  set('range start',-0.1);
  set('range stop',1);
  set('range num points',23);

  fname = sweep_path+slash+project_name+file_tag+'_fwd.ldev';
  save(fname);
  addjob(fname);

}
load(orig_project);

#
# Run jobs
#
runjobs;

#
# Analyze results
#
Ires = cell(ncombos);
Vres = cell(ncombos);
for (icombo = 1:ncombos) {
  ?'Analyzing job for DBR '+num2str(combos(icombo,1))+', contact '+num2str(combos(icombo,2));

  if (combos(icombo,1) < 0) {
    file_tag = '_dark';
  } else {
    file_tag = '_'+num2str(combos(icombo,1))+num2str(combos(icombo,2));
  }

  fname = sweep_path+slash+project_name+file_tag+'_rvs.ldev';
  load(fname);
  anode = getresult('CHARGE','anode');
  Vr = anode.V_anode;
  Vr = Vr(2:length(Vr));
  Ir = pinch(anode.I);
  Ir = -Ir(2:length(Ir)) + 1e-10;

  fname = sweep_path+slash+project_name+file_tag+'_fwd.ldev';
  load(fname);
  anode = getresult('CHARGE','anode');
  Vf = anode.V_anode;
  If = pinch(anode.I) + 1e-10;

  Vres{icombo} = [flip(Vr,1); Vf];
  Ires{icombo} = [flip(Ir,1); If];
}

if (compare_all) {
  plotxy(Vres{1},abs(Ires{1}),
         Vres{2},abs(Ires{2}),
         Vres{3},abs(Ires{3}),
         Vres{4},abs(Ires{4}),
         Vres{5},abs(Ires{5}),
         'voltage (V)','current (A)');
  legend('dark','no contact, no DBR','continuous contact, no DBR',
         'continuous contact, w/ DBR','point contact, w/ DBR');
} else {
  plotxy(Vres{1},abs(Ires{1}),
         Vres{2},abs(Ires{2}),
         'voltage (V)','current (A)');
  legend('dark','point contact, w/ DBR');
}  
load(orig_project);