%***********************************************************************
%     3-D FDTD code with UPML absorbing boundary conditions
%***********************************************************************
%
%     Program author: Keely J. Willis, Graduate Student
%                     UW Computational Electromagnetics Laboratory
%                           Director: Susan C. Hagness
%                     Department of Electrical and Computer Engineering
%                     University of Wisconsin-Madison
%                     1415 Engineering Drive
%                     Madison, WI 53706-1691
%                     kjwillis@wisc.edu

%
%     This MATLAB M-file implements the finite-difference time-domain
%     solution of Maxwell's curl equations over a three-dimensional
%     Cartesian space lattice comprised of uniform cubic grid cells.
%     
%     The dimensions of the computational domain are 8.2 cm
%     (x-direction), 3.4 cm (y-direction), and 3.2 cm (z-direction).  
%     The grid is terminated with UPML absorbing boundary conditions.
%
%     An electric current source comprised of two collinear Jz components
%     (realizing a Hertzian dipole) excites a radially propagating wave.  
%     The current source is located in the center of the grid.  The 
%     source waveform is a differentiated Gaussian pulse given by 
%          J(t)=J0*(t-t0)*exp(-(t-t0)^2/tau^2), 
%     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc-
%     content pulse is approximately 7 GHz. The grid resolution 
%     (dx = 2 mm) was chosen to provide at least 10 samples per 
%     wavelength up through 15 GHz.
%
%
%***********************************************************************

clear

%***********************************************************************
%     Fundamental constants
%***********************************************************************

cc=2.99792458e8;
muz=4.0*pi*1.0e-7;
epsz=1.0/(cc*cc*muz);
etaz=sqrt(muz/epsz);

%***********************************************************************
%     Material parameters 
%***********************************************************************

mur=1.0;
epsr=1.0;
eta=etaz*sqrt(mur/epsr);

%***********************************************************************
%     Grid parameters
%
%     Each grid size variable name describes the number of sampled points 
%     for a particular field component in the direction of that component.
%     Additionally, the variable names indicate the region of the grid 
%     for which the dimension is relevant.  For example, ie_tot is the 
%     number of sample points of Ex along the x-axis in the total 
%     computational grid, and jh_bc is the number of sample points of Hy 
%     along the y-axis in the y-normal UPML regions.
%
%***********************************************************************

ie=41;          % Size of main grid
je=17;
ke=16;
ih=ie+1;
jh=je+1;   
kh=ke+1;   

upml=10;        % Thickness of PML boundaries
ih_bc=upml+1;
jh_bc=upml+1;
kh_bc=upml+1;

ie_tot=ie+2*upml;          % Size of total computational domain
je_tot=je+2*upml;        
ke_tot=ke+2*upml;        
ih_tot=ie_tot+1;
jh_tot=je_tot+1;          
kh_tot=ke_tot+1;          

is=round(ih_tot/2);         % Location of z-directed current source
js=round(jh_tot/2);
ks=round(ke_tot/2);

%***********************************************************************
%     Fundamental grid parameters
%***********************************************************************

delta=0.002;
dt=delta*sqrt(epsr*mur)/(2.0*cc);
nmax=100;

%***********************************************************************
%     Differentiated Gaussian pulse excitation
%***********************************************************************

rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
J0=-1.0*epsz;

%***********************************************************************
%     Initialize field arrays
%***********************************************************************

ex=zeros(ie_tot,jh_tot,kh_tot);
ey=zeros(ih_tot,je_tot,kh_tot);
ez=zeros(ih_tot,jh_tot,ke_tot);
dx=zeros(ie_tot,jh_tot,kh_tot);
dy=zeros(ih_tot,je_tot,kh_tot);
dz=zeros(ih_tot,jh_tot,ke_tot);

hx=zeros(ih_tot,je_tot,ke_tot);
hy=zeros(ie_tot,jh_tot,ke_tot);
hz=zeros(ie_tot,je_tot,kh_tot);
bx=zeros(ih_tot,je_tot,ke_tot);
by=zeros(ie_tot,jh_tot,ke_tot);
bz=zeros(ie_tot,je_tot,kh_tot);

%***********************************************************************
%     Initialize update coefficient arrays
%***********************************************************************

C1ex=zeros(size(ex));
C2ex=zeros(size(ex));
C3ex=zeros(size(ex));
C4ex=zeros(size(ex));
C5ex=zeros(size(ex));
C6ex=zeros(size(ex));

C1ey=zeros(size(ey));
C2ey=zeros(size(ey));
C3ey=zeros(size(ey));
C4ey=zeros(size(ey));
C5ey=zeros(size(ey));
C6ey=zeros(size(ey));

C1ez=zeros(size(ez));
C2ez=zeros(size(ez));
C3ez=zeros(size(ez));
C4ez=zeros(size(ez));
C5ez=zeros(size(ez));
C6ez=zeros(size(ez));

D1hx=zeros(size(hx));
D2hx=zeros(size(hx));
D3hx=zeros(size(hx));
D4hx=zeros(size(hx));
D5hx=zeros(size(hx));
D6hx=zeros(size(hx));

D1hy=zeros(size(hy));
D2hy=zeros(size(hy));
D3hy=zeros(size(hy));
D4hy=zeros(size(hy));
D5hy=zeros(size(hy));
D6hy=zeros(size(hy));

D1hz=zeros(size(hz));
D2hz=zeros(size(hz));
D3hz=zeros(size(hz));
D4hz=zeros(size(hz));
D5hz=zeros(size(hz));
D6hz=zeros(size(hz));

%***********************************************************************
%     Update coefficients, as described in Section 7.8.2.
%
%     In order to simplify the update equations used in the time-stepping
%     loop, we implement UPML update equations throughout the entire
%     grid.  In the main grid, the electric-field update coefficients of 
%     Equations 7.91a-f and the correponding magnetic field update
%     coefficients extracted from Equations 7.89 and 7.90 are simplified
%     for the main grid (free space) and calculated below.
%
%***********************************************************************

C1=1.0;
C2=dt;
C3=1.0;
C4=1.0/2.0/epsr/epsr/epsz/epsz;
C5=2.0*epsr*epsz;
C6=2.0*epsr*epsz;

D1=1.0;
D2=dt;
D3=1.0;
D4=1.0/2.0/epsr/epsz/mur/muz;
D5=2.0*epsr*epsz;
D6=2.0*epsr*epsz;

%***********************************************************************
%     Initialize main grid update coefficients
%***********************************************************************

C1ex(:,jh_bc:jh_tot-upml,:)=C1;     
C2ex(:,jh_bc:jh_tot-upml,:)=C2;
C3ex(:,:,kh_bc:kh_tot-upml)=C3;
C4ex(:,:,kh_bc:kh_tot-upml)=C4;
C5ex(ih_bc:ie_tot-upml,:,:)=C5;
C6ex(ih_bc:ie_tot-upml,:,:)=C6;

C1ey(:,:,kh_bc:kh_tot-upml)=C1;
C2ey(:,:,kh_bc:kh_tot-upml)=C2;
C3ey(ih_bc:ih_tot-upml,:,:)=C3;
C4ey(ih_bc:ih_tot-upml,:,:)=C4;
C5ey(:,jh_bc:je_tot-upml,:)=C5;
C6ey(:,jh_bc:je_tot-upml,:)=C6;

C1ez(ih_bc:ih_tot-upml,:,:)=C1;
C2ez(ih_bc:ih_tot-upml,:,:)=C2;
C3ez(:,jh_bc:jh_tot-upml,:)=C3;
C4ez(:,jh_bc:jh_tot-upml,:)=C4;
C5ez(:,:,kh_bc:ke_tot-upml)=C5;
C6ez(:,:,kh_bc:ke_tot-upml)=C6;

D1hx(:,jh_bc:je_tot-upml,:)=D1;
D2hx(:,jh_bc:je_tot-upml,:)=D2;
D3hx(:,:,kh_bc:ke_tot-upml)=D3;
D4hx(:,:,kh_bc:ke_tot-upml)=D4;
D5hx(ih_bc:ih_tot-upml,:,:)=D5;
D6hx(ih_bc:ih_tot-upml,:,:)=D6;

D1hy(:,:,kh_bc:ke_tot-upml)=D1;
D2hy(:,:,kh_bc:ke_tot-upml)=D2;
D3hy(ih_bc:ie_tot-upml,:,:)=D3;
D4hy(ih_bc:ie_tot-upml,:,:)=D4;
D5hy(:,jh_bc:jh_tot-upml,:)=D5;
D6hy(:,jh_bc:jh_tot-upml,:)=D6;

D1hz(ih_bc:ie_tot-upml,:,:)=D1;
D2hz(ih_bc:ie_tot-upml,:,:)=D2;
D3hz(:,jh_bc:je_tot-upml,:)=D3;
D4hz(:,jh_bc:je_tot-upml,:)=D4;
D5hz(:,:,kh_bc:kh_tot-upml)=D5;
D6hz(:,:,kh_bc:kh_tot-upml)=D6;

%***********************************************************************
%     Fill in PML regions
% 
%     PML theory describes a continuous grading of the material properties
%     over the PML region.  In the FDTD grid it is necessary to discretize
%     the grading by averaging the material properties over a grid cell 
%     width centered on each field component.  As an example of the 
%     implementation of this averaging, we take the integral of the 
%     continuous sigma(x) in the PML region
%   
%         sigma_i = integral(sigma(x))/delta
%   
%     where the integral is over a single grid cell width in x, and is 
%     bounded by x1 and x2.  Applying this to the polynomial grading of 
%     Equation 7.60a produces
%
%         sigma_i = (x2^(m+1)-x1^(m+1))*sigmam/(delta*(m+1)*d^m)
%
%     where sigmam is the maximum value of sigma as described by Equation 
%     7.62. 
%         
%     The definitions of x1 and x2 depend on the position of the field 
%     component within the grid cell.  We have either
%
%         x1 = (i-0.5)*delta,  x2 = (i+0.5)*delta
%  
%     or
%  
%         x1 = (i)*delta,      x2 = (i+1)*delta
%
%     where i varies over the PML region.
%  
%***********************************************************************

rmax=exp(-16);  %desired reflection error, designated as R(0) in Equation 7.62