%Fig_5_6
%gaas_LO_phonon_epsilon.m
%GaAs with LO phonons;
%plotting parameters + fontsizes
clear all;
clf;

FS      = 18;	%label fontsize 18
FSN     = 16;	%number fontsize 16
LW      = 2;	%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}'];
%**************************************************************************

eye=complex(0,1);           %square root of minus one

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

wLO2=wLO^2;
wTO2=wTO^2;

gamma=1.8;                  %Energy broadening GAMMA (meV) 1.8

%[ start:increment:stop];
y=[0:.01:55];               %energy loss OMEGA (meV)
ygamma=y+(eye*gamma);       %energy loss OMEGA with small broadening
epsi=(einf*(((y.*ygamma)-wLO2)./((y.*ygamma)-wTO2))); %with LO phonons only

%**************************************************************************
ttl=['\gamma = ',num2str(gamma,'%3.1f'),' meV'];
figure(1);                  %plot real and imaginary epsilon
plot(y,real(epsi),'b','LineWidth',LW);
hold on;
plot(y,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_{r,0 } = ',num2str(e0,'%3.1f')]);
text(x1,y4,['\epsilon_{r,\infty} = ',num2str(einf,'%3.1f')]);
xttl=['Energy, ',hbar,'\omega (meV)'];
xlabel(xttl);
ylabel('Dielectric function, \epsilon_r(\omega)');
axis([0 max(y) 1.1*min(real(epsi)) 1.1*max(imag(epsi))]);
hold off;

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