% Fig_2_5.m
% Essential Semiconductor Laser Device Physics
% Plot spontaneous emission and gain with Lorentzian and Sech
% broadening as function of photon energy.
% Calculate broadened gain directly and via spontaneous emission with 
% Lorentzian and Sech and note the differences.
% Uses function mu.m and fermi.m
% Carrier density n(m-3), temperature kelvin(K)
% 
function [main]=main
clear all
clf;   

%plotting parameters + fontsizes
FS      = 16;	%label fontsize
FSN     = 14;	%number fontsize
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);

hbartxt=['\fontname{MT Extra}h\fontname{Times New Roman}'];

%*******************************
% Declare Constants and Arrays *
%*******************************
echarge =1.60219e-19;				%electron charge (C)
hbar    =1.0545928e-34;				%Planck's constant (J s)
c       =2.997925e8;				%speed of light in vacuum (m s-1)
kB      =8.617e-5;					%Boltzmann constant (eV K-1)
epsilon0=8.85419e-12;			    %permittivity of free space (F m-1)
 
   m0=9.10956e-31;					%Bare electron mass (kg)
   me=0.065*m0;						%Effective electron mass (kg) 0.07
   mhh=0.65*m0;						%Effective hhole mass (kg) 0.50
   mr=1/(1/me+1/mhh);				%Reduced electron mass (kg)
   rerr=1e-3;						%Relative error
   
   nr=3.3;							%Refractive index
   Eg=0.968;						%Band gap energy (eV)

   n=4.0e18;						%Carrier density (cm-3)
   ncarrier=n*1e6;					%Convert carrier density to (m-3)
   
   kelvin=328.0;					%Temperature (K)
   kBT=kB*kelvin;					%Thermal energy (eV)
   beta=1/kBT;						%Inverse thermal energy (eV-1)
   gamma=0.025						%Broadening (eV) 0.002 0.01 0.02 0.05
   
   %const=-echarge^2*hbar/(2*nr*c*m0^2*epsilon0)*((2*mr/hbar^2)^1.5);
   %const=2.64e4;%constant that gives ~ 330 cm-1 gain at 2e18 cm-3 in GaAs
   const=3.00e4;%constant that gives ~ 1000 cm-1 gain at 4e18 cm-3 in InGaAsP
   

   muhh=mu(mhh,ncarrier,kelvin,rerr);%Call mu chemical potential function for holes
   mue =mu(me, ncarrier,kelvin,rerr);%Call mu chemical potential function for electrons
   deltamu=mue+muhh					%Delta mu
   
   deltae=0.001;		%1 meV energy step, deltae < = gamma
   npoints=500;         %number of points in plot
   %Emin=0.1			%start plot from this energy below bandgap Eg (eV)
   %Emax=0.3			%end plot at this energy above bandgap Eg (eV)
   
   %pre-allocate 
   Energy       =zeros(1,npoints);
   rspon        =zeros(1,npoints);
   gain1        =zeros(1,npoints);
   Ephoton      =zeros(1,npoints);
   rsponbrd     =zeros(1,npoints);
   rsponbrd2    =zeros(1,npoints);
gain_rsp_lore   =zeros(1,npoints);
gain_rsp_sech2  =zeros(1,npoints);
gain_dir_sech2  =zeros(1,npoints);
gain_dir_lore2  =zeros(1,npoints);
   
   	for m=1:npoints
         Energy(m)=0;
         rspon(m)=0;
         gain1(m)=0;
         
    end
     	%Calculate spontaneous emission with no Lorentzian broadening
    for j=150:1:npoints
      	Energy(j)=(j*deltae)-0.15;		%Photon energy - Eg
      	Ehh=(Energy(j))/(1+mhh/me);		%Energy in hole band
      	Ee=(Energy(j))/(1+me/mhh);		%Energy in conduction band
      	fhh=fermi(beta,Ehh,muhh);		%Call Fermi function for holes
      	fe=fermi(beta,Ee,mue);			%Call Fermi function for electrons
         rspon(j)=(const)*(Energy(j)^0.5)*(fe*fhh);
         gain1(j)=rspon(j)*(1-exp((Energy(j)-deltamu)*beta));
   	end
      %Calculate spontaneous emission with Lorentzian broadening
      %Integrate 
      for j=1:npoints
         Ephoton(j)=(Eg-0.15)+(j*deltae);
%
%rsponbrd(j)=QUAD('rspint',0,Eg,.001,0,Eg,gamma,Ephoton,mhh,me,muhh,mue,beta);
%
         rspint=0;
         rspint2=0;
         gain_dir_sech=0;
         gain_dir_lore=0;
         E=0;
         for k=1:1400
E=E+deltae;
Ehh=E/(1+mhh/me);
Ee=E/(1+me/mhh);
fhh=fermi(beta,Ehh,muhh);
fe=fermi(beta,Ee,mue);
rspint=rspint+(sqrt(E))*fhh*fe*deltae*(gamma/2/pi)/((Eg+E-Ephoton(j))^2+(gamma/2)^2);
rspint2=rspint2+(sqrt(E))*fhh*fe*deltae*sech(abs(Eg+E-Ephoton(j))/gamma)/(gamma*pi);
gain_dir_sech=gain_dir_sech+(sqrt(E))*(fhh+fe-1)*deltae*sech((Eg+E-Ephoton(j))/gamma)/(gamma*pi);
gain_dir_lore=gain_dir_lore+(sqrt(E))*(fhh+fe-1)*deltae*(gamma/2/pi)/((Eg+E-Ephoton(j))^2+(gamma/2)^2);
      end
rsponbrd(j)=(const)*rspint;%spontaneous emission with Lorentzian broadening
rsponbrd2(j)=(const)*rspint2;%spontaneous emission with Sech broadening
gain_rsp_lore(j)=rsponbrd(j)*(1-exp(((Ephoton(j)-Eg)-deltamu)*beta));%gain from rspn with Lorentzian broadening
gain_rsp_sech2(j)=rsponbrd2(j)*(1-exp(((Ephoton(j)-Eg)-deltamu)*beta));%gain from rspn with sech broadening
gain_dir_sech2(j)=(const)*gain_dir_sech;%gain with sech broadening
gain_dir_lore2(j)=(const)*gain_dir_lore;%gain with Lorentzian broadening

      end   
%vector AA(1,nTemp) contains legend strings
AA(1,1)={('rsp no broadening')};
AA(1,2)={('rsp Lorentzian')};
AA(1,3)={('rsp Sech')};
%********************  plot figures ************************
   figure(1)
   plot(Energy, rspon,'r','LineWidth',LW);%rspon no broadening
   hold on;
   plot(Ephoton-Eg, rsponbrd,'g','LineWidth',LW);%rspon Lorenzian broadening
   plot(Ephoton-Eg, rsponbrd2,'b','LineWidth',LW);%rspon Sech broadening
   hold off
   axis([-.1,.3,0,1500]);
   xlabel(['Photon energy, ',hbartxt,'\omega (eV)']);
   ylabel('Spontaneous emission, \itrspon\rm (arb.)');
   title(['n=',num2str(n/1e18),'x10^{18} cm^{-3},  T=',num2str(kelvin),...
       ' K,  m_e=',num2str(me/m0),', m_{hh}=',num2str(mhh/m0),...
       ', \gamma=',num2str(gamma),' eV']);
   grid on;
   legend(AA,'Location','NorthEast');   %print legend from string vector AA 
   
%s=char('y','b','r','g','c','m','k');	%plot curves in different colors

%vector A(1,nTemp) contains legend strings
A(1,1)={('g no broadening')};
A(1,2)={('g rsp Lorentzian')};
A(1,3)={('g rsp Sech')};
A(1,4)={('g direct Lorentzian')};
A(1,5)={('g direct Sech')};
   
   figure(2)
   plot(Energy, gain1,'r','LineWidth',2);
   hold on;
   plot(Ephoton-Eg, gain_rsp_lore,'g','LineWidth',LW);
   plot(Ephoton-Eg, gain_rsp_sech2,'b','LineWidth',LW);
   plot(Ephoton-Eg, gain_dir_lore2,'m','LineWidth',LW);
   plot(Ephoton-Eg, gain_dir_sech2,'k','LineWidth',LW);
   hold off;
   axis([-.15,.15,-1200,1200]);
   xlabel(['Photon energy, ',hbartxt,'\omega (eV)']);
   ylabel('gain,  \itg\rm_{opt} (cm^-^1)');
   title(['n=',num2str(n/1e18),'x10^{18} cm^{-3},  T=',num2str(kelvin),...
       ' K,  m_e=',num2str(me/m0),', m_{hh}=',num2str(mhh/m0),...
       ', \Delta\mu=',num2str(deltamu),' eV, \gamma=',num2str(gamma),' eV']);
   grid on;
   legend(A,'Location','NorthWest');   %print legend from string vector A