%Chapt8Exercise7a
%Plot Thomas Fermi scattering rate as function of angle and compare with RPA
clear;
clf;
n=1e18;						%electron carrier density (cm-3)
n=n*1e6;					%convert to (m-3)
m0=9.10956e-31;				%bare electron mass (kg)
echarge=1.602192e-19;		%electron charge (C)
epsilon0=8.85419e-12;		%Permittivity of free space (F m-1)
hbar=1.054592e-34;			%Planck's constant (J s)
hbar3=hbar^3;

m=0.07*m0;					%effective electron mass
epsilonr0=13.2;				%relative low frequency dielectric constant
epsilon=epsilon0*epsilonr0;
kF=(3*(pi^2)*n)^(1/3)		%Fermi wave vector (m-1)

E=-0.1*echarge;
for j=1:1:2
E=E+0.2*echarge;			%Electron energy (eV)   

k=sqrt(2*m*E)/hbar;			%Electron wave vector (m-1)
k3=k^3;

%RPA model constants
   rs0=((3/(4*pi*n))^(1/3))*(m*(echarge^2)/(4*pi*epsilon0*(hbar^2)));
   xi=((32*(pi^2)/9)^(1/3))/(pi^2);
%Thomas-Fermi model constants  
qTF=sqrt(kF*m*echarge^2/(epsilon*(pi^2)*(hbar^2)));

%plot scattering rate as function of scattered angle
theta=[pi/1800:pi/1800:pi];
	q=2*k*sin(theta/2);			%scattered wave vector (m-1)
   eta=sin(theta/2);			%eta calculated using theta in radians
   deta=pi*cos(theta/2)./2/180;	%differential value of eta in degrees
   eta3=eta.^3;
   
   %RPA model
   x=q/kF;
   absx=abs(x+2)./abs(x-2);
   RPAepsilonr=epsilonr0+((rs0./x.^3)*xi.*(x+((1-((x.^2)/4)).*log(absx))));
   %Thomas-Fermi model
   TFepsilon=epsilon*(1+qTF^2./q.^2);
   
   RPArate=2*pi*m/hbar3/k3*n*(echarge^2/4/pi/epsilon0)^2.*deta./RPAepsilonr.^2./eta3;
   TFrate =2*pi*m/hbar3/k3*n*(echarge^2/4/pi)^2.*deta./TFepsilon.^2./eta3;
   
  % plot(theta*180/pi, RPArate);
  hold on;
  figure(1);
  plot(theta*180/pi, TFrate,'r');
  xlabel('angle, theta (degree)');
  ylabel('scattering rate, 1/tau (s-1 degree-1)');
  
  plot(theta*180/pi, RPArate,'b');
  xlabel('Angle, theta (degree)');
  ylabel('Elastic scattering rate, 1/tau (s-1 degree-1)');
 end 
ttl=sprintf('Chapt8Exercise7a, n=%4.1e cm^{-3}, m*=%5.2f m_0, epsilon0r=%5.1f',(n/1e6),(m/m0),epsilonr0);
title(ttl);
 hold off;

