%Fig_5_5.m
%optical Lorentz dielectric constant and optical reflectivity;
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);
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=20000;          %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 complex permittivity of dielectric
%*************************************************************************
eps=1-(Ep2./((y1.*y1-E02)+(eye*y1*gamma))); 
x  =1./((y1.*y1-E02)+(eye*y1*gamma));       % position x
x1  =1./((y1.*y1-E02)+(eye*y1*gamma*2));    %position x with 2x damping
x2  =1./((y1.*y1-E02)+(eye*y1*gamma*4));    %position x with 4x damping

%*************************************************************************
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_{p}}\rm=',...
    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
k=n_index.*omega/clight;%complex wave vector [m^-1]

vg_clight   =omega./(y1*echarge/(clight*hbar));%velocity of light
vp          =omega./real(k);%phase velocity in medium is vp=omega/Re(k)
%*************************************************************************
%******* group velocity is v_g = d omega / d k ***************************
%*************************************************************************
for i=1:npoints-1
vg(i)=(y1max*echarge/(clight*hbar*(npoints+1)))/(real(k(i+1))-real(k(i)));
end
vg(npoints)=vg(npoints-1);

%*************************************************************************
figure(1);
plot(real(k),omega./1e15,'b','LineWidth',LW);
hold on;
%grid on;
plot(imag(k),omega./1e15,'r','LineWidth',LW);
plot(omega/clight,omega./1e15,'r');
plot(0,w0*1e-15,'*');
plot(0,wp*1e-15,'o');
axis([0 1.1*max(real(k)) 0 max(omega./1e15)]);
title('Real and imaginary part of \it{k}\rm')
ylabel('Frequency, \omega (\times 10^{15} rad s^{-1})');
xlabel('Wave number, \it{k}\rm (m^{-1})');
hold off;