% Chapt7Exercise3.m
% uses function fermi.m
% carrier density n(cm-3), temperature kelvin(K), relative error rerr
% returns chemical potential mu1 measured in units of meV

n=1.e18;					%carrier density (cm-3)
kelvin=300.0;				%absolute temperature (K)
rerr=1.e-3;					%relative error

m0=9.10956;					%bare electron mass (kg x 10^31)
m1=0.07;					%effective electron mass
echarge=1.60219;			%electron charge (C x 10^19)
hbar=1.0545928;				%Planck's constant (J s x 10^34)
kB=8.617e-5;				%Boltzmann's constant (eV K-1)
   
kF1=(3.0*(pi^2.0)*n)^(1/3);	%Fermi wave vector (cm-1)
lambdaF=2.0*pi/kF1;
eF=1.e-11*((hbar*kF1)^2)/(2.0*m0*m1*echarge);%Fermi energy (meV)
kBT=1000.*kelvin*kB;		%thermal energy (meV)
beta=1./kBT;				%inverse thermal energy (meV-1)
   
mumax=eF;					%maximum possible value of chemical potential
x=((n*1.e6)/2.)*(((2.e-15*pi*beta*(hbar^2))/(echarge*m0*m1))^1.5);
mumin=(+1./beta)*log(x);
   
emax=eF+(15./beta);			%maximum limit of integration
de=emax/1000.;				%energy step
const=1.e16*((20.*echarge*m1*m0)^.5)*echarge*m1*m0/((pi^2)*(hbar^3));
   
for j=1:25
   mu1=mumin+((mumax-mumin)/2.);
	energy=0.0;
	ainter=0.0;
   
%	calculate carrier density n'
	for i=1:1000;
		energy=energy+de;
		ainter=ainter+(((sqrt(energy))*de)*fermi(beta,energy,mu1));
   end;
   nprime=const*ainter;
%  delta is relative error in carrier density
   delta=(n-nprime)/n;
   	if((abs(delta)) < rerr)
      	break;
     			elseif(delta < 0.) 
         	mumax=mu1;
			else 
			mumin=mu1;
   	end;
end;
%	print output: chemical potential mu1(meV), Fermi energy eF(meV), temperature kelvin(K) 
%	kBT (meV), carrier density nprime(cm-3), number of iterations j
ttl1=['chemical potential 			= ',num2str(mu1),' meV'];
ttl2=['Fermi energy 				= ',num2str(eF),' meV'];
ttl3=['Fermi wave vector 			= ',num2str(kF1),' cm-1'];
ttl4=['Fermi wavelength 			= ',num2str(lambdaF),' cm'];
ttl5=['Temperature                  = ',num2str(kelvin),' K'];
ttl6=['kBT                          = ',num2str(kBT),' meV'];
ttl7=['carrier density              = ',num2str(nprime),' cm-3'];
ttl8=['number of iterations 		= ',num2str(j)'];
%	compare with second-order mu2(meV) and fourth-order mu4(meV) expansion
ttl9=['chemical potential (second-order) 	= ',num2str(eF-(((pi*kBT)^2)/(12.*eF))),' meV'];
ttl10=['chemical potential (fourth-order) 	= ',num2str(eF-(((pi*kBT)^2)/(12.*eF))-(7*((pi*kBT)^4)/(960.*(eF^3)))),' meV'];

Solution =strvcat(ttl1,ttl2,ttl3,ttl4,ttl5,ttl6,ttl7,ttl8,ttl9,ttl10)
	if j >= 25
   	'check convergence!'
	end;


