%Fig2_11.m
%Essential Electron Transport for Device Physics
%plots chemical potential as function of temperature for
%carrier density n(cm-3), temperature kelvin(K), effective mass m1, 
%relative error rerr, uses function fermi.m
clear;
clf;
    nmin=1.e16;				%min. carrier density (cm-3)
nstepmax=100;               %number of curves for different carrier density
    nmax=nmin*nstepmax;
    
      m0=9.109382;			%bare electron mass (x10^31 kg)
      m1=0.07;				%effective electron mass
 echarge=1.6021764;         %electron charge (x10^19 C)
	hbar=1.05457159;		%Planck's constant (x10^34 J s)
      kB=8.61734e-5;		%Boltzmann constant (eV K^-1)
   
    rerr=1.e-3;				%relative error

nplotpoints=100;            %number of data points in each plot
  minKelvin=1;				%minimum temperature (K)
  maxKelvin=400;            %maximum temperature (K)
kelvinscale=linspace(minKelvin,maxKelvin,nplotpoints+1);
   ndensity=linspace(nmin,nmax,nstepmax);

  muminscale=-150;          %set minimum on plot scale
   		if minKelvin < 1
            error('check minKelvin >=1 K')
         end;
         if minKelvin > maxKelvin
   			error('check minKelvin < maxKelvin')
			end;

stepKelvin=(maxKelvin-minKelvin)/nplotpoints;	%Step in temperature (K)
kF1=(3.0*(pi^2.0)*nmax)^(1/3);

Emax=1.1*(1.e-11*((hbar*kF1)^2)/(2.0*m0*m1*echarge));%Maximum value on z-axis   
   
   const=1.e16*((20.*echarge*m1*m0)^.5)*echarge*m1*m0/((pi^2)*(hbar^3));
   for nstep=1:1:nstepmax
       n=ndensity(nstep);           %carrier density in plot (cm-3)								
     kF1=(3.0*(pi^2.0)*n)^(1/3);    %Fermi wave vector
   	  eF=1.e-11*((hbar*kF1)^2)/(2.0*m0*m1*echarge);%Fermi energy (meV)
   mumax=eF;			%The maxiumum value of the chemical potential

  for k = 1:1:nplotpoints+1
      kelvin(k)=((k-1)*stepKelvin)+minKelvin;
		kBT=1000.*kelvin(k)*kB;				%Thermal energy (meV)
   	beta=1./kBT;
   
		x=((n*1.e6)/2.)*(((2.e-15*pi*beta*(hbar^2))/(echarge*m0*m1))^1.5);
	mumin=(+1./beta)*log(x);	%Minimum value of chemical potential
   
   	 emax=eF+(15./beta);		%Maximum energy in the integral
   	   de=emax/1000.;			%Energy integration step
		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=(n-nprime)/n;				%relative error in carrier density
   	if((abs(delta)) < rerr)
      	break;
     			elseif(delta < 0.) 
         	mumax=mu1;
			else 
			mumin=mu1;
   	end;
end;

			if j >= 25
   			'check convergence!'
			end;
      
      mu(nstep,k)=mu1;
      
end
%hold on;
ncarrier(nstep)=n;
   end
   
 for i=1:nstepmax
     for j=1:(nplotpoints+1)
         if mu(i,j) <= muminscale
             mu(i,j)=muminscale;
         end
     end
 end
 
%**************************************************************************
%******** plot figure **********
%**************************************************************************
%ttl=sprintf('Fig2.11, m* = %4.2fxm_0',m1);
ttl=['\rmFig2.11, \itm\rm* = ',num2str(m1,'%4.2f'),'\times\itm\rm_0'];
figure(1);
surf(kelvinscale,log10(ndensity),mu);
axis([0 maxKelvin log10(nmin) log10(nmax) muminscale Emax]);
colorbar;
colormap(jet);%jet copper bone gray winter spring
shading interp;
grid on;
xlabel('Temperature, \itT\rm (K)');
ylabel('log_{10}(\itn\rm_0) (cm^{-3})');
zlabel('Chemical potential,  \mu (meV)');
title(ttl);


