%Fig_6_8.m
%
%Drude dielectric
%
clear all;
clf;
t=cputime;

%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)
hbar=['\fontname{MT Extra}h\fontname{Times New Roman}'];

eye         =complex(0.,1.); %square root of -1
clight      =2.99792458e8; %speed of light in vacuum [m s^-1]
echarge     =1.60219e-19;  %electron charge [C]
epsilon0    =8.8541878e-12;%permittivity of free-space [F m^-1]
mu0         =4*pi*1e-7;    %periability of free space [H m^-1]
m0          =9.10956e-31;  %bare electron mass [kg]
meff        =1.01*m0;      %effective electron mass in Cu [kg]
nelectrons  =8.46e28;      %electron density in Cu [m^-3]
sigma00     =1.0*5.9e7;    %DC conductity of Cu at room temperature [S m^-1]

tscale=[0.01, 0.1, 1];%tau scale
for ii=1:1:3
%DC conductity of Cu at room temperature [S m^-1]
sigma0      =sigma00*tscale(ii);

tau         =sigma0*meff/(nelectrons*(echarge^2));%scattering time [s]
delta       =eye*1.0/tau;   %broadening 

omegap2     =nelectrons*(echarge^2)/(meff*epsilon0);%plasma frequency squared
omegap      =sqrt(omegap2);%plasma frequency [rad s^-1]
omegapeV    =sqrt(omegap2)*6.58211889e-16;%plasma energy [eV]
omegap2tau  =omegap2*tau;
omegap2tau2 =omegap2*tau^2;

omegamin=2*pi*1000e6;     %minimum frequency in plot 2*pi*1000e6 2*pi*1000e12 [rad s^-1]
omegamax=2*pi*10000e12;   %maximum frequency in plot [rad s^-1] 10000e12

npoints=1000000;             %100000
omega  =linspace(omegamin,omegamax,npoints);	%frequency [rad s^-1]

omega2      =omega.^2;
omega2tau2  =omega2*tau^2;
klight      =omega/clight;

epsilon     = 1-(sigma0*tau./(epsilon0*(1+omega2tau2)));
epsilon     =epsilon+(eye*sigma0./(epsilon0*omega.*(1+omega2tau2)));
x           =tau./((omega.*(eye+omega.*tau)));% position x
if ii==1
    xconst=max(abs(x));%normalization constant
end
x           =x./xconst;

kx        =sqrt(omega2.*epsilon)/clight;
kkx       =sqrt(omega2-omegap2)/clight;

n_index     =clight*kx./omega;%complex refractive index
nr          =real(n_index);%real part of refractive index
nrimag      =imag(n_index);%imaginary part of refractive index

aclight     =omega./klight;

vgroup=zeros(npoints,1);
for i=1:1:npoints-1 
    vgroup(i)      =(omega(i+1)-omega(i))/(real(kx(i+1))-real(kx(i)));
end
vgroup(npoints)=vgroup(npoints-1);

ttl=['\sigma_{0} = ',...
    num2str(sigma0,'%3.2e'),' S m^{-1}, \it{n}\rm_0 = ',...
num2str(nelectrons/1e28,'%3.2f'),'\times10^{28} m^{-3}, \tau = ',...
num2str(tau*1e15,'%2.0f'),' fs, \it{E}\rm_p = ',...
num2str(omegapeV,'%2.1f'),' eV'];
%*************************************************************************
figure(1)
%reflectivity from vacuum at normal incidence
reflect=(abs((n_index-1)./(n_index+1))).^2;  
plot(log10(omega./(2e9*pi)),reflect,'b','LineWidth',LW);
hold on;
%grid on;
% plot(E0,0,'*');
% plot(Ep,0,'o');
title(ttl);
xlabel('Frequency, log_{10}(\it{f}\rm) (GHz)');
ylabel('\it{R}\rm');
%hold off;
axis([log10(omegamin/(2e9*pi)),log10(omegamax/(2e9*pi)), 0, 1.1]);

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

figure(2)
ttl10a=['Amplitude, log_{10}|\it{x}\rm_0(\omega)|'];
ttl10b=['Phase, \phi (radians/\pi)'];
displacement=abs(x);
phase=atan2(-imag(x),-real(x));

subplot(2,1,1);plot(log10(omega./(2e9*pi)),...
    log10(displacement),'b','LineWidth',LW);
hold on;
%grid on;
title(ttl);
ylabel(ttl10a);
axis([log10(omegamin/(2e9*pi)),log10(omegamax/(2e9*pi)), -8, 2]);
subplot(2,1,2);plot(log10(omega./(2e9*pi)),phase/pi,'r','LineWidth',LW);
hold on;
%grid on;

xlabel('Frequency, log_{10}(\it{f}\rm) (GHz)');
ylabel(ttl10b);
axis([log10(omegamin/(2e9*pi)),log10(omegamax/(2e9*pi)), 0, 1.1]);

end
hold off