####################################################################
#
# Analyze the results from the transient simulation of the CIS pixel
#
#
#
# (c)2012 Lumerical Solutions, Inc.
####################################################################

Ccg = 1e-15; # conversion gain capacitance from sense contact
pwr = [10;linspace(20,100,5)]; # illumination intensity (W/m2) from sweep

dark_current = 1e5;  # e/s
OEb = 0.1734/0.25;   # Blue OE, scaled to max (25%)

Texp = 2e-4;  # illumination time (s)

#
# Load the sweep results
#
origname = currentfilename;
if (1) {
  count = matrix(6,1);

  t = matrix(2000,6);
  Vs = matrix(2000,6);
 
  cd("cis-pulse-1fF-200us_illumination");
  for (j = 1:6) {
    fname = "illumination_" + num2str(j) + ".ldev";
    load(fname);
  
    tj = pinch(getdata("CHARGE","sense.t"));
    st = size(tj); 
    count(j) = st(1);
  
    if (count(j) > 2000) { 
      message("Only 2000 entries in " + num2str(j) + " are recorded");
      count(j) = 2000;
    }

    t(1:count(j),j) = tj(1:count(j));
  
    Vsj = pinch(getdata("CHARGE","sense.Vc"));
    Vs(1:count(j),j) = Vsj(1:count(j));

  }
  load(origname);
}

#
# Interpolate the transient response onto the same time points
#
tmin = 0;
tmax = max(t);
ni = max(count);

ti = linspace(tmin,tmax,ni);
Vsi = matrix(ni,6);
for (j = 1:6) {
  Vsi(1:ni,j) = interp(Vs(1:count(j),j),t(1:count(j),j),ti);
}
plot(ti*1e6,Vsi(1:ni,1),Vsi(1:ni,2),Vsi(1:ni,3),Vsi(1:ni,4),Vsi(1:ni,5),Vsi(1:ni,6),"Time (us)","Vsense (V)");

ioff = find(ti >= Texp);
ioff = ioff(1);

Tsim = tmax - tmin;

#
# Calculate integrated charge
#
qQ = Ccg*(3.3 - Vsi(ioff,1:6))/1.602e-19;
qQ2 = Ccg*(3.3 - Vsi(ni,1:6))/1.602e-19;

plot(pwr,qQ/1000,qQ2/1000,"Illumination Power (W/m2)","Charge count (ke-)");
plot(pwr,(qQ2-qQ)/1000,"Illumination Power (W/m2)","Full-well capacity (ke-)");

#
# Calculate SNR
#
dark_count = dark_current*Texp; 
Nsig = pwr*OEb*(2.07e-6)^2*(.46e-6/6.63e-34/3e8)*Tsim; 

dark_noise_dB = matrix(6,1);
dark_noise_dB(1:6) = 10*log10(dark_count);

snr = 20*log10(qQ2) - 10*log10(dark_count) - 10*log10(Nsig);
plot(log10(pwr),snr,"Illumination Power (W/m2, log scale)","SNR (dB)");