# Export a composite data type PML

export PML

# Export a function mesti_build_fdfd_matrix()

export mesti_build_fdfd_matrix

mutable struct PML

npixels::Integer

sigma_max_over_omega::Real

power_sigma::Real

alpha_max_over_omega::Real

power_alpha::Real

kappa_max::Real

power_kappa::Real

# Below are only used in mesti() and mesti2s()

direction::String

side::String

# Below is only used in mesti2s()

npixels_spacer::Union{Integer,Nothing}

# Construct default parameters from Table 7.1 of Taflove & Hagness's 2005 FDTD book

# Without specifying pml.npixels

PML() = (pml = new(); return pml)

# With specifying pml.npixels

PML(n) = (pml = new(); pml.npixels=n; return pml)

end

"""

MESTI_BUILD_FDFD_MATRIX The finite-difference frequency-domain operator in 3D.

A = mesti_build_fdfd_matrix(epsilon_xx, epsilon_xy, epsilon_xz, epsilon_yx, epsilon_yy, epsilon_yz, epsilon_zx, epsilon_zy, epsilon_zz, k0dx, xBC, yBC, zBC, xPML, yPML, zPML, use_UPML) returns A as a sparse matrix representing

wave operator [curl(curl) - k0^2*epsilon_r(x,y,z)]*(dx^2)

discretized on a square grid with grid size dx through center difference.

for the electric field: (Ex, Ey, Ez).

=== Input Arguments ===

epsilon_xx (array; required):

epsilon_xx is a nx_Ex-by-ny_Ex-by-nz_Ex numeric array (real or complex)

discretizing the relative permittivity profile epsilon_xx on Ex-sites.

The notation ni_Ej is the total number of grids for Ej in i-direction.

epsilon_xy (array; optional):

epsilon_xy is a nx_Ez-by-ny_Ez-by-nz_Ex numeric array (real or complex)

discretizing the relative permittivity profile epsilon_xy on lower corner of the Yee lattice.

epsilon_xz (array; optional):

epsilon_xz is a nx_Ey-by-ny_Ex-by-nz_Ey numeric array (real or complex)

discretizing the relative permittivity profile epsilon_xz on lower corner of the Yee lattice.

epsilon_yx (array; optional):

epsilon_yx is a nx_Ez-by-ny_Ez-by-nz_Ey numeric array (real or complex)

discretizing the relative permittivity profile epsilon_yx on lower corner of the Yee lattice.

epsilon_yy (array; required):

epsilon_yy is a nx_Ey-by-ny_Ey-by-nz_Ey numeric array (real or complex)

discretizing the relative permittivity profile epsilon_yy on Ey-sites.

epsilon_yz (array; optional):

epsilon_yz is a nx_Ey-by-ny_Ex-by-nz_Ex numeric array (real or complex)

discretizing the relative permittivity profile epsilon_yz on lower corner of the Yee lattice.

epsilon_zx (array; optional):

epsilon_zx is a nx_Ey-by-ny_Ez-by-nz_Ey numeric array (real or complex)

discretizing the relative permittivity profile epsilon_zx on lower corner of the Yee lattice.

epsilon_zy (array; optional):

epsilon_zy is a nx_Ez-by-ny_Ex-by-nz_Ex numeric array (real or complex)

discretizing the relative permittivity profile epsilon_zy on lower corner of the Yee lattice.

epsilon_zz (array; required):

epsilon_zz is a nx_Ez-by-ny_Ez-by-nz_Ez numeric array (real or complex)

discretizing the relative permittivity profile epsilon_zz on Ez-sites.

k0dx (numeric scalar, real or complex; required):

Normalized frequency k0*dx = (2*pi/vacuum_wavelength)*dx.

xBC (character vector or numeric scalar; required):

Boundary condition (BC) at the two ends in x direction, effectively

specifying Ex(n,m,l) at n=0 and n=nx_Ex 1,

Ey(n,m,l) at n=0 and n=nx_Ey 1,

Ez(n,m,l) at n=0 and n=nx_Ez 1.

one pixel beyond the computation domain. Available choices are:

"periodic" - Ex(n nx_Ex,m,l) = Ex(n,m,l),

Ey(n nx_Ey,m,l) = Ey(n,m,l),

Ez(n nx_Ez,m,l) = Ez(n,m,l).

"PEC" - Ex(0,m,l) = Ex(1,m,l); Ex(nx_Ex 1,m,l) = Ex(nx_Ex,m,l),

Ey(0,m,l) = Ey(nx_Ey 1,m,l) = 0,

Ez(0,m,l) = Ez(nx_Ez 1,m,l) = 0.

"PMC" - Ex(0,m,l) = Ex(nx_Ex 1,m,l) = 0,

Ey(0,m,l) = Ey(1,m,l); Ey(nx_Ey 1,m,l) = Ey(nx_Ey,m,l),

Ez(0,m,l) = Ez(1,m,l); Ez(nx_Ez 1,m,l) = Ez(nx_Ez,m,l).

"PECPMC" - Ex(0,m,l) = Ex(1,m,l); Ex(nx_Ex 1,m,l) = 0,

Ey(0,m,l) = 0; Ey(nx_Ey 1,m,l) = Ey(nx_Ey,m,l),

Ez(0,m,l) = 0; Ez(nx_Ez 1,m,l) = Ez(nx_Ez,m,l),

"PMCPEC" - Ex(0,m,l) = 0; Ex(nx_Ex 1,m,l) = Ex(nx_Ex,m,l),

Ey(0,m,l) = Ey(1,m,l); Ey(nx_Ey 1,m,l) = 0,

Ez(0,m,l) = Ez(1,m,l); Ez(nx_Ez 1,m,l) = 0.

where PEC stands for perfect electric conductor and PMC stands for perfect

magnetic conductor.

When xBC is a numeric scalar, the Bloch periodic boundary condition is

used with p(n nx,m,l) = p(n,m,l)*exp(1i*xBC), p = Ex,Ey,Ez; in other words,

xBC = kx_B*nx*dx = kx_B*Lambda where kx_B is the Bloch wave number and

Lambda = nx*dx is the periodicity in x. In this case, nx = nx_Ex = nx_Ey = nx_Ez.

yBC (character vector or numeric scalar; required):

Boundary condition in y direction, analogous to xBC.

zBC (character vector or numeric scalar; required):

Boundary condition in z direction, analogous to xBC.

xPML (a two-element PML vector; optional):

Parameters for perfectly matched layer (PML) in x direction.

xPML = [PML_left, PML_right]

If users do not want to put PML on either side, just set the field "npixels" = 0.

For example, PML_left.npixels = 0: no PML on the left.

In each case, PML is a PML structure with the following fields:

npixels (non-negative integer scalar; required): Number of PML pixels.

Note this is within syst.epsilon or syst.inv_epsilon.

power_sigma (non-negative scalar; optional):

Power of the polynomial grading for the conductivity sigma;

defaults to 3.

sigma_max_over_omega (non-negative scalar; optional):

Conductivity at the end of the PML; By default, it is

set to an optimized value based on resolution and the

background refractive index. This is used to attenuate propagating waves.

power_kappa (non-negative scalar; optional):

Power of the polynomial grading for the real-coordinate-stretching

factor kappa; defaults to 3.

kappa_max (real scalar no smaller than 1; optional):

Real-coordinate-stretching factor at the end of the PML; By default,

it is set to an optimized value based on resolution and the

background refractive index. This is used to accelerate the attenuation

of evanescent waves. kappa_max = 1 means no real-coordinate stretching.

power_alpha (non-negative scalar; optional):

Power of the polynomial grading for the CFS alpha factor;

defaults to 1.

alpha_max_over_omega (non-negative scalar; optional):

Complex-frequency-shifting (CFS) factor at the beginning

of the PML. This is typically used in time-domain simulations

to suppress late-time (low-frequency) reflections.

We don't use it by default (alpha_max_over_omega = 0)

since we are in frequency domain.

We use the following PML coordinate-stretching factor:

s(p) = kappa(p) sigma(p)./(alpha(p) - i*omega)

with

sigma(p)/omega = sigma_max_over_omega*(p.^power_sigma),

kappa(p) = 1 (kappa_max-1)*(p.^power_kappa),

alpha(p)/omega = alpha_max_over_omega*((1-p).^power_alpha),

where omega is frequency, and p goes linearly from 0 at the beginning of

the PML to 1 at the end of the PML.

yPML (a two-element vector; optional):

Parameters for PML in y direction, analogous to xPML.

zPML (a two-element vector; optional):

Parameters for PML in z direction, analogous to xPML.

use_UPML (logical scalar; optional, defaults to true):

Whether to use uniaxial PML (UPML) or not. If not, stretched-coordinate

PML (SC-PML) will be used.

=== Output Arguments ===

A (sparse matrix):

(nt_Ex nt_Ey nt_Ez)-by-(nt_Ex nt_Ey nt_Ez) sparse matrix representing

the 3D FDFD operator. nt_Ex/nt_Ey/nt_Ez is the total number of grids

for Ex/Ey/Ez.

is_symmetric_A (logical scalar):

Whether matrix A is symmetric or not.

xPML (two-element vector):

PML parameters used on the low and high sides of x direction, if any.

yPML (two-element vector):

PML parameters used on the low and high sides of y direction, if any.

zPML (two-element vector):

PML parameters used on the low and high sides of z direction, if any.

"""

function mesti_build_fdfd_matrix(epsilon_xx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Matrix{Int64},Matrix{Float64},Matrix{ComplexF64}}, epsilon_xy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_xz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_yx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_yy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_yz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_zx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_zy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_zz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, k0dx::Union{Real,Complex}, xBC::Union{String,Real,Complex,Nothing}, yBC::Union{String,Real,Complex}, zBC::Union{String,Real,Complex}, xPML::Union{Vector{PML},Nothing} = [PML(0), PML(0)], yPML::Vector{PML} = [PML(0), PML(0)], zPML::Vector{PML} = [PML(0), PML(0)], use_UPML::Bool=true)

# Make deepcopy of them to avoid mutating input argument

xPML = deepcopy(xPML); yPML = deepcopy(yPML); zPML = deepcopy(zPML)

# Take care of the 2D TM case

if ndims(epsilon_xx) == 2

use_2D_TM = true

if ~isa(epsilon_yy, Nothing) || ~isa(epsilon_zz, Nothing) || ~isa(epsilon_xy, Nothing) || ~isa(epsilon_xz, Nothing) || ~isa(epsilon_yx, Nothing) || ~isa(epsilon_yz, Nothing) || ~isa(epsilon_zx, Nothing) || ~isa(epsilon_zy, Nothing)

throw(ArgumentError("Only epsilon_xx is required for 2D TM fields Ex(y,z), but other components should not be given or should be nothing"))

end

if xBC != nothing

@warn "Only yBC and zBC are required for 2D TM fields Ex(y,z). xBC will be ignored."

end

if xPML != nothing

@warn "Only yPML and zPML are required for 2D TM fields Ex(y,z). xPML will be ignored."

end

else

use_2D_TM = false

end

# Check the presence off-diagonal part in epsilon

if ~isa(epsilon_xy, Nothing) || ~isa(epsilon_xz, Nothing) || ~isa(epsilon_yx, Nothing) || ~isa(epsilon_yz, Nothing) || ~isa(epsilon_zx, Nothing) || ~isa(epsilon_zy, Nothing)

include_off_diagonal_epsilon = true

else

include_off_diagonal_epsilon = false

end

if use_2D_TM # 2D TM field

# Number of sites in y and z for Ex

(ny_Ex, nz_Ex) = size(epsilon_xx)

# Convert BC to take care of lowercase or uppercase

yBC = convert_BC(yBC, "y")

zBC = convert_BC(zBC, "z")

# Total number of grid points for Ex

nt_Ex = ny_Ex*nz_Ex

# Checking the whether the permittivity is lossless or not

if (maximum(imag(epsilon_xx)) >= sqrt(eps()))

@warn("syst.epsilon_xx contains imaginary part; the permittivity profile is not lossless.")

end

else # 3D case

# Number of grid points in x, y, and z for Ex, Ey, and Ez

(nx_Ex, ny_Ex, nz_Ex) = size(epsilon_xx)

(nx_Ey, ny_Ey, nz_Ey) = size(epsilon_yy)

(nx_Ez, ny_Ez, nz_Ez) = size(epsilon_zz)

# Convert BC to take care of lowercase or uppercase

xBC = convert_BC(xBC, "x")

yBC = convert_BC(yBC, "y")

zBC = convert_BC(zBC, "z")

# Check number of grid points with the boundary conditions

if nx_Ey != nx_Ez; throw(ArgumentError("Number of grids along x provided by epsilon_yy and epsilon_zz should be same.")); end

if ny_Ex != ny_Ez; throw(ArgumentError("Number of grids along y provided by epsilon_xx and epsilon_zz should be same.")); end

if nz_Ex != nz_Ey; throw(ArgumentError("Number of grids along z provided by epsilon_xx and epsilon_yy should be same.")); end

check_BC_and_grid(xBC, nx_Ex, nx_Ey, nx_Ez, "x")

check_BC_and_grid(yBC, ny_Ex, ny_Ey, ny_Ez, "y")

check_BC_and_grid(zBC, nz_Ex, nz_Ey, nz_Ez, "z")

if (~isa(epsilon_xy, Nothing) && ~(size(epsilon_xy) == (nx_Ez, ny_Ez, nz_Ex)))

throw(ArgumentError("The size of epsilon_xy should be should be (size(epsilon_zz, 1), size(epsilon_zz, 2), size(epsilon_xx, 3)) = ($(size(epsilon_zz, 1)), $(size(epsilon_zz, 2)), $(size(epsilon_xx, 3)))."))

end

if (~isa(epsilon_xz, Nothing) && ~(size(epsilon_xz) == (nx_Ey, ny_Ex, nz_Ey)))

throw(ArgumentError("The size of epsilon_xz should be should be (size(epsilon_yy, 1), size(epsilon_xx, 2), size(epsilon_yy, 3)) = ($(size(epsilon_yy, 1)), $(size(epsilon_xx, 2)), $(size(epsilon_yy, 3)))."))

end

if (~isa(epsilon_yx, Nothing) && ~(size(epsilon_yx) == (nx_Ez, ny_Ez, nz_Ey)))

throw(ArgumentError("The size of epsilon_yx should be should be (size(epsilon_zz, 1), size(epsilon_zz, 2), size(epsilon_yy, 3)) = ($(size(epsilon_zz, 1)), $(size(epsilon_zz, 2)), $(size(epsilon_yy, 3)))."))

end

if (~isa(epsilon_yz, Nothing) && ~(size(epsilon_yz) == (nx_Ey, ny_Ex, nz_Ex)))

throw(ArgumentError("The size of epsilon_yz should be should be (size(epsilon_yy, 1), size(epsilon_xx, 2), size(epsilon_xx, 3)) = ($(size(epsilon_yy, 1)), $(size(epsilon_xx, 2)), $(size(epsilon_xx, 3)))."))

end

if (~isa(epsilon_zx, Nothing) && ~(size(epsilon_zx) == (nx_Ey, ny_Ez, nz_Ey)))

throw(ArgumentError("The size of epsilon_zx should be should be (size(epsilon_yy, 1), size(epsilon_zz, 2), size(epsilon_yy, 3)) = ($(size(epsilon_yy, 1)), $(size(epsilon_zz, 2)), $(size(epsilon_yy, 3)))."))

end

if (~isa(epsilon_zy, Nothing) && ~(size(epsilon_zy) == (nx_Ez, ny_Ex, nz_Ex)))

throw(ArgumentError("The size of epsilon_zy should be should be (size(epsilon_zz, 1), size(epsilon_xx, 2), size(epsilon_xx, 3)) = ($(size(epsilon_zz, 1)), $(size(epsilon_xx, 2)), $(size(epsilon_xx, 3)))."))

end

# Checking the whether the permittivity is lossless or not

if (maximum(imag(epsilon_xx)) >= sqrt(eps()))

@warn("syst.epsilon_xx contains imaginary part; the permittivity profile is not lossless.")

end

if (maximum(imag(epsilon_yy)) >= sqrt(eps()))

@warn("syst.epsilon_yy contains imaginary part; the permittivity profile is not lossless.")

end

if (maximum(imag(epsilon_zz)) >= sqrt(eps()))

@warn("syst.epsilon_zz contains imaginary part; the permittivity profile is not lossless.")

end

if (~isa(epsilon_xy, Nothing) && ~isa(epsilon_yx, Nothing) && maximum(abs.(epsilon_xy - conj(epsilon_yx))) >= sqrt(eps()))

@warn("syst.epsilon_xy != conj(syst.epsilon_yx); the permittivity profile is not lossless.")

end

if (~isa(epsilon_xz, Nothing) && ~isa(epsilon_zx, Nothing) && maximum(abs.(epsilon_xz - conj(epsilon_zx))) >= sqrt(eps()))

@warn("syst.epsilon_xz != conj(syst.epsilon_zx); the permittivity profile is not lossless.")

end

if (~isa(epsilon_yz, Nothing) && ~isa(epsilon_zy, Nothing) && maximum(abs.(epsilon_yz - conj(epsilon_zy))) >= sqrt(eps()))

@warn("syst.epsilon_yz != conj(syst.epsilon_zy); the permittivity profile is not lossless.")

end

# Total number of grid points for Ex, Ey, and Ez

nt_Ex = nx_Ex*ny_Ex*nz_Ex

nt_Ey = nx_Ey*ny_Ey*nz_Ey

nt_Ez = nx_Ez*ny_Ez*nz_Ez

end

# Estimate background permittivity, used to assign default sigma_max_over_omega for PML

if use_2D_TM # 2D TM field

epsilon_bg_y = [1, 1]

epsilon_bg_z = [1, 1]

if yPML[1].npixels != 0 || yPML[2].npixels != 0

epsilon_bg_y = real(vcat(mean(epsilon_xx[1,:]), mean(epsilon_xx[end,:])))

# Make sure that the two sides are the same when the system is periodic; this ensures that the s-factor is continuous across the periodic boundary.

if isa(yBC, Number) || yBC == "periodic"

epsilon_bg_y = mean(epsilon_bg_y)*[1, 1]

end

end

if zPML[1].npixels != 0 || zPML[2].npixels != 0

epsilon_bg_z = real(vcat(mean(epsilon_xx[:,1]), mean(epsilon_xx[:,end])))

if isa(zBC, Number) || zBC == "periodic"

epsilon_bg_z = mean(epsilon_bg_z)*[1, 1]

end

end

# Set default values for PML parameters

yPML = mesti_set_PML_params(yPML, k0dx, epsilon_bg_y, "y")

zPML = mesti_set_PML_params(zPML, k0dx, epsilon_bg_z, "z")

else # 3D case

epsilon_bg_x_Ex = [1, 1]

epsilon_bg_y_Ey = [1, 1]

epsilon_bg_z_Ez = [1, 1]

if xPML[1].npixels != 0 || xPML[2].npixels != 0

epsilon_bg_x_Ex = real(vcat(mean(epsilon_xx[1,:,:]), mean(epsilon_xx[end,:,:])))

# Make sure that the two sides are the same when the system is periodic; this ensures that the s-factor is continuous across the periodic boundary.

if isa(xBC, Number) || xBC == "periodic"

epsilon_bg_x_Ex = mean(epsilon_bg_x_Ex)*[1, 1]

end

end

if yPML[1].npixels != 0 || yPML[2].npixels != 0

epsilon_bg_y_Ey = real(vcat(mean(epsilon_yy[:,1,:]), mean(epsilon_yy[:,end,:])))

if isa(yBC, Number) || yBC == "periodic"

epsilon_bg_y_Ey = mean(epsilon_bg_y_Ey)*[1, 1]

end

end

if zPML[1].npixels != 0 || zPML[2].npixels != 0

epsilon_bg_z_Ez = real(vcat(mean(epsilon_zz[:,:,1]), mean(epsilon_zz[:,:,end])))

if isa(zBC, Number) || zBC == "periodic"

epsilon_bg_z_Ez = mean(epsilon_bg_z_Ez)*[1, 1]

end

end

# Set default values for PML parameters

xPML = mesti_set_PML_params(xPML, k0dx, epsilon_bg_x_Ex, "x")

yPML = mesti_set_PML_params(yPML, k0dx, epsilon_bg_y_Ey, "y")

zPML = mesti_set_PML_params(zPML, k0dx, epsilon_bg_z_Ez, "z")

end

# Build the first derivative and use Kronecker outer product to go from 1D to 2D or from 1D to 3D

if use_2D_TM

# Build the first derivative and average matrices on E, such that

# ddx_E*p = dp/dx

# avg_x_E*p = average of p among two neighboring pixels

# where p = Ex is a 1d vector.

# Note that sx_E and sx_H are column vectors.

# Later we will use Kronecker outer product to go from 1D to 2D.

(ddy_E, avg_y_E, sy_E, sy_H, ind_yPML_E) = build_ddx_E(ny_Ex, yBC, yPML, "y") # ddy_E operates on Ex or Ez

(ddz_E, avg_z_E, sz_E, sz_H, ind_zPML_E) = build_ddx_E(nz_Ex, zBC, zPML, "z") # ddz_E operates on Ex or Ey

# The derivative matrices on H

# We need sparse() to force datatype conversion.

ddy_H = sparse(-ddy_E') # ddy_H operates on Hx or Hz

ddz_H = sparse(-ddz_E') # ddz_H operates on Hx or Hy

# Number of grid points in y and z for Hx

ny_Hx = size(ddy_E, 1) # ny_Hx = ny_Hz = ny_Ey

nz_Hx = size(ddz_E, 1) # nz_Hx = nz_Hy = nz_Ez

# Include coordinate streching into the first derivatives

# So far we only implement for 2D TM fields Ex(y,z); will add 2D TE fields later

ddy_E = spdiagm(ny_Hx, ny_Hx, 1 ./sy_H)*ddy_E

ddz_E = spdiagm(nz_Hx, nz_Hx, 1 ./sz_H)*ddz_E

if ~use_UPML

ddy_H = spdiagm(ny_Ex, ny_Ex, 1 ./sy_E)*ddy_H

ddz_H = spdiagm(nz_Ex, nz_Ex, 1 ./sz_E)*ddz_H

end

# Build the operators; use Kronecker outer product to go from 1D to 2D.

# TM: A = [- (d/dy)^2 - (d/dz)^2 - k0^2*epsilon(y,z)]*(dx^2)

if ~use_UPML

# Stretched-coordinate PML (SC-PML); here, both 1/s_E and 1/s_H have already been multiplied onto the first derivatives.

A = - kron(ddz_H*ddz_E, sparse(I,ny_Ex,ny_Ex)) - kron(sparse(I,nz_Ex,nz_Ex), ddy_H*ddy_E) - spdiagm(nt_Ex, nt_Ex, (k0dx^2)*epsilon_xx[:])

else

# Uniaxial PML (UPML); here, the first 1/s has already been multiplied onto the first derivatives.

# A_UPML = (sy_E*sz_E)*A_SCPML for TM fields

syz_E = sy_E.*reshape(sz_E, 1, :)

A = - kron(ddz_H*ddz_E, spdiagm(ny_Ex, ny_Ex, sy_E)) - kron(spdiagm(nz_Ex, nz_Ex, sz_E), ddy_H*ddy_E) - spdiagm(nt_Ex, nt_Ex, (k0dx^2)*(syz_E[:].*epsilon_xx[:]))

end

else

# Build the first derivatives in 1D, such that ddx*p = dp/dx, ddy*p = dp/dy, ddz*p = dp/dz

# Later we will use Kronecker outer product to go from 1D to 3D.

# The operator notations follow

# ddx_HzEy: x-derivative operates on Ey producing Hz

# avg_x_Ey: average of Ey among two neighboring pixels along x-direction

# sx_Ey: coordinate-stretching factor for Ey along x-direction

# sx_Hz: coordinate-stretching factor for Hz along x-direction

# ddx_HyEz = ddx_HzEy, sx_Ez = sx_Ey, sx_Hy = sx_Hz, because Ez(Hy) and Ey(Hz)

# share the same x-coordinate inside same Yee-cell. The same reason applies to others.

# We defer three of the derivatives to later

(ddx_HzEy, avg_x_Ey, sx_Ey, sx_Hz, ind_xPML_E) = build_ddx_E(nx_Ey, xBC, xPML, "x")

#ddx_HyEz = ddx_HzEy

(ddy_HxEz, avg_y_Ez, sy_Ez, sy_Hx, ind_yPML_E) = build_ddx_E(ny_Ez, yBC, yPML, "y")

#ddy_HzEx = ddy_HxEz

(ddz_HyEx, avg_z_Ex, sz_Ex, sz_Hy, ind_zPML_E) = build_ddx_E(nz_Ex, zBC, zPML, "z")

#ddz_HxEy = ddz_HyEx

avg_x_Ex = build_ave_x_Ex(nx_Ex, xBC, "x")

avg_y_Ey = build_ave_x_Ex(ny_Ey, yBC, "y")

avg_z_Ez = build_ave_x_Ex(nz_Ez, zBC, "z")

avg_x_Ez = avg_x_Ey

avg_y_Ex = avg_y_Ez

avg_z_Ey = avg_z_Ex

# The derivative matrices on H

# We need sparse() to force datatype conversion.

ddx_EyHz = sparse(-ddx_HzEy'); #ddx_EzHy = sparse(-ddx_HyEz')

ddy_EzHx = sparse(-ddy_HxEz'); #ddy_ExHz = sparse(-ddy_HzEx')

ddz_ExHy = sparse(-ddz_HyEx'); #ddz_EyHx = sparse(-ddz_HxEy')

# Number of grid points in x, y, and z for Hx, Hy, and Hz

nx_Hz = size(ddx_HzEy, 1); nx_Hy = nx_Hz; nx_Hx = nx_Ey

ny_Hx = size(ddy_HxEz, 1); ny_Hz = ny_Hx; ny_Hy = ny_Ez

nz_Hy = size(ddz_HyEx, 1); nz_Hx = nz_Hy; nz_Hz = nz_Ex

# Total number of grid points for Hx, Hy, and Hz

nt_Hx = nx_Hx*ny_Hx*nz_Hx

nt_Hy = nx_Hy*ny_Hy*nz_Hy

nt_Hz = nx_Hz*ny_Hz*nz_Hz

# Start from stretched-coordinate PML (SC-PML); here, 1/s and 1/s_d are both multiplied onto the gradient.

# Convert coordinate-stretching factor vector to diagonal matrices

inv_sx_Hz = spdiagm(nx_Hz, nx_Hz, 1 ./sx_Hz); #inv_sx_Hy = inv_sx_Hz

inv_sy_Hx = spdiagm(ny_Hx, ny_Hx, 1 ./sy_Hx); #inv_sy_Hz = inv_sy_Hx

inv_sz_Hy = spdiagm(nz_Hy, nz_Hy, 1 ./sz_Hy); #inv_sz_Hx = inv_sz_Hy

inv_sx_Ey = spdiagm(nx_Ey, nx_Ey, 1 ./sx_Ey); #inv_sx_Ez = inv_sx_Ey

inv_sy_Ez = spdiagm(ny_Ez, ny_Ez, 1 ./sy_Ez); #inv_sy_Ex = inv_sy_Ez

inv_sz_Ex = spdiagm(nz_Ex, nz_Ex, 1 ./sz_Ex); #inv_sz_Ey = inv_sz_Ex

# Include coordinate stretching into the first derivatives

ddx_HzEy = inv_sx_Hz*ddx_HzEy; ddx_HyEz = ddx_HzEy #ddx_HyEz = inv_sx_Hy*ddx_HyEz

ddy_HxEz = inv_sy_Hx*ddy_HxEz; ddy_HzEx = ddy_HxEz #ddy_HzEx = inv_sy_Hz*ddy_HzEx

ddz_HyEx = inv_sz_Hy*ddz_HyEx; ddz_HxEy = ddz_HyEx #ddz_HxEy = inv_sz_Hx*ddz_HxEy

ddx_EyHz = inv_sx_Ey*ddx_EyHz; ddx_EzHy = ddx_EyHz #ddx_EzHy = inv_sx_Ez*ddx_EzHy

ddy_EzHx = inv_sy_Ez*ddy_EzHx; ddy_ExHz = ddy_EzHx #ddy_ExHz = inv_sy_Ex*ddy_ExHz

ddz_ExHy = inv_sz_Ex*ddz_ExHy; ddz_EyHx = ddz_ExHy #ddz_EyHx = inv_sz_Ey*ddz_EyHx

# Expand derivative matrices to 3D by Kronecker outer product

# Note that even though ddx_HzEy = ddx_HyEz, Dx_HzEy and Dx_HyEz are not the same

# because Ez and Ey can have different number of points along y and z

Dx_HzEy = kron(sparse(I,nz_Ey,nz_Ey), kron(sparse(I,ny_Ey,ny_Ey), ddx_HzEy))

Dx_HyEz = kron(sparse(I,nz_Ez,nz_Ez), kron(sparse(I,ny_Ez,ny_Ez), ddx_HyEz))

Dy_HxEz = kron(sparse(I,nz_Ez,nz_Ez), kron(ddy_HxEz, sparse(I,nx_Ez,nx_Ez)))

Dy_HzEx = kron(sparse(I,nz_Ex,nz_Ex), kron(ddy_HzEx, sparse(I,nx_Ex,nx_Ex)))

Dz_HyEx = kron(ddz_HyEx, kron(sparse(I,ny_Ex,ny_Ex), sparse(I,nx_Ex,nx_Ex)))

Dz_HxEy = kron(ddz_HxEy, kron(sparse(I,ny_Ey,ny_Ey), sparse(I,nx_Ey,nx_Ey)))

Dx_EyHz = kron(sparse(I,nz_Ey,nz_Ey), kron(sparse(I,ny_Ey,ny_Ey), ddx_EyHz))

Dx_EzHy = kron(sparse(I,nz_Ez,nz_Ez), kron(sparse(I,ny_Ez,ny_Ez), ddx_EzHy))

Dy_EzHx = kron(sparse(I,nz_Ez,nz_Ez), kron(ddy_EzHx, sparse(I,nx_Ez,nx_Ez)))

Dy_ExHz = kron(sparse(I,nz_Ex,nz_Ex), kron(ddy_ExHz, sparse(I,nx_Ex,nx_Ex)))

Dz_ExHy = kron(ddz_ExHy, kron(sparse(I,ny_Ex,ny_Ex), sparse(I,nx_Ex,nx_Ex)))

Dz_EyHx = kron(ddz_EyHx, kron(sparse(I,ny_Ey,ny_Ey), sparse(I,nx_Ey,nx_Ey)))

# Construct curl_E and curl_H

# curl_E: curl operator matrix operates on E-field producing H-field.

# curl_H: curl operator matrix operates on H-field producing E-field.

curl_E = vcat(hcat(spzeros(nt_Hx, nt_Ex), -Dz_HxEy, Dy_HxEz),

hcat( Dz_HyEx, spzeros(nt_Hy, nt_Ey), -Dx_HyEz),

hcat( -Dy_HzEx, Dx_HzEy, spzeros(nt_Hz, nt_Ez)))

curl_H = vcat(hcat(spzeros(nt_Ex, nt_Hx), -Dz_ExHy, Dy_ExHz),

hcat( Dz_EyHx, spzeros(nt_Ey, nt_Hy), -Dx_EyHz),

hcat( -Dy_EzHx, Dx_EzHy, spzeros(nt_Ez, nt_Hz)))

# Construct the vectorial Maxwell matrix A

epsilon_diagonal = spdiagm(nt_Ex nt_Ey nt_Ez, nt_Ex nt_Ey nt_Ez, vcat(epsilon_xx[:], epsilon_yy[:], epsilon_zz[:]))

A = curl_H*curl_E-(k0dx)^2*epsilon_diagonal

# Construct the off-diagonal part from the the relative permittivity tensor epsilon_ij, when i does not equal j

if include_off_diagonal_epsilon

A_off_diagonal_epsilon = spzeros(ComplexF64, nt_Ex nt_Ey nt_Ez, nt_Ex nt_Ey nt_Ez)

# Following Oskooi et al, Optics Letters 34, 2778 (2009), we average two points of Ey along y, multiply by epsilon_xy to get Ex, and then average two such points along direction x. To summarize: avg_x*epsilon_xy*avg_y*Ey. Similarly for the other terms.

if ~isa(epsilon_xy, Nothing)

matrix_epsilon_xy = spdiagm(nx_Ez*ny_Ez*nz_Ex, nx_Ez*ny_Ez*nz_Ex, epsilon_xy[:])

A_off_diagonal_epsilon[1:nt_Ex, (nt_Ex 1):(nt_Ex nt_Ey)] = A_off_diagonal_epsilon[1:nt_Ex, (nt_Ex 1):(nt_Ex nt_Ey)] - (k0dx)^2 * kron(sparse(I,nz_Ex,nz_Ex), kron(sparse(I,ny_Ex,ny_Ex), avg_x_Ex'))*matrix_epsilon_xy*kron(sparse(I,nz_Ey,nz_Ey), kron(avg_y_Ey, sparse(I,nx_Ey,nx_Ey))) # avg_x*epsilon_xy*avg_y*Ey

end

if ~isa(epsilon_xz, Nothing)

matrix_epsilon_xz = spdiagm(nx_Ey*ny_Ex*nz_Ey, nx_Ey*ny_Ex*nz_Ey, epsilon_xz[:])

A_off_diagonal_epsilon[1:nt_Ex, (nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez)] = A_off_diagonal_epsilon[1:nt_Ex, (nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez)] - (k0dx)^2 * kron(sparse(I,nz_Ex,nz_Ex), kron(sparse(I,ny_Ex,ny_Ex), avg_x_Ex'))*matrix_epsilon_xz*kron(avg_z_Ez, kron(sparse(I,ny_Ez,ny_Ez), sparse(I,nx_Ez,nx_Ez))) # avg_x*epsilon_xz*avg_z*Ez

end

if ~isa(epsilon_yx, Nothing)

matrix_epsilon_yx = spdiagm(nx_Ez*ny_Ez*nz_Ey, nx_Ez*ny_Ez*nz_Ey, epsilon_yx[:])

A_off_diagonal_epsilon[(nt_Ex 1):(nt_Ex nt_Ey), 1:nt_Ex] = A_off_diagonal_epsilon[(nt_Ex 1):(nt_Ex nt_Ey), 1:nt_Ex] - (k0dx)^2 * kron(sparse(I,nz_Ey,nz_Ey), kron(avg_y_Ey', sparse(I,nx_Ey,nx_Ey)))*matrix_epsilon_yx*kron(sparse(I,nz_Ex,nz_Ex), kron(sparse(I,ny_Ex,ny_Ex), avg_x_Ex)) # avg_y*epsilon_yx*avg_x*Ex

end

if ~isa(epsilon_yz, Nothing)

matrix_epsilon_yz = spdiagm(nx_Ey*ny_Ex*nz_Ex, nx_Ey*ny_Ex*nz_Ex, epsilon_yz[:])

A_off_diagonal_epsilon[(nt_Ex 1):(nt_Ex nt_Ey), (nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez)] = A_off_diagonal_epsilon[(nt_Ex 1):(nt_Ex nt_Ey), (nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez)] - (k0dx)^2 * kron(sparse(I,nz_Ey,nz_Ey), kron(avg_y_Ey', sparse(I,nx_Ey,nx_Ey)))*matrix_epsilon_yz*kron(avg_z_Ez, kron(sparse(I,ny_Ez,ny_Ez), sparse(I,nx_Ez,nx_Ez))) # avg_y*epsilon_yz*avg_z*Ez

end

if ~isa(epsilon_zx, Nothing)

matrix_epsilon_zx = spdiagm(nx_Ey*ny_Ez*nz_Ey, nx_Ey*ny_Ez*nz_Ey, epsilon_zx[:])

A_off_diagonal_epsilon[(nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez), 1:nt_Ex] = A_off_diagonal_epsilon[(nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez), 1:nt_Ex] - (k0dx)^2 * kron(avg_z_Ez', kron(sparse(I,ny_Ez,ny_Ez), sparse(I,nx_Ez,nx_Ez)))*matrix_epsilon_zx*kron(sparse(I,nz_Ex,nz_Ex), kron(sparse(I,ny_Ex,ny_Ex), avg_x_Ex)) # avg_z*epsilon_zx*avg_x*Ex

end

if ~isa(epsilon_zy, Nothing)

matrix_epsilon_zy = spdiagm(nx_Ez*ny_Ex*nz_Ex, nx_Ez*ny_Ex*nz_Ex, epsilon_zy[:])

A_off_diagonal_epsilon[(nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez), (nt_Ex 1):(nt_Ex nt_Ey)] = A_off_diagonal_epsilon[(nt_Ex nt_Ey 1):(nt_Ex nt_Ey nt_Ez), (nt_Ex 1):(nt_Ex nt_Ey)] - (k0dx)^2 * kron(avg_z_Ez', kron(sparse(I,ny_Ez,ny_Ez), sparse(I,nx_Ez,nx_Ez)))*matrix_epsilon_zy*kron(sparse(I,nz_Ey,nz_Ey), kron(avg_y_Ey, sparse(I,nx_Ey,nx_Ey))) # avg_z*epsilon_zy*avg_y*Ey

end

A = A A_off_diagonal_epsilon

end

if use_UPML

# sx_Ex = sx_Hz, because Ex and Hz share the same x-coordinate (and sx) inside same Yee-cell.

# The same reason applies to others.

sx_Ex = sx_Hz

sy_Ey = sy_Hx

sz_Ez = sz_Hy

sy_Ex = sy_Ez

sz_Ey = sz_Ex

sx_Ez = sx_Ey

# Construct a 3D S_E to transform between A_SCPML to A_UPML

(sx_Ex_3d,sy_Ex_3d,sz_Ex_3d) = ndgrid(sx_Ex,sy_Ex,sz_Ex)

(sx_Ey_3d,sy_Ey_3d,sz_Ey_3d) = ndgrid(sx_Ey,sy_Ey,sz_Ey)

(sx_Ez_3d,sy_Ez_3d,sz_Ez_3d) = ndgrid(sx_Ez,sy_Ez,sz_Ez)

S_E = spdiagm(nt_Ex nt_Ey nt_Ez, nt_Ex nt_Ey nt_Ez, vcat((sy_Ex_3d.*sz_Ex_3d./sx_Ex_3d)[:], (sx_Ey_3d.*sz_Ey_3d./sy_Ey_3d)[:], (sx_Ez_3d.*sy_Ez_3d./sz_Ez_3d)[:]))

A = S_E*A

end

end

# determine the symmetry of matrix A, assuming no spatial symmetry in epsilon or inv_epsilon

# 2D TE fields will be added later

is_symmetric_A = true

if (~use_2D_TM && isa(xBC, Number) && xBC != 0 && xBC != pi && (nx_Ex > 1 || ny_Ex > 1 || nz_Ex > 1)) ||

(isa(yBC, Number) && yBC != 0 && yBC != pi && ((use_2D_TM && ny_Ex > 1) || (~use_2D_TM && (nx_Ey > 1 || ny_Ey > 1 || nz_Ey > 1)))) ||

(isa(zBC, Number) && zBC != 0 && zBC != pi && ((use_2D_TM && nz_Ex > 1) || (~use_2D_TM && (nx_Ez > 1 || ny_Ez > 1 || nz_Ez > 1))))

# Bloch periodic boundary condition with k_B*periodicity != 0 or pi breaks the symmetry of A because its ddx is complex-valued, except when there is only one pixel.

is_symmetric_A = false

elseif (~use_2D_TM && (xPML[1].npixels != 0 || xPML[2].npixels != 0 || yPML[1].npixels != 0 || yPML[2].npixels != 0 || zPML[1].npixels != 0 || zPML[2].npixels != 0)) || (use_2D_TM && ~use_UPML && (yPML[1].npixels != 0 || yPML[2].npixels != 0 || zPML[1].npixels != 0 || zPML[2].npixels != 0))

is_symmetric_A = false

end

if ~use_2D_TM

return (A, is_symmetric_A, xPML, yPML, zPML)

else

return (A, is_symmetric_A, yPML, zPML)

end

end

"""

MESTI_BUILD_FDFD_MATRIX The finite-difference frequency-domain operator for 2D TM waves.

A = mesti_build_fdfd_matrix(epsilon_xx, k0dx, yBC, zBC, yPML, zPML, use_UPML)

returns A as a sparse matrix representing wave operator [- (d/dy)^2 - (d/dz)^2 - k0^2*epsilon(y,z)]*(dx^2)

for the Ex(y,z) component of transverse-magnetic (TM) fields (Ex, Hy, Hz),

discretized on a square grid with grid size dx through center difference.

Matrix A has size [ny_Ex*nz_Ex, ny_Ex*nz_Ex].

Here are the differences between 3D and 2D TM versions

epsilon_xx (matrix; required):

epsilon_xx is a ny_Ex-by-nz_Ex matrix (real or complex)

discretizing the relative permittivity profile on Ex-sites. ny_Ex/nz_Ex

is the total number of grids for Ex in y/z-direction.

A (sparse matrix):

nt_Ex-by-nt_Ex sparse matrix representing the 2D FDFD operator.

nt_Ex = ny_Ex*nz_Ex is the total number of grids for Ex.

=== Input Arguments ===

epsilon_xx (matrix; required):

epsilon_xx is a ny_Ex-by-nz_Ex matrix (real or complex)

discretizing the relative permittivity profile epsilon_xx on Ex-sites.

The notation ni_Ej is the total number of grids for Ej in i-direction.

k0dx (numeric scalar, real or complex; required):

Normalized frequency k0*dx = (2*pi/vacuum_wavelength)*dx.

yBC (character vector or numeric scalar; required):

Boundary condition (BC) at the two ends in y direction, effectively

Ex(n,m,l) at m=0 and m=ny_Ex 1.

one pixel beyond the computation domain. Available choices are:

"periodic" - Ex(m ny_Ex,l) = Ex(m,l),

"PEC" - Ex(0,l) = Ex(ny_Ex,l) = 0,

"PMC" - Ex(0,l) = Ex(1,l); Ex(ny_Ex 1,l) = Ex(ny_Ex,l),

"PECPMC" - Ex(0,l) = 0; Ex(ny_Ex 1,m,l) = Ex(ny_Ex,m,l),

"PMCPEC" - Ex(0,l) = Ex(1,l); Ex(ny_Ex 1,m,l) = 0,

where PEC stands for perfect electric conductor and PMC stands for perfect

magnetic conductor.

When yBC is a numeric scalar, the Bloch periodic boundary condition is

used with Ex(m ny_Ex,l) = Ex(m,l)*exp(1i*yBC); in other words,

yBC = ky_B*ny_Ex*dx = ky_B*Lambda where ky_B is the Bloch wave number and

Lambda = ny_Ex*dx is the periodicity in y.

zBC (character vector or numeric scalar; required):

Boundary condition in z direction, analogous to yBC.

yPML (a two-element PML vector; optional):

Parameters for perfectly matched layer (PML) in y direction.

yPML = [PML_left, PML_right]

If users do not want to put PML on either side, just set the field "npixels" = 0.

For example, PML_left.npixels = 0: no PML on the left.

In each case, PML is a PML structure with the following fields:

npixels (non-negative integer scalar; required): Number of PML pixels.

Note this is within syst.epsilon or syst.inv_epsilon.

power_sigma (non-negative scalar; optional): Power of the polynomial

grading for the conductivity sigma; defaults to 3.

sigma_max_over_omega (non-negative scalar; optional):

Conductivity at the end of the PML; defaults to

0.8*(power_sigma 1)/(k0dx*sqrt(epsilon_bg)).

where epsilon_bg is the average relative permittivity along the

last slice of the PML. This is used to attenuate propagating waves.

power_kappa (non-negative scalar; optional): Power of the polynomial

grading for the real-coordinate-stretching factor kappa; defaults

to 3.

kappa_max (real scalar no smaller than 1; optional):

Real-coordinate-stretching factor at the end of the PML; defaults

to 15. This is used to accelerate the attenuation of evanescent

waves. kappa_max = 1 means no real-coordinate stretching.

power_alpha (non-negative scalar; optional): Power of the polynomial

grading for the CFS alpha factor; defaults to 1.

alpha_max_over_omega (non-negative scalar; optional): Complex-

frequency-shifting (CFS) factor at the beginning of the PML. This

is typically used in time-domain simulations to suppress late-time

(low-frequency) reflections. We do not use it by default

(alpha_max_over_omega = 0) since we are in frequency domain.

We use the following PML coordinate-stretching factor:

s(p) = kappa(p) sigma(p)./(alpha(p) - i*omega)

with

sigma(p)/omega = sigma_max_over_omega*(p.^power_sigma),

kappa(p) = 1 (kappa_max-1)*(p.^power_kappa),

alpha(p)/omega = alpha_max_over_omega*((1-p).^power_alpha),

where omega is frequency, and p goes linearly from 0 at the beginning of

the PML to 1 at the end of the PML.

zPML (a two-element vector; optional):

Parameters for PML in z direction, analogous to yPML.

use_UPML (logical scalar; optional, defaults to true):

Whether to use uniaxial PML (UPML) or not. If not, stretched-coordinate

PML (SC-PML) will be used.

=== Output Arguments ===

A (sparse matrix):

nt_Ex-by-nt_Ex sparse matrix representing the 2D FDFD operator.

nt_Ex = ny_Ex*nz_Ex is the total number of grids for Ex.

is_symmetric_A (logical scalar):

Whether matrix A is symmetric or not.

yPML (two-element vector):

PML parameters used on the low and high sides of y direction, if any.

zPML (two-element vector):

PML parameters used on the low and high sides of z direction, if any.

"""

function mesti_build_fdfd_matrix(epsilon_xx::Union{Matrix{Int64},Matrix{Float64},Matrix{ComplexF64}}, k0dx::Union{Real,Complex}, yBC::Union{String,Real,Complex}, zBC::Union{String,Real,Complex}, yPML::Vector{PML} = [PML(0), PML(0)], zPML::Vector{PML} = [PML(0), PML(0)], use_UPML::Bool=true)

return mesti_build_fdfd_matrix(epsilon_xx, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, k0dx, nothing, yBC, zBC, nothing, yPML, zPML, use_UPML)

end

# When only the diagonal terms of epsilon (i.e. epsilon_xx, epsilon_yy, epsilon_zz) are specified in the 3D case.

function mesti_build_fdfd_matrix(epsilon_xx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Matrix{Int64},Matrix{Float64},Matrix{ComplexF64}}, epsilon_yy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_zz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, k0dx::Union{Real,Complex}, xBC::Union{String,Real,Complex,Nothing}, yBC::Union{String,Real,Complex}, zBC::Union{String,Real,Complex}, xPML::Union{Vector{PML},Nothing} = [PML(0), PML(0)], yPML::Vector{PML} = [PML(0), PML(0)], zPML::Vector{PML} = [PML(0), PML(0)], use_UPML::Bool=true)

return mesti_build_fdfd_matrix(epsilon_xx, nothing, nothing, nothing, epsilon_yy, nothing, nothing, nothing, epsilon_zz, k0dx, xBC, yBC, zBC, xPML, yPML, zPML, use_UPML)

end

"""

CHECK_BC_AND_GRID(BC, n_Ex, n_Ey, n_Ez, direction) is a helper function and checks the consistency of boundary condition and number of grids in n_Ex, n_Ey, and n_Ez.

"""

function check_BC_and_grid(BC::Union{String,Real,Complex}, n_Ex::Int, n_Ey::Int, n_Ez::Int, direction::String)

if isa(BC, Number)

if (n_Ex != n_Ez || n_Ex != n_Ey)

throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx, epsilon_yy, and epsilon_zz should be the same for Bloch periodic boundary condition in $(direction)-direction."))

end

elseif (n_Ex != n_Ez || n_Ex != n_Ey) && (BC == "periodic" || BC == "PECPMC" || BC == "PMCPEC")

throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx, epsilon_yy, and epsilon_zz should be the same for $(BC) boundary condition in $(direction)-direction."))

else

if direction == "x"

if (n_Ex != n_Ez 1 || n_Ex != n_Ey 1) && (BC == "PEC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx should be 1 grid more than epsilon_yy (epsilon_zz) for $(BC) boundary condition in $(direction)-direction.")); end

if (n_Ex != n_Ez-1 || n_Ex != n_Ey-1) && (BC == "PMC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx should be 1 grid less than epsilon_yy (epsilon_zz) for $(BC) boundary condition in $(direction)-direction.")); end

elseif direction == "y"

if (n_Ey != n_Ex 1 || n_Ey != n_Ez 1) && (BC == "PEC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_yy should be 1 grid more than epsilon_xx (epsilon_zz) for $(BC) boundary condition in $(direction)-direction.")); end

if (n_Ey != n_Ex-1 || n_Ey != n_Ez-1) && (BC == "PMC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_yy should be 1 grid less than epsilon_xx (epsilon_zz) for $(BC) boundary condition in $(direction)-direction.")); end

elseif direction == "z"

if (n_Ez != n_Ey 1 || n_Ez != n_Ex 1) && (BC == "PEC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_zz should be 1 grid more than epsilon_xx (epsilon_yy) for $(BC) boundary condition in $(direction)-direction.")); end

if (n_Ez != n_Ey-1 || n_Ez != n_Ex-1) && (BC == "PMC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_zz should be 1 grid less than epsilon_xx (epsilon_yy) for $(BC) boundary condition in $(direction)-direction.")); end

end

end

end

"""

CHECK_BC_AND_GRID(BC, n_Ex, n_Ey, direction) is a helper function and checks the consistency of boundary condition and number of grids in n_Ex, and n_Ey.

"""

function check_BC_and_grid(BC::Union{String,Real,Complex}, n_Ex::Int, n_Ey::Int, direction::String)

if isa(BC, Number)

if (n_Ex != n_Ey)

throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx, epsilon_yy, should be the same for Bloch periodic boundary condition in $(direction)-direction."))

end

elseif (n_Ex != n_Ey) && (BC == "periodic" || BC == "PECPMC" || BC == "PMCPEC")

throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx, epsilon_yy, should be the same for $(BC) boundary condition in $(direction)-direction."))

else

if direction == "x"

if (n_Ex != n_Ey 1) && (BC == "PEC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx should be 1 grid more than epsilon_yy for $(BC) boundary condition in $(direction)-direction.")); end

if (n_Ex != n_Ey-1) && (BC == "PMC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_xx should be 1 grid less than epsilon_yy for $(BC) boundary condition in $(direction)-direction.")); end

elseif direction == "y"

if (n_Ey != n_Ex 1) && (BC == "PEC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_yy should be 1 grid more than epsilon_xx for $(BC) boundary condition in $(direction)-direction.")); end

if (n_Ey != n_Ex-1) && (BC == "PMC"); throw(ArgumentError("Number of grids along $(direction) provided by epsilon_yy should be 1 grid less than epsilon_xx for $(BC) boundary condition in $(direction)-direction.")); end

end

end

end

"""

CONVERT_BC is a helper function to handles the case of string in BC

"""

function convert_BC(BC::Union{String,Real,Complex},direction::String)

if isa(BC, Number)

return BC

elseif lowercase(BC) == "pec"

return "PEC"

elseif lowercase(BC) == "pmc"

return "PMC"

elseif lowercase(BC) == "pecpmc"

return "PECPMC"

elseif lowercase(BC) == "pmcpec"

return "PMCPEC"

elseif lowercase(BC) == "periodic"

return "periodic"

elseif lowercase(BC) == "bloch"

throw(ArgumentError("To use Bloch periodic boundary condition in $(direction)-direction, set the input argument $(direction)BC to k$(direction)_B*p_$(direction) where k$(direction)_B is the Bloch wave number and Lambda_$(direction) is the periodicity along $(direction)-direction."))

else

ethrow(ArgumentError("Input argument $(direction)BC = \"$(BC)\" is not a supported option."))

end

end

"""

BUIlD_DDX_E builds the first-derivative matrix

The following documentation all take direction = x for example and can be analogous to other direction.

ddx: x-direction first-derivative matrix (acting on Ey or Ez)

avg: average matrix (acting on Ey, or Ez)

s_E: x-coordinate-stretching factor for Ey, or Ez (on integer sites)

s_H: x-coordinate-stretching factor for Hy, or Hz (on half-integer sites)

n_E: number of sites in x for Ey or Ez

ind_PML_E: a vector containing indices of the PML pixels for Ey or Ez (on integer sites)

"""

function build_ddx_E(n_E::Int, BC::Union{String,Real,Complex}, pml::Vector{PML}, direction::String)

# Handle periodic and Bloch periodic boundary (adding the feature later) conditions

if isa(BC, Number)

kLambda = BC

BC = "Bloch"

if ~isa(kLambda, Real)

@warn "k$(direction)_B*periodicity = $(real(kLambda)) 1im*$(imag(kLambda)) is a complex number."

end

elseif BC == "periodic"

kLambda = 0

BC = "Bloch"

end

# Build the first-derivative matrix and the average matrix on E

# p = [p(1), ..., p(n_E)].' with p = Ey or Ez, on integer sites

# dp = (dp/dx)*dx, proportional to Hy or Hz, on half-integer sites

# avg_p = average of p between two neighboring sites, on half-integer sites

if BC == "Bloch"

# p(n_E 1) = p(1)*exp(1i*kLambda); p(0) = p(n_E)*exp(-1i*kLambda)

# ddx*p = [dp(1.5), ..., dp(n_E 0.5)].'

ddx = spdiagm(n_E, n_E, 1 => ones(n_E-1), 0 => -ones(n_E), 1-n_E => exp(1im*kLambda)*ones(1))

# avg*p = [avg_p(1.5), ..., avg_p(n_E 0.5)].'

avg = spdiagm(n_E, n_E, 1 => ones(n_E-1)/2, 0 => ones(n_E)/2, 1-n_E => exp(1im*kLambda)*ones(1)/2)

elseif BC == "PEC" # PEC on both sides

# p(0) = p(n_E 1) = 0

# ddx*p = dp = [dp(0.5), ..., dp(n 0.5)].'

ddx = spdiagm(n_E 1, n_E, 0 => ones(n_E), -1 => -ones(n_E))

# avg*p = [avg_p(0.5), ..., avg_p(n_E 0.5)].'

avg = spdiagm(n_E 1, n_E, 0 => ones(n_E)/2, -1 => ones(n_E)/2)

elseif BC == "PMC" # PMC on both sides

# p(0) = p(1); p(n_E 1) = p(n_E)

# ddx*p = [dp(1.5), ..., dp(n_E-0.5)].'; exclude dp(0.5) and dp(n_E 0.5) because they are zero

ddx = spdiagm(n_E-1, n_E, 1 => ones(n_E-1), 0 => -ones(n_E-1))

# avg*p = [avg_p(1.5), ..., avg_p(n_E-0.5)].'; exclude avg_p(0.5) and avg_p(n_E 0.5)

avg = spdiagm(n_E-1, n_E, 1 => ones(n_E-1)/2, 0 => ones(n_E-1)/2)

elseif BC == "PECPMC" # PEC on the low side, PMC on the high side

# p(0) = 0; p(n_E 1) = p(n_E)

# ddx*p = [dp(0.5), ..., dp(n_E-0.5)].'; exclude dp(n_E 0.5) because it is zero

ddx = spdiagm(n_E, n_E, 0 => ones(n_E), -1 => -ones(n_E-1))

# avg*p = [avg_p(0.5), ..., avg_p(n_E-0.5)].'; we exclude avg_p(n_E 0.5)

avg = spdiagm(n_E, n_E, 0 => ones(n_E)/2, -1 => ones(n_E-1)/2)

elseif BC == "PMCPEC" # PMC on the low side, PEC on the high side

# p(0) = p(1); p(n_E 1) = 0

# ddx*p = [dp(1.5), ..., dp(n_E 0.5)].'; exclude dp(0.5) because it is zero

ddx = spdiagm(n_E, n_E, 1 => ones(n_E-1), 0 => -ones(n_E))

# avg*p = [avg_p(1.5), ..., avg_p(n_E 0.5)].'; exclude avg_p(0.5)

avg = spdiagm(n_E, n_E, 1 => ones(n_E-1)/2, 0 => ones(n_E)/2)

else

throw(ArgumentError("Input argument $(direction)BC = \"$(BC)\" is not a supported option."))

end

n_H = size(ddx,1) # number of sites for dp/dx (ie, Hy or Hz)

# Coordinate-stretching factor: s(x) = kappa(x) sigma(x)/(alpha(u) - i*omega)

s_E = ones(ComplexF64, n_E) # s-factor for Ey or Ez (on integer sites)

s_H = ones(ComplexF64, n_H) # s-factor for Hy or Hz (on half-integer sites)

# no coordinate stretching if no PML is specified

if (length(pml) == 2 && pml[1].npixels == 0 && pml[2].npixels == 0)

ind_PML_E = [nothing, nothing]

return ddx, avg, s_E, s_H, ind_PML_E

end

# Number of PML pixels on the low and high sides

npixels = [pml[1].npixels, pml[2].npixels]

n = n_E # Also works for 2D TM fields

# Will add 2D TE case later

# Cannot have more PML pixels than the number of pixels

# In our implementation, the conductivity goes to zero at one pixel before PML. So there needs to be at least one more pixel in addition to PML.

if sum(npixels) >= n

throw(ArgumentError("Total number of pixels = $(n) in $(direction) direction must be greater than the number of PML pixels = $(npixels[1]) $(npixels[2]) = $(sum(npixels)) but is not."))

end

# Below, p(x) is a function that goes linearly from p(x)=0 at one site before PML to p(x)=1 at the "end of PML".

# Note that where the "end of PML" is depends on the boundary condition.

# Let index i=0 be one site before PML, so p(i=0)=0.

# Then, i=1 is the first site of PML and i=npixels is the last site of PML we explicitly simulate. But we do not set p(i=npixels)=1.

# For Dirichlet BC, the BC is such that p=0 at i=npixels 1, so we let the end of PML be i=npixels 1, with an effective PML thickness of npixels 1 pixels.

# For Neumann BC, the BC is such that dp=0 at i=npixels 0.5, so we let the end of PML be i=npixels 0.5, with an effective PML thickness of npixels 0.5 pixels.

# For periodic and Bloch periodic BC, we let p(x) be symmetric on the two sides of p with p(0.5)=p(n 0.5)=1; note that p(0.5) and p(n 0.5) are the same site (with a possible Bloch phase difference). If the PML parameters on the two sides are the same, the s-factor will be continuous across the periodic boundary.

# Construct p(x) and their corresponding indices for dp/dx

# Also works for 2D TM fields; will add 2D TE case later

if BC == "Bloch"

if npixels[1]*npixels[2] == 0 && (npixels[1] != 0 || npixels[2] != 0)

@warn "Bloch periodic boundary condition is used with a single-sided PML in $(direction) direction; transmission through PML will only undergo single-pass attenuation."

end

npixels_effective = [npixels[1] 0.5, npixels[2] 0.5]

# p = Ez; ddx*p = [dp(1.5), ..., dp(n_E 0.5)].'

# no s-factor for dp(0.5) since we only consider dp(n_E 0.5)

p_PML_2 = [((1:npixels[1]).-0.5)/npixels_effective[1], ((1:(npixels[2] 1)).-0.5)/npixels_effective[2]]

ind_PML_2 = [reverse(1:npixels[1]), (n 1).-reverse(1:(npixels[2] 1))]

elseif BC == "PEC" # PEC on both sides

npixels_effective = [npixels[1] 1, npixels[2] 1]

p_PML_2 = [((1:(npixels[1] 1)).-0.5)/npixels_effective[1], ((1:(npixels[2] 1)).-0.5)/npixels_effective[2]]

ind_PML_2 = [reverse(1:(npixels[1] 1)), (n 2).-reverse(1:(npixels[2] 1))]

elseif BC == "PMC" # PMC on both sides

npixels_effective = [npixels[1] 0.5, npixels[2] 0.5]

p_PML_2 = [((1:npixels[1]).-0.5)/npixels_effective[1], ((1:(npixels[2])).-0.5)/npixels_effective[2]]

ind_PML_2 = [reverse(1:npixels[1]), n.-reverse(1:(npixels[2]))]

elseif BC == "PECPMC" # PEC on the low side, PMC on the high side

npixels_effective = [npixels[1] 1, npixels[2] 0.5]

p_PML_2 = [((1:npixels[1] 1).-0.5)/npixels_effective[1], ((1:(npixels[2])).-0.5)/npixels_effective[2]]

ind_PML_2 = [reverse(1:npixels[1] 1), (n 1).-reverse(1:(npixels[2]))]

elseif BC == "PMCPEC" # PMC on the low side, PEC on the high side

npixels_effective = [npixels[1] 0.5, npixels[2] 1]

p_PML_2 = [((1:npixels[1]).-0.5)/npixels_effective[1], ((1:(npixels[2] 1)).-0.5)/npixels_effective[2]]

ind_PML_2 = [reverse(1:npixels[1]), (n 1).-reverse(1:(npixels[2] 1))]

else

throw(ArgumentError("Input argument $(direction)BC = \"$(BC)\" is not a supported option."))

end

# Construct p(x) and their corresponding indices for f and d^2f/dx^2

p_PML_1 = [(1:npixels[1])/npixels_effective[1], (1:npixels[2])/npixels_effective[2]]

ind_PML_1 = [reverse(1:npixels[1]), (n 1).-reverse(1:npixels[2])]

# Also works for 2D TM modes; will add 2D TE case later

p_PML_E = p_PML_1

ind_PML_E = ind_PML_1

p_PML_H = p_PML_2

ind_PML_H = ind_PML_2

# Loop over PML on the two sides

for ii = 1:2

if npixels[ii] > 0

# coordinate-stretching factor s(p) from Eq 7.73 of Taflove & Hagness's 2005 FDTD book

# sigma is the conductivity, equivalent to imag-coordinate stretching, used to attenuate propagating waves.

# kappa is real-coordinate stretching, used to accelerate the attenuation of evanescent waves.

# alpha is used for complex frequency shifting (CFS) to suppress reflection of low-frequency components for time-domain simulations.

# In general, kappa, sigma, alpha can all be arbitrary functions of position.

# To minimize discretization-induced reflection, kappa should start from 1, and sigma should start from 0.

# We use a polynomial grading for all of them: Eqs 7.60 and 7.79 of Taflove & Hagness's 2005 FDTD book.

kappa(p) = 1 . (pml[ii].kappa_max-1)*(p.^(pml[ii].power_kappa))

sigma_over_omega(p) = pml[ii].sigma_max_over_omega*(p.^(pml[ii].power_sigma))

alpha_over_omega(p) = pml[ii].alpha_max_over_omega*((1 .- p).^(pml[ii].power_alpha))

func_s(p) = kappa(p) . sigma_over_omega(p)./(alpha_over_omega(p) .- 1im)

# Evaluate s(p) on integer and half-integer sites

# Note s_E and s_H should be sampled directly from the polynomial functions; if we use the polynomial for one and linear interpolation for the other, the performance will be much worse.

s_E[Int.(ind_PML_E[ii])] = func_s(p_PML_E[ii]) # column vector

s_H[Int.(ind_PML_H[ii])] = func_s(p_PML_H[ii]) # column vector

end

end

return ddx, avg, s_E, s_H, ind_PML_E

end

"""

BUIlD_AVE_X_EX builds the average matrix that acts on Ex along the x-direction and can be analogous to other direction.

avg: average matrix (acting on Ex along x-direction)

"""

function build_ave_x_Ex(n_E::Int, BC::Union{String,Real,Complex}, direction::String)

# Handle periodic and Bloch periodic boundary (adding the feature later) conditions

if isa(BC, Number)

kLambda = BC

BC = "Bloch"

if ~isa(kLambda, Real)

@warn "k$(direction)_B*periodicity = $(real(kLambda)) 1im*$(imag(kLambda)) is a complex number."

end

elseif BC == "periodic"

kLambda = 0

BC = "Bloch"

end

# Build the first-derivative matrix and the average matrix on E

# p = [p(1), ..., p(n_E)].' with p = Ey or Ez, on integer sites

# dp = (dp/dx)*dx, proportional to Hy or Hz, on half-integer sites

# avg_p = average of f between two neighboring sites, on half-integer sites

if BC == "Bloch"

# p(n_E 1) = p(1)*exp(1i*kLambda); p(0) = p(n_E)*exp(-1i*kLambda)

# avg*p = [avg_p(1.5), ..., avg_p(n_E 0.5)].'

avg = spdiagm(n_E, n_E, 1 => ones(n_E-1)/2, 0 => ones(n_E)/2, 1-n_E => exp(1im*kLambda)*ones(1)/2)

elseif BC == "PEC" # PEC on both sides

# p(0) = p(1); p(n_E 1) = p(n_E)

# avg*p = [avg_p(1.5), ..., avg_p(n_E-0.5)].'; exclude avg_p(0.5) and avg_p(n_E 0.5)

avg = spdiagm(n_E-1, n_E, 1 => ones(n_E-1)/2, 0 => ones(n_E-1)/2)

elseif BC == "PMC" # PMC on both sides

# p(0) = p(n_E 1) = 0

# avg*p = [avg_p(0.5), ..., avg_p(n_E 0.5)].'

avg = spdiagm(n_E 1, n_E, 0 => ones(n_E)/2, -1 => ones(n_E)/2)

elseif BC == "PECPMC" # PEC on the low side, PMC on the high side

# p(0) = p(1); p(n_E 1) = 0

# avg*p = [avg_p(1.5), ..., avg_p(n_E 0.5)].'; exclude avg_p(0.5)

avg = spdiagm(n_E, n_E, 1 => ones(n_E-1)/2, 0 => ones(n_E)/2)

elseif BC == "PMCPEC" # PMC on the low side, PEC on the high side

# p(0) = 0; p(n_E 1) = p(n_E)

# avg*p = [avg_p(0.5), ..., avg_p(n_E-0.5)].'; we exclude avg_p(n_E 0.5)

avg = spdiagm(n_E, n_E, 0 => ones(n_E)/2, -1 => ones(n_E-1)/2)

else

throw(ArgumentError("Input argument $(direction)BC = \"$(BC)\" is not a supported option."))

end

return avg

end