%Fig7_11a.m
%Essential Electron Transport for Device Physics
%Zero temperature Lindhard dielectric funtion in InAs with optical phonons;
%calculates inelastic 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);

%colormap (copper);
Emax=400;                   %maximum energy of injected electron (meV) 300
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
eye=complex(0,1);           %square root of minus one
hb=1.05459e-34;             %Plank constant (J s)
%hbar=['\fontname{MT Extra}h\fontname{Arial}'];
echarge=1.60219e-19;        %electron chanrge (C)
a0=0.529177e-10;            %Bohr radius (m)
kf=(3*(pi^2)*n)^(1/3);      %Fermi wave vector (m^-1)
kf1=kf*1e-2;                %Fermi wave vector (cm^-1)
Ef=((hb*kf)^2)/(2*m0*mass); %Fermi energy (eV)
Ef=(Ef*1e3)/echarge;        %Fermi energy (meV)
wLO2=(wLO/Ef)^2;
wTO2=(wTO/Ef)^2;
zeta=mass/(pi*kf*a0);       %Lindhard function prefactor
gamma=3.00;                 %Energy broadening (meV)3, 0.52
gamma1=(eye*gamma)/Ef;      %normalized broadening

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%setup y1=hw/Ef and x1=q/kf arrays
dE=0.35;                    %energy increment (meV)
npoints=int16(Emax/dE);     %number of integration points
dx1=2*(sqrt(Emax/Ef))/(double(npoints)); %increment in x1 
dy1=dE/Ef;                  %increment in y1 

x1=[dx1:dx1:(double(npoints))*dx1];
y1=[dy1:dy1:(double(npoints))*dy1]+gamma1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i1=1:npoints
    x=x1(i1);
    zeta1=zeta/x^3;
    a4=2*x;
   for j=1:length(y1);    
    y=y1(j);
a1=(1-1/4*(x-y/x).^2).*log((y-x*(x+2))./(y-x*(x-2)));
a2=(1-1/4*(x+y/x).^2).*log((-y-x*(x+2))./(-y-x*(x-2)));
a3=zeta1*(a4+a1+a2);
%epsi=((y./(y-i*gamma)).*a3)+einf;          %Electrons only and fix a3    
%epsi=a3+einf;                              %Electrons only    
%epsi=a3+(einf*((y.^2-wLO2)./(y.^2-wTO2))); %Electrons+LO
epsi=((y./(y-gamma1)).*a3)+(einf*((y.^2-wLO2)./(y.^2-wTO2))); %Electrons+LO and fix a3
f1(j,i1)=-imag(1./epsi);   %loss function
f1(j,i1)=f1(j,i1)/x1(i1);  %loss function weighted by 1/q
    end
end

%truncated parabola of integration parameters
Etruncate=Emax-Ef;          %truncation energy (meV)
Emax1=Emax/Ef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up injection energy array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Emaxtmp=Emax1;
nmax=100;           %plot 100 values of electron injection enegy
efactor=Emax/Ef;
for jj=1:1:nmax; 
energyplot(jj)=(double(jj)*Ef*efactor)/(nmax);
Emax1=Emaxtmp*(double(jj))/nmax;
Etruncate1=Emax1-1;  %truncation

for i=1:npoints;
xplot(i)=x1(i);
yplot(i)=Emax1-((x1(i)-sqrt(Emax1))^2); %parabola of integration
    if yplot(i) > Etruncate1            %tuncate parabola
    yplot(i)=Etruncate1;                %truncate parabola
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%set up meshgrid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X1,Y1]=meshgrid(x1,y1);
sumindicator=Y1<=ones(length(y1),1)*yplot; %find valaues of sum indicator for each column
Integral(jj)=sum(sum(f1(sumindicator)))/(sqrt(Emax1)); %this is the integral for current Emax
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
prefactor=(2.0*mass*dx1*dy1*Ef*0.001*echarge)/(pi*a0*kf*hb);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
plot(energyplot,Integral*prefactor/1e12);
xlabel('Energy, \itE_k\rm (meV)');
ylabel('Inelastic scattering rate, 1/ \tau_{in} (ps^{-1})');
ttl2=['\rmFig7.11a, InAs, \itT\rm=0 K, \itm\rm^*_e=',num2str(mass),...
'\times\itm\rm_0, \itn\rm_0=',num2str(n1/1e18),'\times10^{18} cm^{-3}, \itE\rm_F=',...
num2str(Ef,'%5.1f'),' meV, \gamma=',num2str(gamma),...
' meV, \itE\rm_{max}=',num2str(Emax1*Ef),' meV'];
title(ttl2);
return