% Chapt9Exercise1.m
% Plot gain and spontaneous emission as function of photon energy
% for carrier 10 different carrier densities
% uses function mu.m and fermi.m
% carrier density n(m-3), temperature kelvin(K)
% 
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
   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
   
   for k=1:1:10
      
   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 function for holes
   mue=mu(me,ncarrier,kelvin,rerr)   %Call mu chemical potential function for electrons
   deltamu=mue+muhh
      
   const=2.64e4;	%GaAs constant gives gain 330 cm-1 at n = 2E18 cm-3
   
   deltae=0.001;
   	for j=1:300
      	Energy(j)=j*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

      
   figure(1)
   hold on;
   plot(Energy+Eg, gain);
   %axis([1.4 1.7 -7000 7000]);
   xlabel('Photon energy, hw (eV)');
   ylabel('Optical gain, g (cm-1)');
	title(['n(min)=',num2str(1),'x10^1^8 cm^-^3, n(max)=',...
        num2str(k),'x10^1^8 cm^-^3, me=',num2str(me/m0),...
        ', mhh=',num2str(mhh/m0),', T=',num2str(kelvin),...
        ' K, Eg=',num2str(Eg),' eV']);
   grid on;
   hold off;

   figure(2)
   hold on;
   plot(Energy+Eg,rspon,'r');
   xlabel('Photon energy, hw (eV)');
   ylabel('Spontaneous emission, rsp (arb.)');
	title(['n(min)=',num2str(1),'x10^1^8 cm^-^3, n(max)=',...
        num2str(k),'x10^1^8 cm^-^3, me=',num2str(me/m0),...
        ', mhh=',num2str(mhh/m0),', T=',num2str(kelvin),...
        ' K, Eg=',num2str(Eg),' eV']);
   grid on;
   end
   hold off;

   
   