% Fig_2_4.m
% Essential Semiconductor Laser Device Physics
% Plot gain and spontaneous emission as function of photon energy
% for carrier 3 different carrier densities
% uses function mu.m and fermi.m
% carrier density n(m-3), temperature kelvin(K)
% 
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);

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

%********************       constants        ****************************
   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.07*m0;						%Effective electron mass kg
   mhh=0.5*m0;						%Effective hhole mass kg
   mr=1/(1/me+1/mhh);				%Reduced electron mass
   rerr=1e-3;						%Relative error
   
   nr=3.3;							%Refractive index
   Eg=1.4							%Band gap energy eV

   kelvin=300.0;					%Temperature K
   kBT=kB*kelvin;					%Thermal energy eV
   beta=1/kBT;						%Inverse thermal energy eV-1
   
   npoints=300;                     %number of points in plot  
   const=2.64e4;	%GaAs constant gives gain 330 cm-1 at n = 2E18 cm-3
   deltae=0.0004;   %energy increment
   Emax=Eg+(deltae*npoints);        %maximum energy
   
   for k=1:1:3
      
n=k*1.e18;                          %Carrier density cm-3
ncarrier=n*1e6;                     %Convert carrier density to m-3

muhh=mu(mhh,ncarrier,kelvin,rerr)   %Call mu chemical potential for holes
mue=mu(me,ncarrier,kelvin,rerr)     %Call mu chemical potential for electrons
deltamu=mue+muhh
   
   	for j=1:npoints
      	Energy(j)=(j-1)*deltae;			%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
      	gain(j)=const*(Energy(j)^0.5)*(fe+fhh-1);
      	rspon(j)=(const)*(Energy(j)^0.5)*(fe*fhh);
   	end
ttl=(['\itn\rm_{min}=',num2str(1),'\times10^{18}cm^-^3, \itn\rm_{max}=',...
        num2str(k),'\times10^{18}cm^-^3, \itm\rm_e=',num2str(me/m0),...
        '\itm\rm_0, \itm\rm_{hh}=',num2str(mhh/m0),...
        'm_0, \itT\rm=',num2str(kelvin),...
        'K, \itE\rm_g=',num2str(Eg),'eV']);
%********************  plot figures ************************      
figure(1)
hold on;
plot(Energy+Eg, gain, 'LineWidth',LW);
axis([Eg Emax -1500 1000]);
xlabel(['Photon energy, ',hbartxt,'\omega (eV)']);
ylabel('Optical gain,  \itg\rm_{opt} (cm^{-1})');
title(ttl);
grid on;
hold off;

figure(2)
hold on;
plot(Energy+Eg,rspon,'r','LineWidth',LW);
axis([Eg Emax 0 1500]);
xlabel(['Photon energy, ',hbartxt,'\omega (eV)']); 
ylabel('Spontaneous emission, \itrsp\rm (arb.)');
title(ttl);
grid on;
end
hold off;