% Chapter9Exercise3.m
% Plot spontaneous emission and gain with Lorentzian 
% broadening as function of photon energy.
% uses function mu.m and fermi.m
% carrier density n(m-3), temperature kelvin(K)
% 
function [main]=main
clear
clf   

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 (kg)
   rerr=1e-3;						%Relative error
   
   nr=3.3;							%Refractive index
   Eg=1.4;							%Band gap energy (eV)

	n=2.0e18;						%Carrier density (cm-3)
   ncarrier=n*1e6;					%Convert carrier density to (m-3)
   
   kelvin=300.0;					%Temperature (K)
   kBT=kB*kelvin;					%Thermal energy (eV)
   beta=1/kBT;						%Inverse thermal energy (eV-1)
   gamma=0.015						%Broadening (eV)
   
   				%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
   

   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
   %Emin=0.1			%start plot from this energy below bandgap Eg (eV)
   %Emax=0.3			%end plot at this energy above bandgap Eg (eV)
   
   	for m=1:400
         Energy(m)=0;
         rspon(m)=0;
        % gain(m)=0;
         gain1(m)=0;
         
      end
     	%Calculate spontaneous emission with no Lorentzian broadening
      for j=100:1:400
      	Energy(j)=(j*deltae)-0.1;		%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:400
         Ephoton(j)=(Eg-0.1)+(j*deltae);
         %
         %rsponbrd(j)=QUAD('rspint',0,Eg,.001,0,Eg,gamma,Ephoton,mhh,me,muhh,mue,beta);
         %
         rspint=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);
         end
rsponbrd(j)=(const)*rspint;%spontaneous emission with broadening
gain(j)=rsponbrd(j)*(1-exp(((Ephoton(j)-Eg)-deltamu)*beta));%gain with broadening

      end   
        
   figure(1)
   plot(Ephoton-Eg, rsponbrd,'b');
   hold on;
   plot(Energy, rspon,'r');
   hold off
   axis([-.1,.3,0,1000]);
   xlabel('Photon energy, (eV)');
   ylabel('Spontaneous emission, rspon (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;
      
   figure(2)
   plot(Ephoton-Eg, gain,'b');
   hold on;
   plot(Energy, gain1,'r');
   hold off;
   axis([-.1,.1,-500,500]);
	xlabel('Photon energy, (eV)');
   ylabel('gain, (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;
   
   