%Fig7_9.m
%Essential Electron Transport for Device Physics
%Zero temperature Lindhard dielectric funtion in GaAs with optical phonons;
%showing truncated parabola of integration 
%for injected electron energy Emax(meV)
%and log10 loss-function spectral intensity

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}'];

colormap(jet);
Emax=200;                %energy of injected electron (meV) 300
n1=1.e18;                %electron carrier concentration (cm^-3)
n=n1*1.e6;               %convert to m^-3
m0=9.1095e-31;           %bare electron mass
ms=0.07;                 %effective electron mass in conduction band
wLO=36.3;                %longitudinal optic phonon energy (meV)
wTO=33.3;                %transverse optic phonon energy (meV)
einf=11.1;               %high frequency dielectric constant
hb=1.05459e-34;          %Plank constant (J s)

e=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*ms);%Fermi energy (eV)
Ef=(Ef*1e3)/e;           %Fermi energy (meV)
wLO2=(wLO/Ef)^2;
wTO2=(wTO/Ef)^2;
zeta=ms/(pi*kf*a0);      %Lindhard function prefactor
gamma=.01;               %Energy broadening (GAMMA / Ef) 0.01, 0.025,0.055

npoints=600;
dE=Emax/npoints;
y1min=0.001;
y1max=Emax/Ef;
x1=[0.001:.01:6.01];     %{0.01 6.01}, {0.1 60.01}
y1=[0.001:.01:6.01]+i*gamma;
%y1=[0.001:dE:y1max]+i*gamma;
zeta=ms/(pi*kf*a0);
for i1=1:length(x1)
    x=x1(i1)+0.001;
    for j=1:length(y1);    
       y=y1(j)+0.001;
zeta1=zeta/x^3;
a4=2*x;
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-i*gamma)).*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

%Plot truncated parabola of integration
Etruncate=Emax-Ef;  %truncation energy (meV)
Emax1=Emax/Ef;
Etruncate1=Etruncate/Ef;
for i=1:601
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

figure(1);
hold on;
imagesc(x1,real(y1)*Ef,log10(abs(f1)));%plot color scale image
colorbar;
xlabel('Wavevector, \itq\rm (\itk\rm_F)');
yttl=['Energy loss, $\hbar\omega$ (meV)'];
ylabel(yttl,'Interpreter','latex');
axis([0 4 0 4*Ef]);%([0 6 0 6*Ef])
ttl2=['\rmFig7.9, GaAs, log_{10}, \itT\rm=0 K, \itm\rm^*_e=',num2str(ms),...
'\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),...
', \itE\rm_{max}=',num2str(Emax1*Ef),' meV'];
title(ttl2);
plot(xplot,Ef*yplot,'r','LineWidth',2); %plot truncated parabola
hold off