%Fig_6_6.m
%Skin depth
%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]
% epsilon=zeros(1,npoints);
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}, n = ',...
num2str(nelectrons,'%3.2e'),' m^{-3}, \tau = ',...
num2str(tau*1e15,'%2.0f'),' fs, E_p = ',num2str(omegapeV,'%2.1f'),' eV'];

%analytic expression for skin depth at omegamin [m]
deltas=sqrt(2/(omegamin*sigma0*mu0))

deltakxre=real(1/kx(1))
deltakxim=imag(1/kx(1))
deltakxss=1/abs(deltakxre+eye*deltakxim/(deltakxre^2+deltakxim^2))

L0=5*deltas;%maximum length to plot field into metal 5, 15
x=linspace(0,L0,npoints);%set up x-axis
E0=1;
E=E0*exp(eye*kx(1).*x);
Ereal=real(E);
Eimag=imag(E);
Emag =abs(E);

%**********************************************************************
for i=1:npoints
    if Emag(i) >= exp(-1)
        index=i;
    end
end
deltakxs=x(index)%numerical value of position when |E| is e^-1
%**********************************************************************

Emag2=Emag.^2;

ttl2=['\it{f}\rm=',num2str(omegamin/(2e9*pi)),...
'GHz, \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, \delta_s=',...
num2str(deltas*1e9,'%2.1f'),'nm'];

figure(1)
plot(x*1e9,Ereal,'b','LineWidth',LW);  %plot real part 
hold on
plot(x*1e9,Eimag,'r','LineWidth',LW);  %plot imag part 
%grid on;
xlabel('Position, \it{x}\rm (nm)');
ylabel('\it{E}\rm/\it{E}\rm_0 real(blue) imag(red)');
axis([0, L0*1e9, -0.1,1.1])
title(ttl2);

xline=0.3678*ones(1,npoints);%the e^-1 level
figure(2)
plot(x*1e9,Emag,'b','LineWidth',LW);  %plot magnitude of E 
hold on;
plot(x*1e9,xline,':r','LineWidth',LW);  %e^-1 is 18.3nm at 100THz in Cu
%grid on;
xlabel('Position, \it{x}\rm (nm)');
ylabel('|\it{E}\rm/\it{E}\rm_0| ');
axis([0, L0*1e9, 0,1.1])
title(ttl2);
