% Fig_4_8.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=2000;				%number of steps in time integration
curentI =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) 2.5e-16
epsi       = 5.00E-18;	%gain compression (cm3) 5e-18
beta       = 1.00E-04;	%spontaneous emission coefficient 10^-4
gamma_cons = 0.25;		%optical confinement factor 0.25
wavelength = 1.31;		%optical emission wavelength (um) 1.31

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


Imin	   = 0;         %start value of current (mA)
Imax   	   = 20;		%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

temptxt=['\itn\rm_0=',num2str(n0density*1000),...
'\times10^{18} cm^{-3}, \beta=',...
num2str(beta),', \itr\rm_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;
current=0;
dcurrent=(Imax-Imin)/ngplot;

for j=1:ngplot
   current=current+dcurrent;  
   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);
plot(currentI,photonL,'r-','LineWidth',LW);
xlabel(['Current, \itI\rm_{inj} (mA)']);
ylabel(['Light, \itL\rm_{out} (mW/facet)']);
title(temptxt);
grid on;

figure(2);
loglog(currentI,photonL,'r-','LineWidth',LW);
%plot(log10(currentI),log10(photonL),'r-','LineWidth',LW);
title(temptxt)
xlabel(['log_{10}(\itI\rm_{inj}) (mA)']);
ylabel(['log_{10}(\itL\rm_{out}) (mW/facet)']);
grid on;

figure (3);
plot(currentI,carrierN,'b','LineWidth',LW);
xlabel(['Current, \itI\rm_{inj} (mA)']);
ylabel(['Carrier density, \itn\rm (\times 10^{18} cm^{-3})']);
title(temptxt);
grid on;

figure (4);
set(gca,'Fontsize',FSplot);                     %set current gca font size
hl1 = line(currentI,carrierN,'Linewidth',LWplot,'Color','k');
ax1 = gca;                                      %ax1 is first axis
set(ax1,'XColor','k','YColor','k','Linewidth',LWplot,...
    'TickLength',[TLplot TLplot]);
ylabel(['Carrier density, \it n \rm (\times 10^{18} cm^{-3})']);

ax2 = axes('Position',get(ax1,'Position'),...   %ax2 is second axis
'YAxisLocation','right','XAxisLocation','bottom',...
'Color','none',...
'XColor','k','YColor','r','Linewidth',LWplot,'TickLength',[TLplot TLplot]);
set(gca,'Fontsize',FSplot);                     %set current gca font size
hl2 = line(currentI,photonL,'Linewidth',LWplot,'Color','r','Parent',ax2);

xlabel(['Current, \it I \rm_{inj} (mA)']);
ylabel(['Light, \it L \rm_{out} (mW/facet)']);
title(temptxt);
grid on