#############################################################################
# Calculate IV curves and responsivity for the application example for 
# avalanche photodiode (APD) from Huang et al., Optica, Vol. 3, No. 8, 
# Aug. 2016, p. 793.
#
# Steps:
# 1. Load and plot dark current Idark. 
# 2. Load and plot photocurrent Iphoto.
# 3. Calculate responsivity as (Iphoto-Idark)/Pin, where Pin is input optical 
#    power.
# 4. Find unity gain Vmg1 where the second derivative d^2Iphoto/dV^2 is zero.
# 5. Find gain at other bias points V as Iphoto(V)/Iphoto(Vmg1).
##############################################################################


clear;
L = 10; # um  (choose either 10 or 50)
biased_contact = 'collector';
plot_contact = 'collector';

########################## 1. Load and plot dark current #######################

load('apd_electrical_'+num2str(L)+'um_dark.ldev');
#switchtolayout;
#run('CHARGE');

emitter = getresult('CHARGE',plot_contact);
Itotal = pinch(emitter.I);
eval('Vtotal = pinch(emitter.V_'+biased_contact+');');

Idark_wcl = -2*Itotal;
Vdark_wcl = Vtotal;

N = length(Vdark_wcl);
plot(Vdark_wcl(2:N),abs(Idark_wcl(2:N)),"Bias [V]","Current [A]",num2str(L)+"um device","log10y,plot type=line, color=blue, pen=-, linewidth=2");
holdon;

########################## 2. Load and plot photocurrent #######################

load('apd_electrical_'+num2str(L)+'um_photo.ldev');
#switchtolayout;
#run('CHARGE');

emitter = getresult('CHARGE',plot_contact);
Itotal = pinch(emitter.I);
eval('Vtotal = pinch(emitter.V_'+biased_contact+');');

Iphoto_wcl = -2*Itotal;
Vphoto_wcl = Vtotal;

plot(Vphoto_wcl,abs(Iphoto_wcl),"Bias [V]","Current [A]",num2str(L)+"um device","log10y,plot type=line, color=red, pen=-, linewidth=2");

if (L==10) {
    matlabload('apd_10um_IV_data_publication.mat');
    plot(Vpub,abs(Ipub_dark),"Bias [V]","Current [A]",num2str(L)+"um device","log10y,plot type=line, color=blue, pen=--, linewidth=2");
    plot(Vpub,abs(Ipub_photo),"Bias [V]","Current [A]",num2str(L)+"um device","log10y,plot type=line, color=red, pen=--, linewidth=2");
    legend("DEVICE dark current", "DEVICE photo current","Publication dark current","Publication photo current");
    holdoff;
}
else{
    legend("DEVICE dark current", "DEVICE photo current");
    holdoff;
}

############################ 3. Calculate responsivity. #####################

###### 3.1 Interpolate dark and photo currents on the same voltage grid. ###########
Vrwcl = 0;
if(abs(Vdark_wcl(length(Vdark_wcl))) < abs(Vphoto_wcl(length(Vphoto_wcl)))){
    IphotoInterp_wcl = interp(Iphoto_wcl,Vphoto_wcl,Vdark_wcl);
    IdarkInterp_wcl = Idark_wcl;
    Vrwcl = Vdark_wcl;
}else{
    IphotoInterp_wcl = Iphoto_wcl;
    IdarkInterp_wcl = interp(Idark_wcl,Vdark_wcl,Vphoto_wcl);
    Vrwcl = Vphoto_wcl;
}

####### 3.2 Calculate responsivity R ########
Pin = getnamed('CHARGE::generation Ge','scale factor');
Rwcl = abs(IphotoInterp_wcl - IdarkInterp_wcl)/Pin;

# Plot
plot(Vrwcl,Rwcl,"Voltage (V)","Responsivity (A/W)",num2str(L)+"um device","log10y,plot type=line, color=blue, pen=-, linewidth=2");

############################  4. Calculate the unit gain voltage as the root of second derivative ###############

##### 4.1 Calculate first derivative dIphoto_wcl/dV #####
Vd = zeros(length(Vphoto_wcl)-1);
Id = zeros(length(Vphoto_wcl)-1);
for(i=1:length(Vphoto_wcl)-1){
    Vd(i) = (Vphoto_wcl(i)+Vphoto_wcl(i+1))/2.0;
    dV = (Vphoto_wcl(i+1)-Vphoto_wcl(i));
    Id(i) = (Iphoto_wcl(i+1)-Iphoto_wcl(i))/dV;
}

##### 4.2 Calculate second derivative d^2Iphoto_wcl/dV^2 #####
Vd2 = zeros(length(Vd)-1);
Id2 = zeros(length(Vd)-1);
for(i=1:length(Vd)-1){
    Vd2(i) = (Vd(i)+Vd(i+1))/2.0;
    dV = (Vd(i+1)-Vd(i));
    Id2(i) = (Id(i+1)-Id(i))/dV;
}

##### 4.3 Find zero in second derivative #####
zeroInd = find(Id2,0);
Vmg1_eth = Vd2(zeroInd);
Img1_eth = interp(Iphoto_wcl,Vphoto_wcl,Vmg1_eth);

############################## 5. Find gain vs bias ####################################
MG_eth = Iphoto_wcl/Img1_eth;
Vmg_eth = Vphoto_wcl;

#Plot
plot(Vmg_eth,abs(MG_eth),"Bias [V]","Multiplication gain [unitless]",num2str(L)+"um device","log10y,plot type=line, color=blue, pen=-, linewidth=2");