%Fig_5_4.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] 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;     
E02=E0^2;

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

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 is complex permittivity of dielectric
%*************************************************************************
eps=1-(Ep2./((y1.*y1-E02)+(eye*y1*gamma)));
%*************************************************************************
ttl=['\it{n}\rm_0 = ',num2str(n0/1e28,'%5.1f'),...
    ' \times 10^{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'];
%*************************************************************************
omega=y1*echarge/hbar;  %omega [rad s^-1]
n_index=(eps).^0.5;     %complex refractive index
%*************************************************************************
figure(1)
% reflectivity at normal incidence is a ratio of complex index squared
reflect=(abs((n_index-1)./(n_index+1))).^2;
plot(y1,reflect,'r','LineWidth',LW);
hold on;
%grid on;
plot(E0,0,'*');
plot(Ep,0,'o');
title(ttl);
xlabel('Photon energy, \it{E}\rm (eV)');
ylabel('\it{R}\rm');
hold off;
axis([y1min y1max 0 1]);

%*************************************************************************
figure(2)
ttl4=['\it{n}\rm_{index}'];
plot(y1,real(n_index),'b','LineWidth',LW);
hold on;
plot(y1,imag(n_index),'r','LineWidth',LW);
%grid on;
plot(E0,0,'*');
plot(Ep,0,'o');
title(ttl);
xlabel('Photon energy, \it{E}\rm (eV)');
ylabel(ttl4);
hold off;
axis([y1min y1max 0 1.1*max(real(n_index))]);

