%Fig6_1.m
%Essential Electron Transport for Device Physics
%GaAs with LO phonons;
%plotting parameters + fontsizes
clear all;
clf;

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

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

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

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
e0=einf*(wLO/wTO)^2;        %low 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=.035;                  %Energy broadening (GAMMA / Ef) 0.01
x=0.2;                       %Wave vector (k/kf) 0.1
zeta=zeta/x^3;
a4=2*x;
%[ start:increment:stop];
y=[0:.001:1.1];             %energy loss OMEGA/Ef
ygamma=y+(i*gamma);        %energy loss OMEGA/Ef with small broadening
% 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);
% epsi=((y./(y-i*gamma)).*a3)+(einf*((y.^2-wLO2)./(y.^2-wTO2)));%with LO phonons 
%epsi=((y./(y-i*gamma)).*a3)+einf;                              %without phonons
epsi=(einf*(((y.*ygamma)-wLO2)./((y.*ygamma)-wTO2))); %with LO phonons only

%**************************************************************************

ttl=['\rmFig6.1, \itm\rm^*_e = ',num2str(ms),...
    '\times\itm\rm_0, \gamma = ',num2str(gamma*Ef,'%3.1f'),' meV'];
figure(1);                  %plot real and imaginary epsilon
plot(y*Ef,real(epsi),'b','LineWidth',LW);
hold on;
plot(y*Ef,imag(epsi),'r','LineWidth',LW);
grid on

title(ttl);
x1=wLO/10;
yscale=ylim;
y1=yscale(2)/2;
y2=y1*1.2;
y3=y1*1.4;
y4=y1*1.6;
text(x1,y1,[hbar,'\omega_{LO} = ',num2str(wLO),' meV ']);
text(x1,y2,[hbar,'\omega_{TO} = ',num2str(wTO),' meV']);
text(x1,y3,['\epsilon_{\infty} = ',num2str(e0,'%3.1f')]);
text(x1,y4,['\epsilon_{\infty} = ',num2str(einf,'%3.1f')]);
xttl=['Energy loss, ',hbar,'\omega (meV)'];
xlabel(xttl);
ylabel('Dielectric function, \epsilon(\omega)');
axis([0 max(y*Ef) 1.1*min(real(epsi)) 1.1*max(imag(epsi))]);
hold off;

figure(2);                  %plot loss function
plot(y*Ef,-imag(1./epsi),'r','LineWidth',LW);
grid on
title(ttl);
xlabel(xttl);
ylabel('Loss function, Im(-1/\epsilon(\omega))');
x1=wLO/10;
yscale=ylim;
y1=yscale(2)/2;
y2=y1*1.2;
y3=y1*1.4;
y4=y1*1.6;
text(x1,y1,[hbar,'\omega_{LO} = ',num2str(wLO),' meV ']);
text(x1,y2,[hbar,'\omega_{TO} = ',num2str(wTO),' meV']);
text(x1,y3,['\epsilon_{\infty} = ',num2str(e0,'%3.1f')]);
text(x1,y4,['\epsilon_{\infty} = ',num2str(einf,'%3.1f')]);
y1=yscale(2)/2;
axis([0 max(y*Ef) 0 1.1*max(-imag(1./epsi))]);
hold off;

