%Fig7_11b.m
%Essential Electron Transport for Device Physics
%Zero temperature Lindhard dielectric funtion in InAs with optical phonons;
%elastic scattering rate as function of electron 
%injection energy (meV)

clear;
clf;
FS      = 12;	%label fontsize 18
FSN     = 12;	%number fontsize 16
LW      = 1;	%linewidth

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times');
set(0,'DefaultAxesFontSize', FSN);

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times');
set(0,'DefaultTextFontSize', FSN);
%hbar=['\fontname{MT Extra}h\fontname{Times}'];

Emax=400;                   %maximum energy of injected electron (meV) 300
npoints=40;                 %number of points in plot
ntheta=40;                  %number of increments dth in angular integral
n1=1.e17;                   %electron carrier concentration (cm^-3)
n=n1*1.e6;                  %convert to m^-3
m0=9.1095e-31;              %bare electron mass
mass=0.021;                 %effective electron mass in conduction band
wLO=30.2;                   %longitudinal optic phonon energy (meV)
wTO=27.1;                   %transverse optic phonon energy (meV)
einf=11.8;                  %high frequency dielectric constant
hb=1.05459e-34;             %Plank constant (J s)

echarge=1.60219e-19;        %electron chanrge (C)
a0=0.529177e-8;             %Bohr radius (cm)
kf=(3.0*(pi^2)*n)^(1/3);      %Fermi wave vector (m^-1)
Ef=((hb*kf)^2)/(2*m0*mass); %Fermi energy (eV)
Ef=(Ef*1e3)/echarge;        %Fermi energy (meV)

Estep=(Emax-Ef)/npoints;    %electron increment in energy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	ak3=(kf^3);
	a=einf*((wLO/wTO)^2);   %zero frequency contribution from lattice
	const=(2.0*mass*pi*n1*3.0)/(ak3*137.036*a0*1.0e8);
	rs=((mass*.33172)/a0)*((3.0/(4.0*pi*n1))^0.3333333333333);
	ei=(Ef-Estep)/Ef;
	dth=pi/ntheta;
%
%
	for i=1:1:npoints+1
	ei=ei+(Estep/Ef);   %increment electron energy normalized to Ef
	ak=ei^-1.5;
	x1=2.0*sqrt(ei);
	aint1=0.0;          %set angular integral to zero
	theta=0.0;          %set angular integral to zero
	for j=1:1:(ntheta/2)
	theta=theta+dth;
    %********************************************************
    x=x1*sin(theta);
	b=x+(1.-((x^2)/4.))*log(abs((x+2.0)/(x-2.0)));
	g=(tan(theta))*(((sin(theta))*(a+(rs/(x^3))*b))^2);
	gchi=1./g;
    %********************************************************
    aint1=aint1+(dth*gchi);
    end
	y2(i)=const*ak*aint1;
	x2(i)=Ef*ei;
    end
figure(1);
plot(x2,y2*1e12);
xlabel('Energy, \itE_k\rm (meV)');
ylabel('Elastic scattering rate, 1/ \tau_{el} (ps^{-1})');
axis([0 Emax 0 9]);
ttl2=['\rmFig.7.11b, InAs, \itT\rm=0 K, \itm\rm^*_e=',num2str(mass),'\timesm_0, \itn\rm_0=',...
    num2str(n1/1e18),'\times10^{18} cm^{-3}, \itE\rm_F=',...
    num2str(Ef,'%5.1f'),' meV'];
title(ttl2);
return