%Fig_5_3.m
%optical Lorentz dielectric constant and optical reflectivity;
clear all;
clf;
%plotting parameters + fontsizes
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);

eye=complex(0,1);       %square root of minus one
m0=9.1095e-31;          %bare electron mass [kg]
hbar=1.05459e-34;       %Plank constant hbar [J s]
echarge=1.60219e-19;    %electron charge [C]
clight=2.997924e8;      %speed of light [m s^-1]
epsilon0=8.854187e-12;  %permittivity of free space [F m^-1]

n0=12.75e28;            %electron density [m^-3] for Si, 1e25 is 1e19cm^-3
wp2=n0*echarge^2/(epsilon0*m0);
wp=sqrt(wp2);           %plasma frequency [rad s^-1]
% w0=1.00*wp/3.317;       %resonant frequency in terms of wp, set to zero for no atom restoring force 0.80
% w02=w0^2;
Ep=hbar*wp/echarge;     %plasma energy [eV]
Ep2=Ep^2;
% E0=hbar*w0/echarge;     
E0 = 4.0;               %E0 in eV for Si
E02=E0^2;

%gamma=E0/8;             %broadening [eV] E0/8 Ep*8e-2 1.81 meV
gamma = 0.5;            %gamma in eV for Si

npoints=1000;           %number of points in plot
y1min=0.001*Ep;         %minimum energy in plot
y1max=1.5094*Ep;        %maximum energy in plot
y1=linspace(y1min,y1max,npoints);      %energy [eV]
y2=y1+(eye*gamma);
%*************************************************************************
% epsilon
%*************************************************************************
eps=1-(Ep2./((y1.*y1-E02)+(eye*y1*gamma))); %complex permittivity of dielectric
x  =1./((y1.*y1-E02)+(eye*y1*gamma));       % position x
%*************************************************************************
ttl=['\it{n}\rm_0 = ',num2str(n0/1e28,'%5.1f'),...
    '\times10^{28} m^{-3}, \it{E}\rm_{0} = ',...
    num2str(E0,'%5.2f'),' eV, \it{E}\rm_{p} = ',...
    num2str(Ep,'%5.2f'),' eV, \gamma = ',num2str(gamma,'%5.2f'),' eV'];
%*************************************************************************
figure(1);
plot(y1,real(eps),'b','LineWidth',LW);%y1*1000 is meV
hold on;
grid on;
plot(y1,imag(eps),'r','LineWidth',LW);
plot(E0,0,'*');
plot(Ep,0,'o');
axis([y1min y1max 1.1*min(real(eps)) 1.1*max(imag(eps))]);
title(ttl);
xlabel('Photon energy, \it{E}\rm (eV)');
ylabel('\epsilon_r(\it{E}\rm)');
hold off;

%*************************************************************************
figure(2);
displacement=2*pi*abs(x);
phase=atan2(-imag(x),-real(x));
ttl4=['|\it{x}\rm(\omega)| (blue), \theta (red)'];
plot(y1,displacement,'b','LineWidth',LW);
hold on;
grid on;
plot(y1,phase,'r','LineWidth',LW);
title(ttl);
xlabel('Photon energy, \it{E}\rm (eV)');
ylabel(ttl4);
hold off;
axis([y1min y1max 0 1.1*pi]);