% Fig_4_7.m
% Essential Semiconductor 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);

%declare 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=10000;				%number of time-steps

%********************  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)
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)
beta       = 1.00E-04;		%spontaneous emission coefficient 10^-4
gamma_cons = 0.25;			%optical confinement factor
wavelength = 1.31;			%optical emission wavelength (um)

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

pstart1    = 2.5;			%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;

%********************  the main integration loop ************************

for i = 1 : 1 : nsteps
   time(i) = ( double(i-1) +1.0 ) * tincrement ;
   current(i) = rioff1 ;
   
   if time(i) > pstart1
      current(i) = rion1 ;
   end
   
   if time(i) > pstart2
      current(i) = rion2 ;
   end
   
[sdensity, ndensity]=runge4(current(i),sdensity,ndensity,tincrement,data);
         carrier(i) = ndensity * 1.0e21/1.0e18	;
          photon(i) = (lightout*sdensity)	;
end;%main ends here.

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 ************************
%time-evolution of current, light output and carrier density.
%
%Injected current change in time.
figure(1);
hold on;
axis([0 time(nsteps) 0 max(current)*1.05]);

temp=['Injected current as function of time'];
title(temp, 'fontsize', FS)

xlabel(['Time, \itt\rm (ns)'],'fontsize', FS);
ylabel(['Current,  \itI\rm (mA)'],'fontsize', FS);
grid on;
plot(time,current,'r-','LineWidth',LW);

hold off;

%Evolution of light output from mirror one.
figure(2);
hold on;
axis([0 time(nsteps) 0 max(photon)*1.05]);

temp=['Optical power from mirror-one as function of time'];
title(temp, 'fontsize', FS)

xlabel(['Time, \itt\rm (ns)'],'fontsize', FS);
ylabel(['Power,  \itL\rm (mW/facet)'],'fontsize', FS);
grid on;

plot(time,photon,'r-','LineWidth',LW);
hold off;

%Evolution of carrier density.
figure(3);
hold on;
axis([0 time(nsteps) min(carrier) max(carrier)*1.05]);

temp=['Carrier density as function of time'];
title(temp, 'fontsize', FS);

grid on;
xlabel(['Time, \itt\rm (ns)'],'fontsize', FS);
ylabel(['Carrier density,  \itn\rm (10^1^8cm^-^3)'],'fontsize', FS);

plot(time,carrier,'r-','LineWidth',LW);
hold off;   

%Light output as function of carrier density.
figure(4);
hold on;
axis([min(carrier) max(carrier)*1.05 0 max(photon)*1.05]);

temp=['Light output as function of carrier density'];
title(temp, 'fontsize', FS);

grid on;
xlabel(['Carrier density, \itn\rm (10^1^8cm^-^3)'],'fontsize', FS);
ylabel(['Power, \itL\rm_{out} (mW/facet)'],'fontsize', FS);

plot(carrier,photon,'r-','LineWidth',LW);
hold off;   






