%FigG_2
%Zero temperature Lindhard dielectric funtion in GaAs with or without
%longitudinal optical phonons;
clear;
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_txt=['\fontname{MT Extra}h\fontname{Arial}'];

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);
gamma=.01;                  %Energy broadening (GAMMA / Ef) 0.01
x=0.1;                       %Wave vector (k/kf) 0.1
zeta=zeta/x^3;
a4=2*x;
%[ start:increment:stop];
y=[0:.01:2.0]+i*gamma;        %energy loss OMEGA/Ef
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=zeta*(a4+a1+a2);
%with LO phonons 
%epsi=((y./(y-i*gamma)).*a3)+(einf*((y.^2-wLO2)./(y.^2-wTO2)));
%without phonons
epsi=((y./(y-i*gamma)).*a3)+einf;                           

figure(1);                  %plot real and imaginary epsilon
plot(real(y)*Ef,real(epsi),'b');
hold on;
plot(real(y)*Ef,imag(epsi),'r');
axis([0 100 -600 800]);
ttl=['\rmFigG.2, \itn\rm_0=',num2str(n1,'%5.1e'),' cm^{-3}, \itm\rm^*_e=',num2str(ms),...
'\times\itm\rm_0, \itE\rm_F=',num2str(Ef,'%5.1f'),' meV, \itk\rm_F=',num2str(kf, '%5.1e'),...
' cm^{-1}, \itX\rm=',num2str(x),'\times\itk\rm_F, \gamma=',num2str(gamma),'\times\itE\rm_F'];
title(ttl);
x1=Ef;
yscale=ylim;
y1=yscale(2)/2;
y2=y1*1.2;
y3=y1*1.4;
y4=y1*1.6;
text(x1,y1,['\itk\rm_F = ',num2str(kf1, '%5.1e'),' cm^{-1}']);
text(x1,y2,['\itX\rm = ',num2str(x),'\times\itk\rm_F,  \gamma = ',num2str(gamma),'\times\itE\rm_F']);
%text(x1,y3,[hbar,'\omega_{LO} = ',num2str(wLO),' meV, ',hbar,'\omega_{TO} = ',num2str(wTO),' meV']);
text(x1,y4,['\epsilon_{r\infty} = ',num2str(einf)]);
xttl=['Energy loss, $\hbar\omega$ (meV)'];
xlabel(xttl,'Interpreter','latex');
ylabel('Dielectric function, \epsilon');
hold off;

figure(2);                  %plot loss function
plot(real(y)*Ef,-imag(1./epsi),'r');
title(ttl);
xlabel('Energy loss, $\hbar\omega$ (meV)', 'Interpreter', 'latex');
ylabel('Loss function, -Im/\epsilon');
hold on;                    %keep current plot in the figure;
%axis manual
axis([0 100 0 5]);
plot(real(y)*Ef,-imag(100./epsi),'k');  %plot 100x y-axis to see detail structure
yscale=ylim;
y1=yscale(2)/2;
text(x1,y1,['x100']);
hold off;