%Fig_6_3.m
%real (blue) and imaginary (red) permittivity
%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}'];


echarge    =1.60219e-19;  %electron charge [C]
epsilon0   =8.8541878e-12;%permittivity of free-space [F 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*1000e12;
omegamax=2*pi*3000e12;

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

omega2tau2=omega.^2*tau^2;

epsilon1  =1-(omegap2tau2./(1+omega2tau2));%real permittivity
epsilon2  =(omegap2tau)./(omega.*(1+omega2tau2));%imaginary permittivity
%**************************************************************************
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, \omega_p = ',...
num2str(omegap/(2e12*pi),'%2.1f'),' THz'];
%**************************************************************************
omega=omega+delta;
omega2tau2=omega.^2*tau^2;
epsilon=1-((omega./(omega-delta)).*(omegap2tau2./(1+omega2tau2)));
%**************************************************************************
figure(1);
subplot(2,1,1);plot(real(omega./(2e12*pi)),(real(epsilon)),...
    'b','LineWidth',LW);        %plot real part 
hold on;
%grid on;
plot(real(omega./(2e12*pi)),(imag(epsilon)),...
    'r','LineWidth',LW);        %plot imag part
xlabel('Frequency, \it{f}\rm (THz)');
ylabel('\epsilon_r(\it{f}\rm)');
axis([omegamin/(2e12*pi),omegamax/(2e12*pi),...
    min(epsilon1)*1.1,max(epsilon1)*6.1])
title(ttl);
subplot(2,1,2),plot(real(omega./(2e12*pi)),...
    -imag(1./epsilon),'k','LineWidth',LW);    %plot loss function
xlabel('Frequency, \it{f}\rm (THz)');
ylabel('S(\it{f}\rm)');
axis([omegamin/(2e12*pi),omegamax/(2e12*pi),0,1.1*max(-imag(1./epsilon))]);
hold off