%Fig5_14a.m
%Essential Electron Transport for Device Physics
%elastic ionized impurity scattering rate using TF screeing
%m=electron mass, n=number of charge impurities
function TFelasticScattering
clear;
clf;
  
m0=9.109382e-31;				%bare electron mass (kg)
echarge=1.6021764e-19;		%electron charge (C)
hbar=1.05457159e-34;			%Planck's constant (J s)
epsilon0=8.8541878e-12;		%permittivity of free-space (F m-1) 
ezero=13.2;						%low-frequency relative dielectric constant
m1=0.070;						%effective electron mass
emax=300.0d-3;					%maximum particle energy (eV)
n=1.0d17;						%carrier density, n (cm-3)

npoints = 200;					%number of points in energy plot
 
 for k=1:1
   n=n*10;						%plot two curves n=10^17 cm-3 and n=10^18 cm-3

 kF1=(3*(pi^2)*n)^(1/3);	%Fermi wave vector (cm-1)
 kF1=kF1*1.e2;					%convert to units of m-1
    kF2=kF1^2;
    kF3=kF1^3;

 qTF2=kF1*m0*m1*(echarge^2)/(epsilon0*ezero*(pi^2)*(hbar^2));	% Thomas-Fermi screening number squared
 eF=((hbar*kF1)^2)/(2.0*m0*m1*echarge);								% Fermi energy in eV

 const=(pi*n*1.0e6*(echarge^4)*m0*m1);
 const=const/(((4.*pi*epsilon0)^2)*(ezero^2)*(hbar^3));
 
 emin=eF;						%initial particle energy
 demax=(emax-emin)/200.0;	%increment in particle energy (eV)
 dth=pi/180.0;					%increment in scattering angle theta

 ei=emin;

 for i=1:npoints				% the main loop
   ei=ei+demax;
   ak=(sqrt(2.*m0*m1*ei*echarge))/hbar;
      ak3=ak^3;
   aint1=0.0;
   theta=0.0;

  for j=1:180
     theta=theta+dth;
     aint1=aint1+(dth*gchi1(theta,ak,qTF2));		%call function aint1
  end

	y1(i)=const*aint1/ak3;
   x1(i)=ei*1.d3;
 end
%plot figure
ttl=['\rmFig5.14a, \itn\rm_{0}=',num2str(n,'%4.1e'),...
    ' cm^{-3}, \itm\rm^*_e=',num2str(m1,'%5.2f'),...
    '\times\itm\rm_0, \epsilon_{r0}=',num2str(ezero,'%5.1f')]; 
figure(1);
hold on;
plot(x1,y1*10^-12,'r');%inverse ps
grid on
axis([0,emax*1.d3,0,16]);
end;
xlabel('Energy, \itE\rm (meV)');
ylabel('Elastic scattering rate, 1/\tau_{el} (ps^{-1})');

title(ttl);
hold off;

function gchi1 = gchi1(theta,k,qTF2)
 x=2.0*k*sin(theta/2.0);
 x2=x^2;	
 g=(tan(theta/2.0))*((sin(theta/2.))^2)*((1.0+(qTF2/x2))^2);
 gchi1=1.0/g;
return