cef(z,N);

% Relative Permittivity of Gold calculated originally coded by Collin
% Meierbachtol (C)2009 based on the Brendel-Bormann model described by 
% A. D. Rakic, et. al., App. Opt., vol. 37, no. 22, 1998.

% *** Requires the cef (complex error function) file cef.m in order to run.

function w = cef(z,N);

%Visualization Parameters (USER INPUT)
lmin = 1e-7; %Minimum wavelength in visualization
lmax = 1e-6; %Maximum wavelength in visualization
dl = 5e-9; %Visualization wavelength incrament
N = 1000; %Number of iterations when calculating complex error function

%Inherent Physical constants
c = 2.997e8;
h = 6.626e-34;
hbar = h/(2*pi);
q = 1.602e-19;
i = complex(0,1);

for lambda = lmin:dl:lmax;
    omega = 2*pi*c./lambda;

    %Numerical constants of metal (extracted from cited article)
    wp = 9.03;
    f0 = 0.77;
    fj = [0.054 0.050 0.312 0.719 1.648];
    GAMMA0 = 0.050;
    GAMMAj = [0.074 0.035 0.083 0.125 0.179];
    wj = [0.218 2.885 4.069 6.137 27.97];
    sigmaj = [0.742 0.349 0.830 1.246 1.795];

    %Conversion of constants and calculations
    wp = q*wp/hbar;
    GAMMA0 = (q/hbar)*GAMMA0;
    GAMMAj = (q/hbar).*GAMMAj;
    wj = (q/hbar).*wj;
    sigmaj = (q/hbar).*sigmaj;
    OMEGA_P = sqrt(f0).*wp;
    
    %First component of relative permittivity
    epsf = 1 - OMEGA_P^2./(omega.*(omega + i*GAMMA0));

    %Second (summation) component of relative permittivity
    aj = (omega^2 + i*omega.*GAMMAj).^(0.5);
    zplus = (aj+wj)./(sqrt(2).*sigmaj);
    zminus = (aj-wj)./(sqrt(2).*sigmaj);
    epsb = (i*sqrt(pi).*fj.*wp^2./(2*sqrt(2).*aj.*sigmaj)).*(cef(zplus,N)+cef(zminus,N));
    
    %Total permittivity calculation
    eps = epsf + sum(epsb);
    
    %Visualzation over optical wavelengths (compare to Figure 2 of cited article) 
    figure(2)
    loglog(lambda,abs(real(eps)),'o',lambda,abs(imag(eps)),'+')
    hold on
    axis([lmin lmax 0.1 2000])
    xlabel('Wavelength [m]')
    ylabel('|\epsilon_{1r}| & \epsilon_{2r}')
    title('Relative Permittivity of Gold \epsilon_{r} = \epsilon_{1r} + j\epsilon_{2r}')
    legend
end
hold off