%Fig7_8b.m
%Essential Electron Transport for Device Physics
%Zero temperature Lindhard dielectric funtion in GaAs with optical phonons;
%dispersion relation 
%and log10 loss-function spectral intensity

clear all;
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(jet);
xmax=2.9;                %max value of scattered wave vector, q/kf
ymax=2.5;                %max value of energy loss, omega/Ef
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)
hbar=['\fontname{MT Extra}h\fontname{Arial}'];
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

x1=[0.001:.01:xmax+0.01];     %{0.01 6.01}, {0.1 60.01}
y1=[0.001:.01:ymax+0.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
    end
end

figure(1);
hold on;
imagesc(x1,real(y1)*Ef,log10(abs(f1)));%plot color scale image
colorbar;%log10 color scale
xlabel('Wave number, \itq\rm (\itk\rm_F)');
yttl=['Energy loss, $\hbar\omega$ (meV)'];
ylabel(yttl,'Interpreter','latex');


axis([0 xmax 0 ymax*Ef]);%([0 6 0 6*Ef])
ttl2=['\rmFig7.8b, 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)];
title(ttl2);
hold off