%Fig_6_5.m
%
%Drude dielectric
%
clear all;
clf;
t=cputime;
eye=complex(0.,1.); %square root of -1
%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}'];

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]
sigma0      =1.0*5.9e7;    %DC conductity of Cu at room temperature [S m^-1]

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 1000e6 2*pi*1000e12
omegamax=2*pi*10000e12;    %maximum frequency in plot

npoints=100000;
omega  =linspace(omegamin,omegamax,npoints);	%frequency [rad s^-1]
omega2=omega.^2;
omega2tau2=omega2*tau^2;

epsilon = 1-(sigma0*tau./(epsilon0*(1+omega2tau2)));
a=epsilon.*omega2/(clight^2);

epsilon=epsilon+(eye*sigma0./(epsilon0*omega.*(1+omega2tau2)));

b=omega2.*sigma0./((epsilon0*omega.*(1+omega2tau2))*(clight^2));
a(1)
b(1)

theta=atan(b(1)/a(1))
invdeltax=sqrt(a(1)^2+b(1)^2)
sin(theta/2)
deltax=abs((sin(theta/2))*(sqrt(1/(invdeltax))))
kx        =sqrt(omega2.*epsilon)/clight;

ttl=['\sigma_{0} = ',num2str(sigma0,'%3.2e'),...
    ' S m^{-1}, \it{n}\rm_0 = ',...
num2str(nelectrons,'%3.2e'),' m^{-3}, \tau = ',...
num2str(tau*1e15,'%2.0f'),' fs, \it{E}\rm_p = ',...
num2str(omegapeV,'%2.1f'),' eV'];

figure(1);
plot(log10(real(omega./(2e9*pi))),log10(imag(kx)),...
    'r','LineWidth',LW);  %plot imag part
hold on;
%grid on;
plot(log10(real(omega./(2e9*pi))),log10(real(kx)),...
    'b','LineWidth',LW);  %plot real part 

xlabel('Frequency, log_{10}(\it{f}\rm) (GHz)');
ylabel('log_{10}(\it{k}\rm(\omega)) real(blue) imag(red) (m^{-1})');
axis([log10(omegamin/(2e9*pi)),log10(omegamax/(2e9*pi)),...
    min(log10(imag(kx)))*0.95,max(log10(real(kx)))*1.05])
title(ttl);
hold off

figure(2);
plot(log10(real(omega./(2e9*pi))),log10(1./imag(kx)),...
    'r','LineWidth',LW);  %plot imag part
hold on;
%grid on;
plot(log10(real(omega./(2e9*pi))),log10(1./real(kx)),...
    'b','LineWidth',LW);  %plot real part 
xlabel('Frequency, log_{10}(\it{f}\rm) (GHz)');
ylabel('log_{10}(1/\it{k}\rm(\omega)) real(blue) imag(red) (m)');
axis([log10(omegamin/(2e9*pi)),log10(omegamax/(2e9*pi)),...
    min(log10(1./real(kx)))*1.05,max(log10(1./imag(kx)))*0.95])
title(ttl);
hold off