% This is a base for a program intended to investigate properties of
% the PML absorbing boundary condition in one-dimensional FDTD.
% Daniel Sjoberg, 2012-09-19
function refl = pml(Nt, Nz, a, tw, t0, f0, Npml, sigmaPML, R, ...
WhichProfile, DoAnimation);
% Some physical constants
c0 = 299792458; % Speed of light in vacuum
mu0 = 4*pi*1e-7; % Permeability of vacuum
eps0 = 1/c0^2/mu0; % Permittivity of vacuum
eta0 = sqrt(mu0/eps0); % Wave impedance of vacuum
% Computational variables
E = zeros(1, Nz+1); % One extra value for the E-field
H = zeros(1, Nz); % Magnetic field, only internal points
z = linspace(0, a, Nz+1); % Vector to plot the E-field against
dz = a/Nz; % Spatial discretization
dt = R*dz/c0; % Temporal discretization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set up the pulse
%
t = linspace(0, (Nt-1)*dt, Nt);
pulsefunction = exp(-(t-t0).^2/(tw/2).^2).*sin(2*pi*f0*(t-t0));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Design of the PML
%
% The fields are given on a shifted grid, therefore
% we need to compute the A and B parameters for 2*Npml-1 values.
% Note that sigmaPML = sigma_electric/epsilon0 = sigma_magnetic/mu0.
%
if WhichProfile == 0
profile = ones(1, 2*Npml-1); % Constant profile
else
profile = linspace(0, 1, 2*Npml-1); % Linear profile
end;
sigmaPMLvec = sigmaPML*profile;
% Compute parameters for a given sigma distribution
A = (1-sigmaPMLvec*dt/2)./(1+sigmaPMLvec*dt/2);
B = 1./(1+sigmaPMLvec*dt/2);
Ah = A(1:2:end); % Note length(Ah) = Npml
Bh = B(1:2:end); % Note length(Bh) = Npml
Ae = A(2:2:end); % Note length(Ae) = Npml-1
Be = B(2:2:end); % Note length(Be) = Npml-1
Epml = zeros(1, Npml); % Initial conditions
Hpml = zeros(1, Npml); % Initial conditions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The main loop of the simulation. The endpoints of the
% E-field (E(1) and Epml(end)) are never updated, thus
% they are always zero, corresponding to PEC boundary
% condition.
%
for k=1:Nt
H = H - (dt/dz)/mu0*diff(E);
EpmlExtended = [E(end), Epml];
Hpml = Ah.*Hpml - (dt/dz)/mu0*Bh.*diff(EpmlExtended);
E(2:end-1) = E(2:end-1) - (dt/dz)/eps0*diff(H);
Epml(1:end-1) = Ae.*Epml(1:end-1) - (dt/dz)/eps0*Be.*diff(Hpml);
E(end) = E(end) - (dt/dz)/eps0*(Hpml(1)-H(end));
E(1) = pulsefunction(k);
if DoAnimation == 1
pause(0.01);
plot(z, E);
end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Insert your calculation of the reflection coefficient here.
refl = 0;