% Fig_4_9.m
% Essential Semiconductor Laser Device Physics
% solves single mode laser rate equations
% using 4th order Runge-Kutta method for fixed time increment and plots
% light output, L, and carrier density, n, 
% as a function of injection current, I.

% 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

%  start the program and clear all the previous definitions.
clear all
clf;

FS      = 18;	%label FontSize 18
FSN     = 16;	%number FontSize 16
LW      = 2;	%LineWidth
FSplot  = 14;   %font size for plots 
LWplot  = 2;    %line width for plots
TLplot  = 0.03; %tick length for plots

% 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 constant [J s]
echargec=1.6021e-19;	  	%electron charge [Coulomb]
vlightcm=2.997925e10;  		%velocity of light [cm s-1]
pi2=2*pi;			  	  	%value of 2pi

ngplot=100;					%number of points plotted
nsteps=40000;				%number of steps in time integration
currentI=zeros(1,ngplot);
carrierN=zeros(1,ngplot);
photonL =zeros(1,ngplot);
%********************  start the main program ***************************
% input 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-05;	%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)

Imin	  = 1;          %start value of current (mA) 1
Imax   	  = 40;			%end value of current (mA) 20

% 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

s=char('k','b','g','r','m','c','y');	%plot curves in different colors
temptxt=['n_0=',num2str(n0density*1000),...
    '\times10^{18} cm^{-3}, \beta_{min}=',...
num2str(beta*10),', r_1=',num2str(mirrone),', \alpha_i=',...
num2str(alfa_i),' cm^{-1}'];
%********************  the main integration loop ************************

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;

photonold=0;
dcurrent=(Imax-Imin)/ngplot;%increment in current (mA)

for ii=1:1:5
    data(10)=data(10).*10;   %increment value of beta
    current=Imin;            %reset current
for j=1:ngplot
   current=current+dcurrent; %increment current (mA)
   
   for i = 1:nsteps;   
[sdensity, ndensity]=runge4( current, sdensity, ndensity, tincrement, data);
   	carrier = ndensity * 1.0e21/1.0e18	;
      photon = (lightout*sdensity);
      photonold=photon;
   end;
currentI(j)=current;
photonL(j)=photon;
carrierN(j)=carrier;

end
%********************  plot figures ************************

figure(1);
loglog(currentI,photonL,s(ii),'LineWidth',LW);
title(temptxt)
xlabel(['log_{10}(\itI\rm_{inj}) (mA)']);
ylabel(['log_{10}(\itL\rm_{out}) (mW/facet)']);
hold on;
grid on;

figure(2);
loglog(currentI,carrierN,s(ii),'LineWidth',LW);
title(temptxt)
xlabel(['log_{10}(\itI\rm_{inj}) (mA)']);
ylabel(['log_{10}(\itn\rm) (\times 10^{18} cm^{-3})']);
hold on;
grid on;

end
hold off


