#############################################################################
# Calculate IV curves and responsivity for the application example for 
# avalanche photodiode (APD) from H. T. Chen, et. al., IEEE Photonics Journal, 
# vol. 7, no. 4, p. 7902909, Aug. 2015.
#
# 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  (norm length)
biased_contact = 'anode';
plot_contact = 'anode';
##### Load results from the publication. ######
basefilename = 'apd5V_10um';
matlabload(basefilename+'_IV_publication.mat');
matlabload(basefilename+'_gain_publication.mat');

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

load('apd5V_electrical_'+num2str(L)+'um_dark.ldev');

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

Idark = -2*Itotal;
Vdark = Vtotal;

N = length(Vdark);
plot(Vdark(2:N),abs(Idark(2:N)),"Bias [V]","Current [A]","","log10y,plot type=line, color=blue, pen=-, linewidth=2");
holdon;

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

load('apd5V_electrical_'+num2str(L)+'um_photo.ldev');

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

Iphoto = -2*Itotal;
Vphoto = Vtotal;

plot(Vphoto,abs(Iphoto),"Bias [V]","Current [A]","","log10y,plot type=line, color=red, pen=-, linewidth=2");

plot(Vpub,abs(Ipub_dark),"Bias [V]","Current [A]","","log10y,plot type=line, color=blue, pen=--, linewidth=2");
plot(Vpub,abs(Ipub_198dbm),"Bias [V]","Current [A]","","log10y,plot type=line, color=red, pen=--, linewidth=2");
legend("DEVICE dark current", "DEVICE photo current","Publication dark current","Publication photo current");
holdoff;

############################ 3. Calculate and plot responsivity. #####################

###### 3.1 Interpolate dark and photo currents on the same voltage grid. #####
V = 0;
if(abs(Vdark(length(Vdark))) < abs(Vphoto(length(Vphoto)))){
    IphotoInterp = interp(Iphoto,Vphoto,Vdark);
    IdarkInterp = Idark;
    V = Vdark;
}else{
    IphotoInterp = Iphoto;
    IdarkInterp = interp(Idark,Vdark,Vphoto);
    V = Vphoto;
}

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


##### 3.3 Plot responsivity #####
plot(V,abs(R),"Bias [V]","Responsivity [A/W]","","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/dV #####
Vd = zeros(length(Vphoto)-1);
Id = zeros(length(Vphoto)-1);
for(i=1:length(Vphoto)-1){
    Vd(i) = (Vphoto(i)+Vphoto(i+1))/2.0;
    dV = (Vphoto(i+1)-Vphoto(i));
    Id(i) = (Iphoto(i+1)-Iphoto(i))/dV;
}

##### 4.2 Calculate second derivative d^2Iphoto/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,Vphoto,Vmg1_eth);

############################# 5. Calculate gain vs. bias #######################################
 
##### 5.1 Find gain vs bias ########
MG_eth = Iphoto/Img1_eth;
Vmg_eth = Vphoto;

##### 5.2 Plot gain vs bias  #######
plot(Vmg_eth,abs(MG_eth),"Bias [V]","Multiplication gain [unitless]","","log10y,plot type=line, color=blue, pen=-, linewidth=2");
holdon;
plot(V_pub,abs(Gain_pub),"Bias [V]","Multiplication gain [unitless]","","log10y,plot type=line, color=blue, pen=--, linewidth=2");
legend("DEVICE gain", "Publication gain");