% Fig_4_11.m
% Essential Semionductor Laser Device Physics
% solves single-mode laser diode rate equations
% calls runge4.m which is a 4th order Runge-Kutta integrator
% for fixed time increment
% plots carrier density and photon density as function of time.

% gain=g=gamma*gslope*(n-n0)*(1-epsi*s)
% carrier density, n, current, I, and photon density, s, are related via 
% fcn=n'=(I/ev)-(n/tau_n)-(g*s)
% fcs=s'=(g-K)*s+beta*Rsp

clear all
clf

FS      = 18;	%label FontSize 18
FSN     = 16;	%number FontSize 16
LW      = 2;	%LineWidth

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times New Roman');
set(0,'DefaultAxesFontSize', FSN);

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times New Roman');
set(0,'DefaultTextFontSize', FSN);

%********************       constants        ****************************

hconstjs=6.6262e-34;   		%Planck's constant [J s]
echargec=1.6021e-19;	  	%electron charge [Coulomb]
vlightcm=2.997925e10;  		%velocity of light [cm s-1]

nsteps=20000;				%number of time-steps 10000

%********************  start the main program ***************************
%assign values to parameters

ngroup     = 4;				%refractive index

clength    = 3.00E-02;		%cavity length (cm)
thick      = 1.40E-05;		%thickness of active region (cm)
width      = 8.00E-05;		%width of active region (cm)

tincrement = 1.00E-12;		%time increment (s)
initialn   = 0.00E+18;		%initial value of carrier density (cm-3)

Anr        = 2.00E+08;		%non-radiative recombination rate (s-1)2e8
Bcons      = 1.00E-10;		%radiative recombination coefficent (cm3 s-1) 
Ccons      = 1.00E-29;		%non-linear recombination coefficient (cm6 s-1)

n0density  = 1.00E+18;		%transparency carrier density (cm-3)
gslope     = 2.50E-16;		%optical gain-slope coefficient (cm2 s-1)
epsi       = 5.00E-18;		%gain compression (cm3) 5e-18
beta       = 1.00E-05;		%spontaneous emission coefficient 10^-4
gamma_cons = 0.25;			%optical confinement factor 0.25
wavelength = 1.31;			%optical emission wavelength (um)

mirrone    = 0.32;			%optical reflectivity from first mirror  0.32
mirrtwo    = 0.32;			%optical reflectivity from second mirror 0.32
alfa_i     = 40;			%internal optical loss (cm-1) 40

pstart1    = 1.0;			%start of first current step (ns)
pstart2    = 3.5;			%start of second current step (ns)

rioff1     = 0;				%off value of current (mA)
% rion1      = 20;			%first step in current (mA)
% rion2      = 20;			%second step in current (mA)

%initialize sdensity, ndensity and calculate total optical loss [cm-1]
	sdensity = 0.0;
	ndensity = initialn * 1.0e-21;

	alfa_m   = log ( 1.0/(mirrone * mirrtwo) ) / ( 2.0 * clength ) ;
	alfasum  = alfa_m + alfa_i ;

%rescale integrator: time[ns], current[mA], length[nm] and 1nm3=1e-21cm3
echarge   	= echargec*1.0e9*1.0e3  ;% electron charge		
vlight    	= vlightcm/1.0e9/1.0e-7 ;% velocity of light [m/s] to [nm s-1]	

tincrement 	= tincrement * 1.0e9 ;
gslope    	= gslope * 1.0e14 ;
n0density 	= n0density * 1.0e-21 ;
   
%length scale in nm.
   clength  = clength / 1.0e-7 ;
	width   = width / 1.0e-7 ;
   thick	= thick / 1.0e-7 ;
   
%recombination constants   
	Anr     = Anr / 1.0e9 ;
	Bcons   = Bcons / (1.0e9*1.0e-21) ;
   Ccons   	= Ccons / (1.0e9*1.0e-42) ;
   
	epsi    = epsi / 1.0e-21 ;
	gslope  = gslope * vlight / ngroup ;
	alfasum = alfasum * 1.0e-7 ;
	kappa   = alfasum * vlight / ngroup ;
	volume  = clength * width * thick; 

%scale light output to mW from mirror one
tmp = hconstjs * (vlightcm * 1.0e-2)/(wavelength * 1.0e-6);% use m [J]	
tmp	= tmp * alfa_m * vlightcm;						% use cm [s-1]
tmp	= tmp * 1000.0*volume/ngroup;					% use mW and nm [nm3]
lightout=tmp*(1.0-mirrone)/(2.0-mirrone-mirrtwo);	% output mirror one

data(1)=gamma_cons;
data(2)=gslope;
data(3)=n0density;
data(4)=epsi;
data(5)=kappa;
data(6)=Bcons;
data(7)=Anr;
data(8)=Ccons;
data(9)=volume;
data(10)=beta;

Imax=20.0;    %maximum current (mA) 20
nIstep=80.0;  %number of step current values 10
%********************  initialize arrays ************************
current =ones(1,nsteps)*rioff1;
rion1   =ones(1,nIstep);
carrier =zeros(nIstep,nsteps);
photon  =zeros(nIstep,nsteps);
tdelay  =zeros(1,nIstep);
timep=linspace(1,nsteps,nsteps)*tincrement;

temptxt=['\itn\rm_0=',num2str(n0density*1000),...
    '\times10^{18} cm^{-3}, \beta_{min}=',...
    num2str(beta*10),', \itr\rm_1=',num2str(mirrone),', \alpha_i=',...
    num2str(alfa_i),' cm^{-1}'];
%********************  increment maximum current  ************************
for kk=1:1:5
    data(10)=data(10).*10;
    
    for ii = 1:nIstep
sdensity=0;
ndensity=0;
   
rion1(ii)=(Imax*double(ii)/nIstep);
%************************************************************************ 
%********************  the main integration loop ************************
%************************************************************************ 
for i = 1 : 1 : nsteps
   if timep(i) > pstart1
      current(i) = rion1(ii) ;
   end
[sdensity,ndensity]=runge4( current(i),sdensity,ndensity,tincrement,data);
   carrier(ii,i) = ndensity * 1.0e21/1.0e18	;
   photon(ii,i) = (lightout*sdensity)	;
end
%main ends here.
%******************** search for delay time ************************
%threshold defined a half the steady-state value
%find first index gt threshold
k=find(photon(ii,1:nsteps-1)>(photon(ii,nsteps)/2) , 1, 'first');
 tdelay(kk,ii)=double(k)*tincrement - pstart1;     % delay in ns

end
end
% photon_number=photon( nsteps - 1 )*volume / lightout 
% number of photons at last data point
% carrier_number=carrier( nsteps - 1 )*volume/1.0e3	 
% number of carriers at last data point

%******************** plot figures *************************** 
figure(1)
plot(rion1, tdelay,'LineWidth',LW)
xlabel(['Current, \itI\rm_{step} (mA)']);
ylabel(['Time delay, \itt\rm_d (ns)']);
%axis([0 Imax 0 max(tdelay)*1.05]);
title(temptxt,'FontSize',FS);
grid on
%************************************************************************ 
figure(2)
loglog(rion1, tdelay,'LineWidth',LW);
xlabel(['Current, \itI\rm_{step} (mA)']);
ylabel(['Time delay, \itt\rm_d (ns)']);
title(temptxt,'FontSize',FS);
%axis([0 Imax 0 max(tdelay)*1.05]);
grid on
%************************************************************************
figure(3)
plot(rion1, 1./tdelay,'LineWidth',LW);
xlabel(['Current, \itI\rm_{step} (mA)']);
ylabel(['Inverse time delay, 1/ \itt\rm_d (ns^{-1})']);
title(temptxt,'FontSize',FS);
%axis([0 Imax 0 max(tdelay)*1.05]);
grid on
%************************************************************************ 
%check time evolution of light output from mirror one.
figure(4);
hold on;
plot(timep,photon,'b-');
%axis([0 max(timep) 0 max(photon)*1.05]);
title(temptxt,'FontSize',FS)
xlabel(['Time, \itt\rm (ns)']);ylabel(['Power, \itL\rm_{out} (mW/facet)']);
grid on;
%************************************************************************ 
%check time evolution of carrier density.
figure(5);
plot(timep,carrier,'b-');
hold on;
%axis([0 max(timep) min(carrier) max(carrier)*1.05]);
temp=['Carrier density as function of time'];
title(temptxt,'FontSize',FS);
grid on;
xlabel(['Time, \itt\rm (ns)']);
ylabel(['Carrier density, \itn\rm (10^1^8cm^-^3)']);
hold off;  
%************************************************************************