!***********************************************************************
!     3-D FDTD code with CPML absorbing boundary conditions
!***********************************************************************
!
!     Program author: Jamesina J. Simpson, Ph.D. Graduate Student
!                     National Science Foundation Fellow
!                     NU Computational Electromagnetics Laboratory
!                     Director: Allen Taflove
!                     Department of Electrical Engineering and 
!                                                      Computer Science
!                     Northwestern University
!                     2145 Sheridan Road
!                     Evanston, IL 60208
!                     j-simpson@northwestern.edu
!
!     Copyright 2005
!
!     This FORTRAN 90 code implements the finite-difference time-domain
!     solution of Maxwell's curl equations over a three-dimensional
!     Cartesian space lattice.  The grid is terminated by CPML absorbing
!     boundary conditions.  However, in a straightforward manner it can
!     be altered to have PEC planes at any of the outer grid boundaries.  
!     Also, the number of grid cells as well as the thickness of 
!     the PML in each Cartesian direction can be varied independently.
!
!     For illustrative purposes, the code supplied here models the PML 
!     numerical experiment described in Section 7.11.2.  Three output 
!     files are generated having the following data:
!
!       (1) Ez at the source for each time step;
!       (2) Ey at the probe point for each time step;  and
!       (3) The plane of Ez values 1 mm above the source in the k-
!           direction recorded at the time step, "record_grid".  This 
!           data can be viewed in Matlab using the following commands:
!
!             >>  load view_grid.dat
!             >>  for j = 1:Jmax
!             >>      image(1:Imax,j) = view_grid(1+(j-1)*Imax:j*Imax);
!             >>  end
!             >>  pcolor(image'); shading flat
!
!     The relative error (equation 7.135) graphed in Fig. 7.6 on 
!     page 318 can be reproduced by comparing the output file 
!     "probe.dat" with that of a much larger benchmark grid having 
!     Imax, Jmax, and Kmax increased to the values mentioned in the text.
!
!     This code has been tested using the Intel Fortran Compiler 8.0 for
!     Linux.  The executable can be created by typing
!     "ifort fdtd3D_CPML.f90" at the command prompt.
!
!     This program is not guaranteed to be free of defects or bugs.
!     Please report any bugs that may exist to:
!                                          j-simpson@northwestern.edu
!
!
!***********************************************************************


   PROGRAM fdtd3D_CPML

   IMPLICIT NONE

!  ..................................
!  Input Fundamental Constants (MKS units)
   REAL, PARAMETER ::                            &
      pi = 3.14159265358979, C = 2.99792458E8, &
      muO = 4.0 * pi * 1.0E-7, epsO = 1.0/(C*C*muO)

!  ..................................
!  Specify Material Relative Permittivity and Conductivity
   REAL, PARAMETER::                      &
      epsR = 1.0, sigM1 = 0.0   ! free space

!  ..................................
!  Specify Number of Time Steps and Grid Size Parameters
   INTEGER, PARAMETER ::                                     &
      nmax = 2100, &  ! total number of time steps
      Imax = 51, Jmax = 126, Kmax = 26  ! grid size corresponding to the
                                        ! number of Ez field components 
!  ..................................
!  Specify Grid Cell Size in Each Direction and Calculate the 
!  Resulting Courant-Stable Time Step
   REAL, PARAMETER ::                                        &
      dx = 1.0E-3, dy = 1.0E-3, dz = 1.0E-3,  &  ! cell size in each direction
      dt = 0.99 / (C*(1.0/dx**2+1.0/dy**2+1.0/dz**2)**0.5)
                                                 ! time step increment

!  ..................................
!  Specify the Impulsive Source (See Equation 7.134)
   REAL, PARAMETER ::                                        &
      tw = 53.0E-12, tO = 4.0*tw  

!  ..................................
!  Specify the Time Step at which the Grid is Recorded for Visualization
   INTEGER, PARAMETER ::                                        &
      record_grid = 300 

!  ..................................
!  Specify the PEC Plate Boundaries and the Source/Recording Points
   INTEGER, PARAMETER ::                                    &
      istart = (Imax-1)/2-11, iend = istart+24, jstart = istart,   &
      jend = jstart + 99, kstart = Kmax/2, kend = kstart,      &
      isource = istart, jsource = jstart, ksource = kend+1,     &
      irecv1 = iend, jrecv1 = jend+1, krecv1 = kend  ! Ey at probe point

!  ..................................
!  Specify the CPML Thickness in Each Direction (Value of Zero 
!  Corresponds to No PML, and the Grid is Terminated with a PEC)
   INTEGER, PARAMETER ::                        &
      ! PML thickness in each direction 
      nxPML_1 = 11, nxPML_2 = nxPML_1, nyPML_1 = nxPML_1,      &
      nyPML_2 = nxPML_1, nzPML_1 = nxPML_1, nzPML_2 = nxPML_2
!  ..................................
!  Specify the CPML Order and Other Parameters
   INTEGER, PARAMETER ::                        & 
      m = 3, ma = 1 
   REAL, PARAMETER  ::                     &
      sig_x_max = 0.75 * (0.8*(m+1)/(dx*(muO/epsO*epsR)**0.5)),   &
      sig_y_max = 0.75 * (0.8*(m+1)/(dy*(muO/epsO*epsR)**0.5)),   &
      sig_z_max = 0.75 * (0.8*(m+1)/(dz*(muO/epsO*epsR)**0.5)),   &
      alpha_x_max = 0.24,   &
      alpha_y_max = alpha_x_max, alpha_z_max = alpha_x_max, &
      kappa_x_max = 15.0, &
      kappa_y_max = kappa_x_max, kappa_z_max = kappa_x_max

   INTEGER ::                                                &
	i,j,ii,jj,k,kk,n

   REAL  ::                                                   &
      source, P1, P2

!     TM components
   REAL,DIMENSION(Imax, Jmax, Kmax)  ::                      &
      Ez, CA, CB, sig, eps

   REAL,DIMENSION(Imax-1, Jmax, Kmax)  ::                      &
      Hy

   REAL,DIMENSION(Imax,Jmax-1, Kmax)  ::                      &
      Hx

   REAL  ::                        &
      DA, DB

!     TE components
   REAL,DIMENSION(Imax-1, Jmax-1, Kmax-1)  ::                      &
      Hz

   REAL,DIMENSION(Imax-1, Jmax, Kmax-1)  ::                      &
      Ex

   REAL,DIMENSION(Imax,Jmax-1, Kmax-1)  ::                      &
      Ey

!  PML
   REAL ,DIMENSION(nxPML_1,Jmax,Kmax) ::                       &
      psi_Ezx_1

   REAL ,DIMENSION(nxPML_2,Jmax,Kmax) ::                       &
      psi_Ezx_2

   REAL ,DIMENSION(nxPML_1-1,Jmax,Kmax) ::                       &
      psi_Hyx_1

   REAL ,DIMENSION(nxPML_2-1,Jmax,Kmax) ::                       &
      psi_Hyx_2

   REAL ,DIMENSION(Imax,nyPML_1,Kmax) ::                       &
      psi_Ezy_1                               

   REAL ,DIMENSION(Imax,nyPML_2,Kmax) ::                       &
      psi_Ezy_2

   REAL ,DIMENSION(Imax,nyPML_1-1,Kmax) ::                       &
      psi_Hxy_1                               

   REAL ,DIMENSION(Imax,nyPML_2-1,Kmax) ::                       &
      psi_Hxy_2

   REAL ,DIMENSION(Imax,Jmax-1,nzPML_1-1) ::                       &
      psi_Hxz_1

   REAL ,DIMENSION(Imax,Jmax-1,nzPML_2-1) ::                       &
      psi_Hxz_2

   REAL ,DIMENSION(Imax-1,Jmax,nzPML_1-1) ::                       &
      psi_Hyz_1

   REAL ,DIMENSION(Imax-1,Jmax,nzPML_2-1) ::                       &
      psi_Hyz_2

   REAL ,DIMENSION(Imax-1,Jmax,nzPML_1) ::                       &
      psi_Exz_1

   REAL ,DIMENSION(Imax-1,Jmax,nzPML_2) ::                       &
      psi_Exz_2

   REAL ,DIMENSION(Imax,Jmax-1,nzPML_1) ::                       &
      psi_Eyz_1

   REAL ,DIMENSION(Imax,Jmax-1,nzPML_2) ::                       &
      psi_Eyz_2

   REAL ,DIMENSION(nxPML_1-1,Jmax-1,Kmax-1) ::                       &
      psi_Hzx_1
   REAL ,DIMENSION(nxPML_1,Jmax-1,Kmax-1) ::                       &
      psi_Eyx_1

   REAL ,DIMENSION(nxPML_2-1,Jmax-1,Kmax-1) ::                       &
      psi_Hzx_2
   REAL ,DIMENSION(nxPML_2,Jmax-1,Kmax-1) ::                       &
      psi_Eyx_2

   REAL ,DIMENSION(Imax-1,nyPML_1-1,Kmax-1) ::                       &
      psi_Hzy_1                               
   REAL ,DIMENSION(Imax-1,nyPML_1,Kmax-1) ::                       &
      psi_Exy_1                               

   REAL ,DIMENSION(Imax-1,nyPML_2-1,Kmax-1) ::                       &
      psi_Hzy_2
   REAL ,DIMENSION(Imax-1,nyPML_2,Kmax-1) ::                       &
      psi_Exy_2

   REAL ,DIMENSION(nxPML_1) ::                       &
      be_x_1, ce_x_1, alphae_x_PML_1, sige_x_PML_1, kappae_x_PML_1
   REAL ,DIMENSION(nxPML_1-1) ::                       &
      bh_x_1, ch_x_1, alphah_x_PML_1, sigh_x_PML_1, kappah_x_PML_1

   REAL ,DIMENSION(nxPML_2) ::                       &
      be_x_2, ce_x_2, alphae_x_PML_2, sige_x_PML_2, kappae_x_PML_2
   REAL ,DIMENSION(nxPML_2-1) ::                       &
      bh_x_2, ch_x_2, alphah_x_PML_2, sigh_x_PML_2, kappah_x_PML_2

   REAL ,DIMENSION(nyPML_1) ::                       &
      be_y_1, ce_y_1, alphae_y_PML_1, sige_y_PML_1, kappae_y_PML_1
   REAL ,DIMENSION(nyPML_1-1) ::                       &
      bh_y_1, ch_y_1, alphah_y_PML_1, sigh_y_PML_1, kappah_y_PML_1

   REAL ,DIMENSION(nyPML_2) ::                       &
      be_y_2, ce_y_2, alphae_y_PML_2, sige_y_PML_2, kappae_y_PML_2
   REAL ,DIMENSION(nyPML_2-1) ::                       &
      bh_y_2, ch_y_2, alphah_y_PML_2, sigh_y_PML_2, kappah_y_PML_2

   REAL ,DIMENSION(nzPML_1) ::                       &
      be_z_1, ce_z_1, alphae_z_PML_1, sige_z_PML_1, kappae_z_PML_1
   REAL ,DIMENSION(nzPML_1-1) ::                       &
      bh_z_1, ch_z_1, alphah_z_PML_1, sigh_z_PML_1, kappah_z_PML_1

   REAL ,DIMENSION(nzPML_2) ::                       &
      be_z_2, ce_z_2, alphae_z_PML_2, sige_z_PML_2, kappae_z_PML_2
   REAL ,DIMENSION(nzPML_2-1) ::                       &
      bh_z_2, ch_z_2, alphah_z_PML_2, sigh_z_PML_2, kappah_z_PML_2

!     denominators for the update equations
   REAL,DIMENSION(Imax-1)  ::                      &
      den_ex, den_hx

   REAL,DIMENSION(Jmax-1)  ::                      &
      den_ey, den_hy

   REAL,DIMENSION(Kmax-1)  ::                      &
      den_ez, den_hz


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  OPEN OUTPUT FILES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   OPEN (UNIT = 30, FILE = "source.dat")
   OPEN (UNIT = 31, FILE = "probe.dat")
   OPEN (UNIT = 33, FILE = "view_grid.dat")

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  INITIALIZE VARIABLES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   P1 = 0.0
   P2 = 0.0

   Ez(:,:,:) = 0.0
   Hz(:,:,:) = 0.0
   Ex(:,:,:) = 0.0
   Hx(:,:,:) = 0.0
   Ey(:,:,:) = 0.0
   Hy(:,:,:) = 0.0
   sig(:,:,:) = sigM1
   eps(:,:,:) = epsR*epsO

   psi_Exy_1(:,:,:) = 0.0
   psi_Exy_2(:,:,:) = 0.0
   psi_Exz_1(:,:,:) = 0.0
   psi_Exz_2(:,:,:) = 0.0
   psi_Eyx_1(:,:,:) = 0.0
   psi_Eyx_2(:,:,:) = 0.0
   psi_Eyz_1(:,:,:) = 0.0
   psi_Eyz_2(:,:,:) = 0.0
   psi_Ezy_1(:,:,:) = 0.0
   psi_Ezy_2(:,:,:) = 0.0
   psi_Ezx_1(:,:,:) = 0.0
   psi_Ezx_2(:,:,:) = 0.0
   psi_Hxy_1(:,:,:) = 0.0
   psi_Hxy_2(:,:,:) = 0.0
   psi_Hxz_1(:,:,:) = 0.0
   psi_Hxz_2(:,:,:) = 0.0
   psi_Hyx_1(:,:,:) = 0.0
   psi_Hyx_2(:,:,:) = 0.0
   psi_Hyz_1(:,:,:) = 0.0
   psi_Hyz_2(:,:,:) = 0.0
   psi_Hzy_1(:,:,:) = 0.0
   psi_Hzy_2(:,:,:) = 0.0
   psi_Hzx_1(:,:,:) = 0.0
   psi_Hzx_2(:,:,:) = 0.0

   write(*,*)"Imax: ", Imax
   write(*,*)"Jmax: ", Jmax
   write(*,*)"Kmax: ", Kmax
   write(*,*)"dt: ", dt
   write(*,*)"nmax: ", nmax
   write(*,*)"max time: ", nmax*dt
   write(*,*)"record grid after ", record_grid, "dt"

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  SET CPML PARAMETERS IN EACH DIRECTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO i = 1,nxPML_1
      sige_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i) / (nxPML_1 - 1.0) )**m
      alphae_x_PML_1(i) = alpha_x_max*((i-1.0)/(nxPML_1-1.0))**ma
      kappae_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*   &