text:
export HelmholtzMap, HelmholtzPreconditioner, ModeSolver, update_k, update_k!
export update_ε⁻¹, ε⁻¹_ω, mag_m_n, mag_m_n2, mag_m_n!, mag_mn, kx_ct, kx_tc, zx_tc, zx_ct
export ε⁻¹_dot, ε⁻¹_dot_t, _M!, _P!, kx_ct!, kx_tc!, zx_tc!, kxinv_ct!
export kxinv_tc!, ε⁻¹_dot!, ε_dot_approx!, HMₖH, HMH, tc, ct, ng_z, eid!

"""
################################################################################
#																			   #
#			Function Definitions Implementing Non-Mutating Operators		   #
#																			   #
################################################################################
"""

# 3D

"""
    tc: v⃗ (transverse vector) → a⃗ (cartesian vector)
"""
function tc(H::AbstractArray{T,4},mn) where T<:Union{Real,Complex}
    @tullio h[a,i,j,k] := H[b,i,j,k] * mn[a,b,i,j,k]
end

"""
    ct: a⃗ (cartesian vector) → v⃗ (transverse vector)
"""
function ct(h::AbstractArray{T,4},mn) where T<:Union{Real,Complex}
    @tullio H[a,i,j,k] := h[b,i,j,k] * mn[b,a,i,j,k]
end

"""
    kx_tc: a⃗ (cartesian vector) = k⃗ × v⃗ (transverse vector)
"""
function kx_tc(H::AbstractArray{T,4},mn,mag) where T
	kxscales = [-1.; 1.]
	kxinds = [2; 1]
    @tullio d[a,i,j,k] := kxscales[b] * H[kxinds[b],i,j,k] * mn[a,b,i,j,k] * mag[i,j,k] nograd=(kxscales,kxinds) # fastmath=false
	# @tullio d[a,i,j,k] := H[2,i,j,k] * m[a,i,j,k] * mag[i,j,k] - H[1,i,j,k] * n[a,i,j,k] * mag[i,j,k]  # nograd=(kxscales,kxinds) fastmath=false
end

"""
    kx_c2t: v⃗ (transverse vector) = k⃗ × a⃗ (cartesian vector)
"""
function kx_ct(e⃗::AbstractArray{T,4},mn,mag) where T
	# mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])))
	kxscales = [-1.; 1.]
    kxinds = [2; 1]
    @tullio H[b,i,j,k] := kxscales[b] * e⃗[a,i,j,k] * mn[a,kxinds[b],i,j,k] * mag[i,j,k] nograd=(kxinds,kxscales) # fastmath=false
end

"""
    zx_t2c: a⃗ (cartesian vector) = ẑ × v⃗ (transverse vector)
"""
function zx_tc(H::AbstractArray{T,4},mn) where T
	zxinds = [2; 1; 3]
	zxscales = [-1.; 1.; 0.]
	@tullio zxH[a,i,j,k] := zxscales[a] * H[b,i,j,k] * mn[zxinds[a],b,i,j,k] nograd=(zxscales,zxinds) # fastmath=false
end

"""
    zx_c2t: v⃗ (transverse vector) = ẑ × a⃗ (cartesian vector)
"""
function zx_ct(e⃗::AbstractArray{T,4},mn) where T
	zxinds = [2; 1; 3]
	zxscales = [-1.; 1.; 0.]
	@tullio zxe⃗[b,i,j,k] := zxscales[a] * e⃗[a,i,j,k] * mn[b,zxinds[a],i,j,k] nograd=(zxscales,zxinds) # fastmath=false
end

"""
    ε⁻¹_dot_t: e⃗  = ε⁻¹ ⋅ d⃗ (transverse vectors)
"""
function ε⁻¹_dot_t(d⃗::AbstractArray{T,4},ε⁻¹) where T
	# eif = flat(ε⁻¹)
	@tullio e⃗[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * fft(d⃗,(2:4))[b,i,j,k]  #fastmath=false
	return ifft(e⃗,(2:4))
end

"""
    ε⁻¹_dot: e⃗  = ε⁻¹ ⋅ d⃗ (cartesian vectors)
"""
function ε⁻¹_dot(d⃗::AbstractArray{T,4},ε⁻¹) where T
	# eif = flat(ε⁻¹)
	@tullio e⃗[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * d⃗[b,i,j,k]  #fastmath=false
end

function HMₖH(H::AbstractArray{Complex{T},4},ε⁻¹,mag,m,n)::T where T<:Real
	mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])))
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( zx_tc(H,mn), (2:4) ), real(flat(ε⁻¹))), (2:4)),mn,mag) ) )
end

function HMₖH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,3},m::AbstractArray{T,4},n::AbstractArray{T,4})::T where T<:Real
	Nx,Ny,Nz = size(mag)
	Ha = reshape(H,(2,Nx,Ny,Nz))
	HMₖH(Ha,ε⁻¹,mag,m,n)
end

function HMₖH(H::AbstractArray{Complex{T},4},ε⁻¹,mag,mn)::T where T<:Real
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( zx_tc(H,mn), (2:4) ), real(flat(ε⁻¹))), (2:4)),mn,mag) ) )
end

function HMₖH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,3},mn::AbstractArray{T,5})::T where T<:Real
	Nx,Ny,Nz = size(mag)
	Ha = reshape(H,(2,Nx,Ny,Nz))
	HMₖH(Ha,ε⁻¹,mag,mn)
end

function HMH(H::AbstractArray{Complex{T},4},ε⁻¹,mag,m,n)::T where T<:Real
	mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])))
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( kx_tc(H,mn,mag), (2:4) ), real(flat(ε⁻¹))), (2:4)),mn,mag) ) )
end

function HMH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,3},m::AbstractArray{T,4},n::AbstractArray{T,4})::T where T<:Real
	Nx,Ny,Nz = size(mag)
	Ha = reshape(H,(2,Nx,Ny,Nz))
	HMH(Ha,ε⁻¹,mag,m,n)
end

function HMH(H::AbstractArray{Complex{T},4},ε⁻¹,mag,mn)::T where T<:Real
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( kx_tc(H,mn,mag), (2:4) ), real(flat(ε⁻¹))), (2:4)),mn,mag) ) )
end

function HMH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,3},m::AbstractArray{T,5})::T where T<:Real
	Nx,Ny,Nz = size(mag)
	Ha = reshape(H,(2,Nx,Ny,Nz))
	HMH(Ha,ε⁻¹,mag,mn)
end

function ng_z(Hₜ::AbstractArray{Complex{T},4},ω,ε⁻¹,nng,mag,m,n)::T where T<:Real
	mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])))
	E = 1im * ε⁻¹_dot( fft( kx_tc(Hₜ,mn,mag), (2:4) ), real(ε⁻¹))
	H = (-1im * ω) * fft( tc(Hₜ,mn), (2:4) )
	W = real(dot(E,_dot(nng,E))) + ( ω^2 * size(H,2) * size(H,3) * size(H,4) )
	@tullio P_z := conj(E)[1,ix,iy,iz] * H[2,ix,iy,iz] - conj(E)[2,ix,iy,iz] * H[1,ix,iy,iz]
	return W / (2*real(P_z))
end

function ng_z(Hₜ::AbstractVector{Complex{T}},ω,ε⁻¹,nng,mag::AbstractArray{T,3},m::AbstractArray{T,4},n::AbstractArray{T,4})::T where T<:Real
	Nx,Ny,Nz = size(mag)
	Ha = reshape(Hₜ,(2,Nx,Ny,Nz))
	ng_z(Ha,ω,ε⁻¹,nng,mag,m,n)
end

function ng_z(Hₜ::AbstractArray{Complex{T},4},ω,ε⁻¹,nng,mag,mn)::T where T<:Real
	E = 1im * ε⁻¹_dot( fft( kx_tc(Hₜ,mn,mag), (2:4) ), real(ε⁻¹))
	H = (-1im * ω) * fft( tc(Hₜ,mn), (2:4) )
	W = real(dot(E,_dot(nng,E))) + ( ω^2 * size(H,2) * size(H,3) * size(H,4) )
	@tullio P_z := conj(E)[1,ix,iy,iz] * H[2,ix,iy,iz] - conj(E)[2,ix,iy,iz] * H[1,ix,iy,iz]
	return W / (2*real(P_z))
end

function ng_z(Hₜ::AbstractVector{Complex{T}},ω,ε⁻¹,nng,mag::AbstractArray{T,3},m::AbstractArray{T,5})::T where T<:Real
	Nx,Ny,Nz = size(mag)
	Ha = reshape(Hₜ,(2,Nx,Ny,Nz))
	ng_z(Ha,ω,ε⁻¹,nng,mag,mn)
end

# 2D

# """
#     tc: v⃗ (transverse vector) → a⃗ (cartesian vector)
# """
# function tc(H::AbstractArray{T,3},mn) where T<:Union{Real,Complex}
#     @tullio h[a,i,j] := H[b,i,j] * mn[a,b,i,j]
# end
#
# """
#     ct: a⃗ (cartesian vector) → v⃗ (transverse vector)
# """
# function ct(h::AbstractArray{T,3},mn) where T<:Union{Real,Complex}
#     @tullio H[a,i,j] := h[b,i,j] * mn[b,a,i,j]
# end
#
# """
#     kx_tc: a⃗ (cartesian vector) = k⃗ × v⃗ (transverse vector)
# """
# function kx_tc(H::AbstractArray{T,3},mn,mag) where T
# 	kxscales = [-1.; 1.]
# 	kxinds = [2; 1]
#     @tullio d[a,i,j] := kxscales[b] * H[kxinds[b],i,j] * mn[a,b,i,j] * mag[i,j] nograd=(kxscales,kxinds) # fastmath=false
# 	# @tullio d[a,i,j] := H[2,i,j] * m[a,i,j] * mag[i,j] - H[1,i,j] * n[a,i,j] * mag[i,j]  # nograd=(kxscales,kxinds) fastmath=false
# end
#
# """
#     kx_c2t: v⃗ (transverse vector) = k⃗ × a⃗ (cartesian vector)
# """
# function kx_ct(e⃗::AbstractArray{T,3},mn,mag) where T
# 	# mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])))
# 	kxscales = [-1.; 1.]
#     kxinds = [2; 1]
#     @tullio H[b,i,j] := kxscales[b] * e⃗[a,i,j] * mn[a,kxinds[b],i,j] * mag[i,j] nograd=(kxinds,kxscales) # fastmath=false
# end
#
# """
#     zx_t2c: a⃗ (cartesian vector) = ẑ × v⃗ (transverse vector)
# """
# function zx_tc(H::AbstractArray{T,3},mn) where T
# 	zxinds = [2; 1; 3]
# 	zxscales = [-1.; 1.; 0.]
# 	@tullio zxH[a,i,j] := zxscales[a] * H[b,i,j] * mn[zxinds[a],b,i,j] nograd=(zxscales,zxinds) # fastmath=false
# end
#
# """
#     zx_c2t: v⃗ (transverse vector) = ẑ × a⃗ (cartesian vector)
# """
# function zx_ct(e⃗::AbstractArray{T,3},mn) where T
# 	zxinds = [2; 1; 3]
# 	zxscales = [-1.; 1.; 0.]
# 	@tullio zxe⃗[b,i,j] := zxscales[a] * e⃗[a,i,j] * mn[b,zxinds[a],i,j] nograd=(zxscales,zxinds) threads=false # fastmath=false
# end

"""
    tc: v⃗ (transverse vector) → a⃗ (cartesian vector)
"""
function tc(H::AbstractArray{T,3},mn) where T<:Union{Real,Complex}
    @tullio h[a,i,j] := H[b,i,j] * mn[a,b,i,j]
end

"""
    ct: a⃗ (cartesian vector) → v⃗ (transverse vector)
"""
function ct(h::AbstractArray{T,3},mn) where T<:Union{Real,Complex}
    @tullio H[a,i,j] := h[b,i,j] * mn[b,a,i,j]
end

"""
    kx_tc: a⃗ (cartesian vector) = k⃗ × v⃗ (transverse vector)
"""
function kx_tc(H::AbstractArray{T,3},mn,mag) where T
	kxscales = [1.; -1.] #[-1.; 1.]
	kxinds = [2; 1]
    @tullio d[a,i,j] := kxscales[b] * H[kxinds[b],i,j] * mn[a,b,i,j] * mag[i,j] nograd=(kxscales,kxinds) # fastmath=false
	# @tullio d[a,i,j] := H[2,i,j] * m[a,i,j] * mag[i,j] - H[1,i,j] * n[a,i,j] * mag[i,j]  # nograd=(kxscales,kxinds) fastmath=false
end

"""
    kx_c2t: v⃗ (transverse vector) = k⃗ × a⃗ (cartesian vector)
"""
function kx_ct(e⃗::AbstractArray{T,3},mn,mag) where T
	# mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3],size(m)[4])))
	kxscales = [1.; -1.] #[-1.; 1.]
    kxinds = [2; 1]
    @tullio H[b,i,j] := kxscales[b] * e⃗[a,i,j] * mn[a,kxinds[b],i,j] * mag[i,j] nograd=(kxinds,kxscales) # fastmath=false
end

"""
    zx_t2c: a⃗ (cartesian vector) = ẑ × v⃗ (transverse vector)
"""
function zx_tc(H::AbstractArray{T,3},mn) where T
	zxinds = [2; 1; 3]
	zxscales = [-1.; 1.; 0.]
	@tullio zxH[a,i,j] := zxscales[a] * H[b,i,j] * mn[zxinds[a],b,i,j] nograd=(zxscales,zxinds) # fastmath=false
end

"""
    zx_c2t: v⃗ (transverse vector) = ẑ × a⃗ (cartesian vector)
"""
function zx_ct(e⃗::AbstractArray{T,3},mn) where T
	zxinds = [2; 1; 3]
	zxscales = [-1.; 1.; 0.]
	@tullio zxe⃗[b,i,j] := zxscales[a] * e⃗[a,i,j] * mn[b,zxinds[a],i,j] nograd=(zxscales,zxinds) threads=false # fastmath=false
end

"""
    ε⁻¹_dot_t: e⃗  = ε⁻¹ ⋅ d⃗ (transverse vectors)
"""
function ε⁻¹_dot_t(d⃗::AbstractArray{T,3},ε⁻¹) where T
	# eif = flat(ε⁻¹)
	@tullio e⃗[a,i,j] :=  ε⁻¹[a,b,i,j] * fft(d⃗,(2:4))[b,i,j]  #fastmath=false
	return ifft(e⃗,(2:4))
end

"""
    ε⁻¹_dot: e⃗  = ε⁻¹ ⋅ d⃗ (cartesian vectors)
"""
function ε⁻¹_dot(d⃗::AbstractArray{T,3},ε⁻¹) where T
	# eif = flat(ε⁻¹)
	@tullio e⃗[a,i,j] :=  ε⁻¹[a,b,i,j] * d⃗[b,i,j]  #fastmath=false
end

function HMₖH(H::AbstractArray{Complex{T},3},ε⁻¹,mag,m,n)::T where T<:Real
	mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3])))
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( zx_tc(H,mn), (2:3) ), real(ε⁻¹)), (2:3)),mn,mag) ) )
end

function HMₖH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,2},m::AbstractArray{T,3},n::AbstractArray{T,3})::T where T<:Real
	Nx,Ny = size(mag)
	Ha = reshape(H,(2,Nx,Ny))
	HMₖH(Ha,ε⁻¹,mag,m,n)
end

function HMₖH(H::AbstractArray{Complex{T},3},ε⁻¹,mag,mn)::T where T<:Real
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( zx_tc(H,mn), (2:3) ), real(ε⁻¹)), (2:3)),mn,mag) ) )
end

function HMₖH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,2},mn::AbstractArray{T,4})::T where T<:Real
	Nx,Ny = size(mag)
	Ha = reshape(H,(2,Nx,Ny))
	HMₖH(Ha,ε⁻¹,mag,mn)
end

function HMH(H::AbstractArray{Complex{T},3},ε⁻¹,mag,m,n)::T where T<:Real
	mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3])))
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( kx_tc(H,mn,mag), (2:3) ), real(ε⁻¹)), (2:3)),mn,mag) ) )
end

function HMH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,2},m::AbstractArray{T,3},n::AbstractArray{T,3})::T where T<:Real
	Nx,Ny = size(mag)
	Ha = reshape(H,(2,Nx,Ny))
	HMH(Ha,ε⁻¹,mag,m,n)
end

function HMH(H::AbstractArray{Complex{T},3},ε⁻¹,mag,mn)::T where T<:Real
	real( dot(H, kx_ct( ifft( ε⁻¹_dot( fft( kx_tc(H,mn,mag), (2:3) ), real(ε⁻¹)), (2:3)),mn,mag) ) )
end

function HMH(H::AbstractVector{Complex{T}},ε⁻¹,mag::AbstractArray{T,2},mn::AbstractArray{T,4})::T where T<:Real
	Nx,Ny = size(mag)
	Ha = reshape(H,(2,Nx,Ny))
	HMH(Ha,ε⁻¹,mag,mn)
end

function ng_z(Hₜ::AbstractArray{Complex{T},3},ω,ε⁻¹,nng,mag,m,n)::T where T<:Real
	mn = vcat(reshape(m,(1,size(m)[1],size(m)[2],size(m)[3])),reshape(n,(1,size(m)[1],size(m)[2],size(m)[3])))
	E = 1im * ε⁻¹_dot( fft( kx_tc(Hₜ,mn,mag), (2:3) ), real(ε⁻¹))
	H = (-1im * ω) * fft( tc(Hₜ,mn), (2:3) )
	W = real(dot(E,_dot(nng,E))) + ( ω^2 * size(H,2) * size(H,3) )
	@tullio P_z := conj(E)[1,ix,iy] * H[2,ix,iy] - conj(E)[2,ix,iy] * H[1,ix,iy]
	return W / (2*real(P_z))
end

function ng_z(Hₜ::AbstractVector{Complex{T}},ω,ε⁻¹,nng,mag::AbstractArray{T,2},m::AbstractArray{T,3},n::AbstractArray{T,3})::T where T<:Real
	Nx,Ny = size(mag)
	Ha = reshape(Hₜ,(2,Nx,Ny))
	ng_z(Ha,ω,ε⁻¹,nng,mag,m,n)
end

function ng_z(Hₜ::AbstractArray{Complex{T},3},ω,ε⁻¹,nng,mag,mn)::T where T<:Real
	E = 1im * ε⁻¹_dot( fft( kx_tc(Hₜ,mn,mag), (2:3) ), real(ε⁻¹))
	H = (-1im * ω) * fft( tc(Hₜ,mn), (2:3) )
	W = real(dot(E,_dot(nng,E))) + ( ω^2 * size(H,2) * size(H,3) )
	@tullio P_z := conj(E)[1,ix,iy] * H[2,ix,iy] - conj(E)[2,ix,iy] * H[1,ix,iy]
	return W / (2*real(P_z))
end

function ng_z(Hₜ::AbstractVector{Complex{T}},ω,ε⁻¹,nng,mag::AbstractArray{T,2},mn::AbstractArray{T,4})::T where T<:Real
	Nx,Ny = size(mag)
	Ha = reshape(Hₜ,(2,Nx,Ny))
	ng_z(Ha,ω,ε⁻¹,nng,mag,mn)
end

"""
################################################################################
#																			   #
#			  Function Definitions Implementing Mutating Operators			   #
#																			   #
################################################################################
"""

# 3D _M! and _P! subroutines

function kx_tc!(d::AbstractArray{Complex{T},4},H::AbstractArray{Complex{T},4},mn::AbstractArray{T,5},mag::AbstractArray{T,3})::AbstractArray{Complex{T},4} where T<:Real
    # @assert size(Y) === size(X)
    # @assert size(d,4) == 3
    # @assert size(H,4) === 2
    @avx for k ∈ axes(d,4), j ∈ axes(d,3), i ∈ axes(d,2), l in 0:0
	# @avx for i ∈ axes(d,1), j ∈ axes(d,2), k ∈ axes(d,3), l in 0:0
		# scale = -mag[i,j,k]
		d[1+l,i,j,k] = ( H[1,i,j,k] * mn[1+l,2,i,j,k] - H[2,i,j,k] * mn[1+l,1,i,j,k] ) * -mag[i,j,k]
        d[2+l,i,j,k] = ( H[1,i,j,k] * mn[2+l,2,i,j,k] - H[2,i,j,k] * mn[2+l,1,i,j,k] ) * -mag[i,j,k]
        d[3+l,i,j,k] = ( H[1,i,j,k] * mn[3+l,2,i,j,k] - H[2,i,j,k] * mn[3+l,1,i,j,k] ) * -mag[i,j,k]
    end
    return d
end

function zx_tc!(d::AbstractArray{Complex{T},4},H::AbstractArray{Complex{T},4},mn::AbstractArray{T,5})::AbstractArray{Complex{T},4} where T<:Real
    @avx for k ∈ axes(d,4), j ∈ axes(d,3), i ∈ axes(d,2), l in 0:0
		d[1+l,i,j,k] = -H[1,i,j,k] * mn[2+l,1,i,j,k] - H[2,i,j,k] * mn[2+l,2,i,j,k]
        d[2+l,i,j,k] =  H[1,i,j,k] * mn[1+l,1,i,j,k] + H[2,i,j,k] * mn[1+l,2,i,j,k]
    end
    return d
end

function kx_ct!(H::AbstractArray{Complex{T},4},e::AbstractArray{Complex{T},4},mn::AbstractArray{T,5},mag::AbstractArray{T,3},Ninv::T)::AbstractArray{Complex{T},4} where T<:Real
    @avx for k ∈ axes(H,4), j ∈ axes(H,3), i ∈ axes(H,2), l in 0:0
        scale = mag[i,j,k] * Ninv
        H[1+l,i,j,k] =  (	e[1+l,i,j,k] * mn[1+l,2,i,j,k] + e[2+l,i,j,k] * mn[2+l,2,i,j,k] + e[3+l,i,j,k] * mn[3+l,2,i,j,k]	) * -scale  # -mag[i,j,k] * Ninv
		H[2+l,i,j,k] =  (	e[1+l,i,j,k] * mn[1+l,1,i,j,k] + e[2+l,i,j,k] * mn[2+l,1,i,j,k] + e[3+l,i,j,k] * mn[3+l,1,i,j,k]	) * scale   # mag[i,j,k] * Ninv
    end
    return H
end

# function eid!(e::AbstractArray{Complex{T},4},ε⁻¹,d::AbstractArray{Complex{T},4})::AbstractArray{Complex{T},4} where T<:Real
#     @avx for k ∈ axes(e,4), j ∈ axes(e,3), i ∈ axes(e,2), l in 0:0, h in 0:0
#         e[1+h,i,j,k] =  ε⁻¹[1+h,1+l,i,j,k]*d[1+l,i,j,k] + ε⁻¹[2+h,1+l,i,j,k]*d[2+l,i,j,k] + ε⁻¹[3+h,1+l,i,j,k]*d[3+l,i,j,k]
#         e[2+h,i,j,k] =  ε⁻¹[1+h,2+l,i,j,k]*d[1+l,i,j,k] + ε⁻¹[2+h,2+l,i,j,k]*d[2+l,i,j,k] + ε⁻¹[3+h,2+l,i,j,k]*d[3+l,i,j,k]
#         e[3+h,i,j,k] =  ε⁻¹[1+h,3+l,i,j,k]*d[1+l,i,j,k] + ε⁻¹[2+h,3+l,i,j,k]*d[2+l,i,j,k] + ε⁻¹[3+h,3+l,i,j,k]*d[3+l,i,j,k]
#     end
#     return e
# end

function eid!(e::AbstractArray{Complex{T},4},ε⁻¹,d::AbstractArray{Complex{T},4})::AbstractArray{Complex{T},4} where T<:Real
    @avx for k ∈ axes(e,4), j ∈ axes(e,3), i ∈ axes(e,2), l in 0:0, h in 0:0
        e[1+h,i,j,k] =  ε⁻¹[1+h,1+l,i,j,k]*d[1+l,i,j,k] + ε⁻¹[2+h,1+l,i,j,k]*d[2+l,i,j,k] + ε⁻¹[3+h,1+l,i,j,k]*d[3+l,i,j,k]
        e[2+h,i,j,k] =  ε⁻¹[1+h,2+l,i,j,k]*d[1+l,i,j,k] + ε⁻¹[2+h,2+l,i,j,k]*d[2+l,i,j,k] + ε⁻¹[3+h,2+l,i,j,k]*d[3+l,i,j,k]
        e[3+h,i,j,k] =  ε⁻¹[1+h,3+l,i,j,k]*d[1+l,i,j,k] + ε⁻¹[2+h,3+l,i,j,k]*d[2+l,i,j,k] + ε⁻¹[3+h,3+l,i,j,k]*d[3+l,i,j,k]
    end
    return e
end

function eid!(e::AbstractArray{Complex{T},4},ε⁻¹::AbstractArray{TA,3},d::AbstractArray{Complex{T},4})::AbstractArray{Complex{T},4} where {T<:Real,TA<:SMatrix{3,3}}
    er = reinterpret(reshape,SVector{3,Complex{T}},e)
	dr = reinterpret(reshape,SVector{3,Complex{T}},d)
	map!(*,er,ε⁻¹,dr)
	# map!(*,er,ε⁻¹,er)
    return e
end

function kxinv_tc!(e::AbstractArray{Complex{T},4},H::AbstractArray{Complex{T},4},mn::AbstractArray{T,5},inv_mag::AbstractArray{T,3})::AbstractArray{Complex{T},4} where T<:Real
    @avx for k ∈ axes(e,4), j ∈ axes(e,3), i ∈ axes(e,2), l in 0:0
		e[1+l,i,j,k] = ( H[1,i,j,k] * mn[1+l,2,i,j,k] - H[2,i,j,k] * mn[1+l,1,i,j,k] ) * inv_mag[i,j,k]
        e[2+l,i,j,k] = ( H[1,i,j,k] * mn[2+l,2,i,j,k] - H[2,i,j,k] * mn[2+l,1,i,j,k] ) * inv_mag[i,j,k]
        e[3+l,i,j,k] = ( H[1,i,j,k] * mn[3+l,2,i,j,k] - H[2,i,j,k] * mn[3+l,1,i,j,k] ) * inv_mag[i,j,k]
    end
    return e
end

function kxinv_ct!(H::AbstractArray{Complex{T},4},d::AbstractArray{Complex{T},4},mn::AbstractArray{T,5},inv_mag::AbstractArray{T,3},N::T)::AbstractArray{Complex{T},4} where T<:Real
    @avx for k ∈ axes(H,4), j ∈ axes(H,3), i ∈ axes(H,2), l in 0:0
        scale = inv_mag[i,j,k] * N
        H[1+l,i,j,k] =  (	d[1+l,i,j,k] * mn[1+l,2,i,j,k] + d[2+l,i,j,k] * mn[2+l,2,i,j,k] + d[3+l,i,j,k] * mn[3+l,2,i,j,k]	) * scale # inv_mag[i,j,k] * N
		H[2+l,i,j,k] =  (	d[1+l,i,j,k] * mn[1+l,1,i,j,k] + d[2+l,i,j,k] * mn[2+l,1,i,j,k] + d[3+l,i,j,k] * mn[3+l,1,i,j,k]	) * -scale # inv_mag[i,j,k] * N
    end
    return H
end

function ed_approx!(d::AbstractArray{Complex{T},4},ε_ave::AbstractArray{T,3},e::AbstractArray{Complex{T},4})::AbstractArray{Complex{T},4} where T<:Real
    @avx for k ∈ axes(e,4), j ∈ axes(e,3), i ∈ axes(e,2), l in 0:0
        d[1+l,i,j,k] =  ε_ave[i,j,k]*e[1+l,i,j,k]
        d[2+l,i,j,k] =  ε_ave[i,j,k]*e[2+l,i,j,k]
        d[3+l,i,j,k] =  ε_ave[i,j,k]*e[3+l,i,j,k]
    end
    return d
end

# 2D _M! and _P! subroutines

function kx_tc!(d::AbstractArray{Complex{T},3},H::AbstractArray{Complex{T},3},mn::AbstractArray{T,4},mag::AbstractArray{T,2})::AbstractArray{Complex{T},3} where T<:Real
    # @assert size(Y) === size(X)
    # @assert size(d,4) == 3
    # @assert size(H,4) === 2
    @avx for j ∈ axes(d,3), i ∈ axes(d,2), l in 0:0
	# @avx for i ∈ axes(d,1), j ∈ axes(d,2), l in 0:0
		# scale = -mag[i,j,k]
		d[1+l,i,j] = ( H[1,i,j] * mn[1+l,2,i,j] - H[2,i,j] * mn[1+l,1,i,j] ) * -mag[i,j]
        d[2+l,i,j] = ( H[1,i,j] * mn[2+l,2,i,j] - H[2,i,j] * mn[2+l,1,i,j] ) * -mag[i,j]
        d[3+l,i,j] = ( H[1,i,j] * mn[3+l,2,i,j] - H[2,i,j] * mn[3+l,1,i,j] ) * -mag[i,j]
    end
    return d
end

function zx_tc!(d::AbstractArray{Complex{T},3},H::AbstractArray{Complex{T},3},mn::AbstractArray{T,4})::AbstractArray{Complex{T},3} where T<:Real
    @avx for j ∈ axes(d,3), i ∈ axes(d,2), l in 0:0
		d[1+l,i,j] = -H[1,i,j] * mn[2+l,1,i,j] - H[2,i,j] * mn[2+l,2,i,j]
        d[2+l,i,j] =  H[1,i,j] * mn[1+l,1,i,j] + H[2,i,j] * mn[1+l,2,i,j]
    end
    return d
end

function kx_ct!(H::AbstractArray{Complex{T},3},e::AbstractArray{Complex{T},3},mn::AbstractArray{T,4},mag::AbstractArray{T,2},Ninv::T)::AbstractArray{Complex{T},3} where T<:Real
    @avx for j ∈ axes(H,3), i ∈ axes(H,2), l in 0:0
        scale = mag[i,j] * Ninv
        H[1+l,i,j] =  (	e[1+l,i,j] * mn[1+l,2,i,j] + e[2+l,i,j] * mn[2+l,2,i,j] + e[3+l,i,j] * mn[3+l,2,i,j]	) * -scale  # -mag[i,j] * Ninv
		H[2+l,i,j] =  (	e[1+l,i,j] * mn[1+l,1,i,j] + e[2+l,i,j] * mn[2+l,1,i,j] + e[3+l,i,j] * mn[3+l,1,i,j]	) * scale   # mag[i,j] * Ninv
    end
    return H
end

function eid!(e::AbstractArray{Complex{T},3},ε⁻¹,d::AbstractArray{Complex{T},3}) where T<:Real
    @avx for j ∈ axes(e,3), i ∈ axes(e,2), l in 0:0, h in 0:0
        e[1+h,i,j] =  ε⁻¹[1+h,1+l,i,j]*d[1+l,i,j] + ε⁻¹[2+h,1+l,i,j]*d[2+l,i,j] + ε⁻¹[3+h,1+l,i,j]*d[3+l,i,j]
        e[2+h,i,j] =  ε⁻¹[1+h,2+l,i,j]*d[1+l,i,j] + ε⁻¹[2+h,2+l,i,j]*d[2+l,i,j] + ε⁻¹[3+h,2+l,i,j]*d[3+l,i,j]
        e[3+h,i,j] =  ε⁻¹[1+h,3+l,i,j]*d[1+l,i,j] + ε⁻¹[2+h,3+l,i,j]*d[2+l,i,j] + ε⁻¹[3+h,3+l,i,j]*d[3+l,i,j]
    end
    return e
end

function eid!(e::AbstractArray{Complex{T},3},ε⁻¹::AbstractArray{TA,2},d::AbstractArray{Complex{T},3})::AbstractArray{Complex{T},3} where {T<:Real,TA<:SMatrix{3,3}}
    er = reinterpret(reshape,SVector{3,Complex{T}},e)
	dr = reinterpret(reshape,SVector{3,Complex{T}},d)
	map!(*,er,ε⁻¹,dr)
	# map!(*,er,ε⁻¹,er)
    return e
end

function kxinv_tc!(e::AbstractArray{Complex{T},3},H::AbstractArray{Complex{T},3},mn::AbstractArray{T,4},inv_mag::AbstractArray{T,2})::AbstractArray{Complex{T},3} where T<:Real
    @avx for j ∈ axes(e,3), i ∈ axes(e,2), l in 0:0
		e[1+l,i,j] = ( H[1,i,j] * mn[1+l,2,i,j] - H[2,i,j] * mn[1+l,1,i,j] ) * inv_mag[i,j]
        e[2+l,i,j] = ( H[1,i,j] * mn[2+l,2,i,j] - H[2,i,j] * mn[2+l,1,i,j] ) * inv_mag[i,j]
        e[3+l,i,j] = ( H[1,i,j] * mn[3+l,2,i,j] - H[2,i,j] * mn[3+l,1,i,j] ) * inv_mag[i,j]
    end
    return e
end

function kxinv_ct!(H::AbstractArray{Complex{T},3},d::AbstractArray{Complex{T},3},mn::AbstractArray{T,4},inv_mag::AbstractArray{T,2},N::T)::AbstractArray{Complex{T},3} where T<:Real
    @avx for j ∈ axes(H,3), i ∈ axes(H,2), l in 0:0
        scale = inv_mag[i,j] * N
        H[1+l,i,j] =  (	d[1+l,i,j] * mn[1+l,2,i,j] + d[2+l,i,j] * mn[2+l,2,i,j] + d[3+l,i,j] * mn[3+l,2,i,j]	) * scale # inv_mag[i,j] * N
		H[2+l,i,j] =  (	d[1+l,i,j] * mn[1+l,1,i,j] + d[2+l,i,j] * mn[2+l,1,i,j] + d[3+l,i,j] * mn[3+l,1,i,j]	) * -scale # inv_mag[i,j] * N
    end
    return H
end

function ed_approx!(d::AbstractArray{Complex{T},3},ε_ave::AbstractArray{T,2},e::AbstractArray{Complex{T},3})::AbstractArray{Complex{T},3} where T<:Real
    @avx for j ∈ axes(e,3), i ∈ axes(e,2), l in 0:0
        d[1+l,i,j] =  ε_ave[i,j]*e[1+l,i,j]
        d[2+l,i,j] =  ε_ave[i,j]*e[2+l,i,j]
        d[3+l,i,j] =  ε_ave[i,j]*e[3+l,i,j]
    end
    return d
end

# _M! and _P!

function _P!(Hout::AbstractArray{Complex{T},N}, Hin::AbstractArray{Complex{T},N},
	e::AbstractArray{Complex{T},N}, d::AbstractArray{Complex{T},N}, ε_ave::AbstractArray{T},
	mn::AbstractArray{T}, inv_mag::AbstractArray{T},
	𝓕!::FFTW.cFFTWPlan, 𝓕⁻¹!::FFTW.cFFTWPlan,
	Ninv::T)::AbstractArray{Complex{T},N} where {T<:Real,N}
	kxinv_tc!(e,Hin,mn,inv_mag);
	mul!(e,𝓕⁻¹!,e);
    ed_approx!(d,ε_ave,e);
    mul!(d,𝓕!,d);
    kxinv_ct!(Hout,d,mn,inv_mag,Ninv)
end

function _M!(Hout::AbstractArray{Complex{T},N}, Hin::AbstractArray{Complex{T},N},
	e::AbstractArray{Complex{T},N}, d::AbstractArray{Complex{T},N}, ε⁻¹,
	mn::AbstractArray{T}, mag::AbstractArray{T},
	𝓕!::FFTW.cFFTWPlan, 𝓕⁻¹!::FFTW.cFFTWPlan,
	Ninv::T)::AbstractArray{Complex{T},N} where {T<:Real,N}
    kx_tc!(d,Hin,mn,mag);
    mul!(d,𝓕!,d);
    eid!(e,ε⁻¹,d);
    mul!(e,𝓕⁻¹!,e);
    kx_ct!(Hout,e,mn,mag,Ninv)
end

"""
################################################################################
#																			   #
#			  Utility Function Definitions Needed for Constructors 			   #
#																			   #
################################################################################
"""

function mag_m_n!(mag,m,n,k⃗::SVector{3,T},g⃗) where T <: Real
	# for iz ∈ axes(g⃗,3), iy ∈ axes(g⃗,2), ix ∈ axes(g⃗,1) #, l in 0:0
	local ẑ = SVector(0.,0.,1.)
	local ŷ = SVector(0.,1.,0.)
	@fastmath @inbounds for i ∈ eachindex(g⃗)
		@inbounds kpg::SVector{3,T} = k⃗ - g⃗[i]
		@inbounds mag[i] = norm(kpg)
		@inbounds n[i] =  ( ( abs2(kpg[1]) + abs2(kpg[2]) ) > 0. ) ?  normalize( cross( ẑ, kpg ) ) : SVector(-1.,0.,0.) #ŷ
		@inbounds m[i] =  normalize( cross( n[i], kpg )  )
	end
	return mag,m,n
end

mag_m_n!(mag,m,n,kz::T,g⃗) where T <: Real = mag_m_n!(mag,m,n,SVector{3,T}(0.,0.,kz),g⃗)

function mag_m_n2(k⃗::SVector{3,T},g⃗::AbstractArray) where T <: Real
	# g⃗ₜ_zero_mask = Zygote.@ignore(  sum(abs2,g⃗[1:2,:,:,:];dims=1)[1,:,:,:] .> 0. );
	g⃗ₜ_zero_mask = Zygote.@ignore(  sum(abs2,g⃗[1:2,:,:,:];dims=1)[1,:,:,:] .> 0. );
	g⃗ₜ_zero_mask! = Zygote.@ignore( .!(g⃗ₜ_zero_mask) );
	local ŷ = [0.; 1. ;0.]
	local zxinds = [2; 1; 3]
	local zxscales = [-1; 1. ;0.]
	local xinds1 = [2; 3; 1]
	local xinds2 = [3; 1; 2]
	@tullio kpg[ix,iy,iz] := k⃗[a] - g⃗[a,ix,iy,iz] fastmath=false
	@tullio mag[ix,iy,iz] := sqrt <| kpg[a,ix,iy,iz]^2 fastmath=false
	@tullio nt[ix,iy,iz,a] := zxscales[a] * kpg[zxinds[a],ix,iy,iz] * g⃗ₜ_zero_mask[ix,iy,iz] + ŷ[a] * g⃗ₜ_zero_mask![ix,iy,iz]  nograd=(zxscales,zxinds,ŷ,g⃗ₜ_zero_mask,g⃗ₜ_zero_mask!) fastmath=false
	@tullio nmag[ix,iy,iz] := sqrt <| nt[a,ix,iy,iz]^2 fastmath=false
	@tullio n[a,ix,iy,iz] := nt[a,ix,iy,iz] / nmag[ix,iy,iz] fastmath=false
	@tullio mt[a,ix,iy,iz] := n[xinds1[a],ix,iy,iz] * kpg[xinds2[a],ix,iy,iz] - kpg[xinds1[a],ix,iy,iz] * n[xinds2[a],ix,iy,iz] nograd=(xinds1,xinds2) fastmath=false
	@tullio mmag[ix,iy,iz] := sqrt <| mt[a,ix,iy,iz]^2 fastmath=false
	@tullio m[a,ix,iy,iz] := mt[a,ix,iy,iz] / mmag[ix,iy,iz] fastmath=false
	return mag, m, n
end

function mag_m_n2(kz::T,g⃗::AbstractArray) where T <: Real
	mag_m_n2(SVector{3,T}(0.,0.,kz),g⃗)
end

function mag_m_n(k⃗::SVector{3,T},g⃗::AbstractArray{SVector{3,T2}}) where {T<:Real,T2<:Real}
	# for iz ∈ axes(g⃗,3), iy ∈ axes(g⃗,2), ix ∈ axes(g⃗,1) #, l in 0:0
	local ẑ = SVector(0.,0.,1.)
	local ŷ = SVector(0.,1.,0.)
	n = Buffer(g⃗,size(g⃗))
	m = Buffer(g⃗,size(g⃗))
	mag = Buffer(zeros(T,size(g⃗)),size(g⃗))
	# n = bufferfrom(zeros(SVector{3,T},size(g⃗)))
	# m = bufferfrom(zeros(SVector{3,T},size(g⃗)))
	# mag = bufferfrom(zeros(T,size(g⃗)))
	@fastmath @inbounds for i ∈ eachindex(g⃗)
		@inbounds kpg::SVector{3,T} = k⃗ - g⃗[i]
		@inbounds mag[i] = norm(kpg)
		@inbounds n[i] =   ( ( abs2(kpg[1]) + abs2(kpg[2]) ) > 0. ) ?  normalize( cross( ẑ, kpg ) ) : SVector(-1.,0.,0.) #ŷ
		@inbounds m[i] =  normalize( cross( n[i], kpg )  )
	end
	return copy(mag), copy(m), copy(n) # HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},T}(reinterpret(reshape,Float64,copy(m))), HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},T}(reinterpret(reshape,Float64,copy(n)))
end

function mag_m_n(kz::T,g⃗::AbstractArray{SVector{3,T2}}) where {T<:Real,T2<:Real}
	mag_m_n(SVector{3,T}(0.,0.,kz),g⃗)
end

mag_m_n(k::Real,grid::Grid) = mag_m_n(k, g⃗(grid))

"""
(māg,m̄,n̄) → k̄ map
"""
function ∇ₖmag_m_n(māg,m̄,n̄,mag,m⃗,n⃗;dk̂=SVector(0.,0.,1.))
	kp̂g_over_mag = cross.(m⃗,n⃗)./mag
	k̄_mag = sum( māg .* dot.( kp̂g_over_mag, (dk̂,) ) .* mag )
	k̄_m = -sum( dot.( m̄ , cross.(m⃗, cross.( kp̂g_over_mag, (dk̂,) ) ) ) )
	k̄_n = -sum( dot.( n̄ , cross.(n⃗, cross.( kp̂g_over_mag, (dk̂,) ) ) ) )
	return +( k̄_mag, k̄_m, k̄_n )
end

function ∇ₖmag_m_n(ΔΩ,Ω;dk̂=SVector(0.,0.,1.))
	māg,m̄,n̄ = ΔΩ
	mag,m⃗,n⃗ = Ω
	kp̂g_over_mag = cross.(m⃗,n⃗)./mag
	k̄_mag = sum( māg .* dot.( kp̂g_over_mag, (dk̂,) ) .* mag )
	k̄_m = -sum( dot.( m̄ , cross.(m⃗, cross.( kp̂g_over_mag, (dk̂,) ) ) ) )
	k̄_n = -sum( dot.( n̄ , cross.(n⃗, cross.( kp̂g_over_mag, (dk̂,) ) ) ) )
	return +( k̄_mag, k̄_m, k̄_n )
end

function rrule(::typeof(mag_m_n),k⃗::SVector{3,T},g⃗::AbstractArray{SVector{3,T}}) where T <: Real
	local ẑ = SVector(0.,0.,1.)
	local ŷ = SVector(0.,1.,0.)
	n_buf = Buffer(g⃗,size(g⃗))
	m_buf = Buffer(g⃗,size(g⃗))
	kpg_buf = Buffer(g⃗,size(g⃗))
	mag_buf = Buffer(zeros(T,size(g⃗)),size(g⃗))
	@fastmath @inbounds for i ∈ eachindex(g⃗)
		@inbounds kpg_buf[i] = k⃗ - g⃗[i]
		@inbounds mag_buf[i] = norm(kpg_buf[i])
		@inbounds n_buf[i] =   ( ( abs2(kpg_buf[i][1]) + abs2(kpg_buf[i][2]) ) > 0. ) ?  normalize( cross( ẑ, kpg_buf[i] ) ) : SVector(-1.,0.,0.) # ŷ
		@inbounds m_buf[i] =  normalize( cross( n_buf[i], kpg_buf[i] )  )
	end
	mag_m⃗_n⃗ = (copy(mag_buf), copy(m_buf), copy(n_buf))
	kp⃗g = copy(kpg_buf)
	mag_m_n_pullback(ΔΩ) = let Ω=mag_m⃗_n⃗, kp⃗g=kp⃗g, dk̂=normalize(k⃗)
		māg,m̄,n̄ = ΔΩ
		mag,m⃗,n⃗ = Ω
		ê_over_mag = cross.( kp⃗g, (dk̂,) ) ./ mag.^2
		k̄ = sum( māg .* dot.( kp⃗g, (dk̂,) ) ./ mag )
		k̄ -= sum( dot.( m̄ , cross.(m⃗, ê_over_mag ) ) )
		k̄ -= sum( dot.( n̄ , cross.(n⃗, ê_over_mag ) ) )
		return ( NoTangent(), k̄*dk̂, ZeroTangent() )
	end
    return (mag_m⃗_n⃗ , mag_m_n_pullback)
end

# function mag_mn(k::T1,g::AbstractArray{SVector{3,T2},3}) where {T1<:Real,T2<:Real}
# 	mag, m⃗, n⃗ = mag_m_n(k,g)
# 	HybridArray{Tuple{3,2,Dynamic(),Dynamic(),Dynamic()},T1}(copy(reshape(reinterpret(Val(T1),hcat.(m⃗,n⃗)),(3,2,size(g)...))))
# end
#
# function mag_mn(k::T1,g::Grid{3,T2}) where {T1<:Real,T2<:Real}
# 	mag, m⃗, n⃗ = mag_m_n(k,g)
# 	HybridArray{Tuple{3,2,Dynamic(),Dynamic(),Dynamic()},T1}(copy(reshape(reinterpret(Val(T1),hcat.(m⃗,n⃗)),(3,2,size(g)...))))
# end

function mag_mn(k::SVector{3,T1},g::Grid{3,T2}) where {T1<:Real,T2<:Real}
	mag, m⃗, n⃗ = mag_m_n(k,g⃗(g))
	mn = copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...)))
	# return mag, HybridArray{Tuple{3,2,Dynamic(),Dynamic(),Dynamic()},T1}(copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...))))
	return mag, mn
end

function mag_mn(k::SVector{3,T1},g::Grid{2,T2}) where {T1<:Real,T2<:Real}
	mag, m⃗, n⃗ = mag_m_n(k,g⃗(g))
	mn = copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...)))
	# return mag, HybridArray{Tuple{3,2,Dynamic(),Dynamic()},T1}(copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...))))
	return mag, mn
end

function mag_mn(k::SVector{3,T1},g::AbstractArray{SVector{3,T2},3}) where {T1<:Real,T2<:Real}
	mag, m⃗, n⃗ = mag_m_n(k,g)
	mn = copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...)))
	# return mag, HybridArray{Tuple{3,2,Dynamic(),Dynamic(),Dynamic()},T1}(copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...))))
	return mag, mn
end

function mag_mn(k::SVector{3,T1},g::AbstractArray{SVector{3,T2},2}) where {T1<:Real,T2<:Real}
	mag, m⃗, n⃗ = mag_m_n(k,g)
	mn = copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...)))
	# return mag, HybridArray{Tuple{3,2,Dynamic(),Dynamic()},T1}(copy(reshape(reinterpret(T1,hcat.(m⃗,n⃗)),(3,2,size(g)...))))
	return mag, mn
end

function mag_mn(k::T,g::TG) where {T<:Real,TG}
	return mag_mn(SVector{3,T}(0.,0.,k),g)
end

function mag_mn!(mag,mn::AbstractArray{T1,NDp2},k⃗::SVector{3,T2},g⃗) where {T1<:Real,T2<:Real,NDp2}
	local ẑ = SVector(0.,0.,1.)
	local ŷ = SVector(0.,1.,0.)
	# mv = view(mn,1:3,1,eachindex(g⃗)...)
	# nv = view(mn,1:3,2,eachindex(g⃗)...)
	# mvs = reinterpret(reshape,SVector{3,T1},mv)
	# nvs = reinterpret(reshape,SVector{3,T1},nv)
	kpg = zero(k⃗)
	@fastmath @inbounds for i ∈ eachindex(g⃗)
		@inbounds kpg = k⃗ - g⃗[i]
		@inbounds mag[i] = norm(kpg)
		@inbounds mn[1:3,2,i] .=  ( ( abs2(kpg[1]) + abs2(kpg[2]) ) > 0. ) ?  normalize( cross( ẑ, kpg ) ) : SVector(-1.,0.,0.) #[-1.,0.,0.] #ŷ
		@inbounds mn[1:3,1,i] .=  normalize( cross( mn[1:3,2,i], kpg )  )
	end
	# return mag,m,n
	return mag, mn
end

mag_mn!(mag,mn,kz::T,g⃗) where T <: Real = mag_mn!(mag,mn,SVector{3,T}(0.,0.,kz),g⃗)

"""
(māg,m̄n̄) → k̄ map
"""
function ∇ₖmag_mn(māg::AbstractArray{T1,2},mn̄,mag::AbstractArray{T2,2},mn) where {T1<:Real,T2<:Number}
	m = view(mn,:,1,:,:)
	n = view(mn,:,2,:,:)
	@tullio kp̂g_over_mag[i,ix,iy] := m[mod(i-2),ix,iy] * n[mod(i-1),ix,iy] / mag[ix,iy] - m[mod(i-1),ix,iy] * n[mod(i-2),ix,iy] / mag[ix,iy] (i in 1:3)
	kp̂g_over_mag_x_dk̂ = _cross(kp̂g_over_mag,dk̂)
	@tullio k̄_mag := māg[ix,iy] * mag[ix,iy] * kp̂g_over_mag[j,ix,iy] * dk̂[j]
	@tullio k̄_mn := -conj(mn̄)[imn,i,ix,iy] * mn[imn,mod(i-2),ix,iy] * kp̂g_over_mag_x_dk̂[mod(i-1),ix,iy] + conj(mn̄)[imn,i,ix,iy] * mn[imn,mod(i-1),ix,iy] * kp̂g_over_mag_x_dk̂[mod(i-2),ix,iy] (i in 1:3)
	k̄_magmn = k̄_mag + k̄_mn
	return k̄_magmn
end

function ∇ₖmag_mn(māg::AbstractArray{T1,3},mn̄,mag::AbstractArray{T2,3},mn) where {T1<:Real,T2<:Number}
	m = view(mn,:,1,:,:,:)
	n = view(mn,:,2,:,:,:)
	@tullio kp̂g_over_mag[i,ix,iy,iz] := m[mod(i-2),ix,iy,iz] * n[mod(i-1),ix,iy,iz] / mag[ix,iy,iz] - m[mod(i-1),ix,iy,iz] * n[mod(i-2),ix,iy,iz] / mag[ix,iy,iz] (i in 1:3)
	kp̂g_over_mag_x_dk̂ = _cross(kp̂g_over_mag,dk̂)
	@tullio k̄_mag := māg[ix,iy,iz] * mag[ix,iy,iz] * kp̂g_over_mag[j,ix,iy,iz] * dk̂[j]
	@tullio k̄_mn := -conj(mn̄)[imn,i,ix,iy,iz] * mn[imn,mod(i-2),ix,iy,iz] * kp̂g_over_mag_x_dk̂[mod(i-1),ix,iy,iz] + conj(mn̄)[imn,i,ix,iy,iz] * mn[imn,mod(i-1),ix,iy,iz] * kp̂g_over_mag_x_dk̂[mod(i-2),ix,iy,iz] (i in 1:3)
	k̄_magmn = k̄_mag + k̄_mn
	return k̄_magmn
end

# function rrule(::typeof(mag_mn),k⃗::SVector{3,T},g⃗::AbstractArray{SVector{3,T}}) where T <: Real
# 	local ẑ = SVector(0.,0.,1.)
# 	local ŷ = SVector(0.,1.,0.)
# 	n_buf = Buffer(g⃗,size(g⃗))
# 	m_buf = Buffer(g⃗,size(g⃗))
# 	kpg_buf = Buffer(g⃗,size(g⃗))
# 	mag_buf = Buffer(zeros(T,size(g⃗)),size(g⃗))
# 	@fastmath @inbounds for i ∈ eachindex(g⃗)
# 		@inbounds kpg_buf[i] = k⃗ - g⃗[i]
# 		@inbounds mag_buf[i] = norm(kpg_buf[i])
# 		@inbounds n_buf[i] =   ( ( abs2(kpg_buf[i][1]) + abs2(kpg_buf[i][2]) ) > 0. ) ?  normalize( cross( ẑ, kpg_buf[i] ) ) : SVector(-1.,0.,0.) # ŷ
# 		@inbounds m_buf[i] =  normalize( cross( n_buf[i], kpg_buf[i] )  )
# 	end
# 	mag_m⃗_n⃗ = (copy(mag_buf), copy(m_buf), copy(n_buf))
# 	kp⃗g = copy(kpg_buf)
# 	mag_m_n_pullback(ΔΩ) = let Ω=mag_m⃗_n⃗, kp⃗g=kp⃗g, dk̂=normalize(k⃗)
# 		māg,m̄,n̄ = ΔΩ
# 		mag,m⃗,n⃗ = Ω
# 		ê_over_mag = cross.( kp⃗g, (dk̂,) ) ./ mag.^2
# 		k̄ = sum( māg .* dot.( kp⃗g, (dk̂,) ) ./ mag )
# 		k̄ -= sum( dot.( m̄ , cross.(m⃗, ê_over_mag ) ) )
# 		k̄ -= sum( dot.( n̄ , cross.(n⃗, ê_over_mag ) ) )
# 		return ( NoTangent(), k̄*dk̂, ZeroTangent() )
# 	end
#     return (mag_m⃗_n⃗ , mag_m_n_pullback)
# end

"""
################################################################################
#																			   #
#							  Struct Definitions 							   #
#																			   #
################################################################################
"""

mutable struct HelmholtzMap{ND,T,NDp1,NDp2} <: LinearMap{T}
    k⃗::SVector{3,T}
	Nx::Int
	Ny::Int
	Nz::Int
	N::Int
	Ninv::T
	g⃗::Array{SVector{3,T},ND}
	mag::Array{T,ND}
    # m⃗::Array{SVector{3,T},ND}
	# n⃗::Array{SVector{3,T},ND}
	# m::Base.ReinterpretArray{T, NDp1, SVector{3, T}, Array{SVector{3, T}, ND}, true}
	# n::Base.ReinterpretArray{T, NDp1, SVector{3, T}, Array{SVector{3, T}, ND}, true}
	mn::Array{T,NDp2}
    e::Array{Complex{T},NDp1}
    d::Array{Complex{T},NDp1}
    𝓕!::FFTW.cFFTWPlan{Complex{T}, -1, true, NDp1, UnitRange{Int64}}
	𝓕⁻¹!::FFTW.cFFTWPlan{Complex{T}, 1, true, NDp1, UnitRange{Int64}}
	𝓕::FFTW.cFFTWPlan{Complex{T}, -1, false, NDp1, UnitRange{Int64}}
	𝓕⁻¹::FFTW.cFFTWPlan{Complex{T}, 1, false, NDp1, UnitRange{Int64}}
	ε⁻¹::Array{T,NDp2}
	ε_ave::Array{T,ND}  # for preconditioner
	inv_mag::Array{T,ND} # for preconditioner
	shift::T
end

mutable struct HelmholtzPreconditioner{ND,T,NDp1,NDp2} <: LinearMap{T}
	M̂::HelmholtzMap{ND,T,NDp1,NDp2}
end

mutable struct ModeSolver{ND,T,NDp1,NDp2}
	# geom::Geometry
	# materials #::Vector{<:AbstractMaterial}
	grid::Grid{ND,T}
	M̂::HelmholtzMap{ND,T,NDp1,NDp2}
	P̂::HelmholtzPreconditioner{ND,T,NDp1,NDp2}
	# eigs_itr::IterativeSolvers.LOBPCGIterator
	H⃗::Matrix{Complex{T}}
	ω²::Vector{Complex{T}}
	∂ω²∂k::Vector{T}
	# λ⃗::Vector{Complex{T}}
	# b⃗::Vector{Complex{T}}
	# λd::HybridArray
	# λẽ::HybridArray
	# ε⁻¹_bar::Array{SMatrix{3,3,T,9}, ND}
	# kx̄_m⃗::Array{SVector{3, T}, ND}
	# kx̄_n⃗::Array{SVector{3, T}, ND}
	# māg::Array{T,ND}
	# k̄_kx::SVector{3,T}
	# ω̄::T
	# adj_itr::IterativeSolvers.BiCGStabIterable
	# corner_sinds::Array{Int,ND}
	# sinds_proc::Array #{NTuple{8,Int},ND}
	# Srvol::Array{Tuple{SMatrix{3,3,T,9},T},ND}
	# minds::Vector{Int}
end

"""
################################################################################
#																			   #
#							  Constructor Methods 							   #
#																			   #
################################################################################
"""

function HelmholtzMap(k⃗::AbstractVector{T}, ε⁻¹, gr::Grid{3,T}; shift=0. ) where {ND,T<:Real}
	g⃗s = g⃗(gr)
	# mag, m⃗, n⃗ = mag_m_n(k⃗,g⃗s)
	mag, mn = mag_mn(k⃗,g⃗s)
	d0 = randn(Complex{T}, (3,size(gr)...))
	fftax = _fftaxes(gr)
	return HelmholtzMap{3,T,4,5}(
			SVector{3,T}(k⃗),
			gr.Nx,
			gr.Ny,
			gr.Nz,
			N(gr),
			1. / N(gr),
			g⃗s,
			mag,
			mn,
			copy(d0),
			copy(d0),
			# m⃗,
			# n⃗,
			# HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},T}(reinterpret(reshape,T,m⃗)),
			# HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},T}(reinterpret(reshape,T,n⃗)),
		    # HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},Complex{T}}(d0),# (Array{T}(undef,(Nx,Ny,Nz,3))),
		    # HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},Complex{T}}(d0),# (Array{T}(undef,(Nx,Ny,Nz,3))),
			plan_fft!(d0,fftax,flags=FFTW.PATIENT), # planned in-place FFT operator 𝓕!
			plan_bfft!(d0,fftax,flags=FFTW.PATIENT), # planned in-place iFFT operator 𝓕⁻¹!
			plan_fft(d0,fftax,flags=FFTW.PATIENT), # planned in-place FFT operator 𝓕!
			plan_bfft(d0,fftax,flags=FFTW.PATIENT), # planned in-place iFFT operator 𝓕⁻¹!
			ε⁻¹,
			[ 3. * inv(sum(diag(ε⁻¹[:,:,I]))) for I in eachindex(gr)], #[ 3. * inv(sum(diag(einv))) for einv in ε⁻¹],
			[ inv(mm) for mm in mag ], # inverse |k⃗+g⃗| magnitudes for precond. ops
			shift,
		)
end

function HelmholtzMap(k⃗::AbstractVector{T}, ε⁻¹, gr::Grid{2,T}; shift=0. ) where {ND,T<:Real}
	g⃗s = g⃗(gr)
	# mag, m⃗, n⃗ = mag_m_n(k⃗,g⃗s)
	mag, mn = mag_mn(k⃗,g⃗s)
	d0 = randn(Complex{T}, (3,size(gr)...))
	fftax = _fftaxes(gr)
	return HelmholtzMap{2,T,3,4}(
			SVector{3,T}(k⃗),
			gr.Nx,
			gr.Ny,
			gr.Nz,
			N(gr),
			1. / N(gr),
			g⃗s,
			mag,
			mn,
			copy(d0),
			copy(d0),
			# m⃗,
			# n⃗,
			# HybridArray{Tuple{3,Dynamic(),Dynamic()},T}(reinterpret(reshape,T,m⃗)),
			# HybridArray{Tuple{3,Dynamic(),Dynamic()},T}(reinterpret(reshape,T,n⃗)),
		    # HybridArray{Tuple{3,Dynamic(),Dynamic()},Complex{T}}(d0),# (Array{T}(undef,(Nx,Ny,Nz,3))),
		    # HybridArray{Tuple{3,Dynamic(),Dynamic()},Complex{T}}(d0),# (Array{T}(undef,(Nx,Ny,Nz,3))),
			plan_fft!(d0,fftax,flags=FFTW.PATIENT), # planned in-place FFT operator 𝓕!
			plan_bfft!(d0,fftax,flags=FFTW.PATIENT), # planned in-place iFFT operator 𝓕⁻¹!
			plan_fft(d0,fftax,flags=FFTW.PATIENT), # planned in-place FFT operator 𝓕!
			plan_bfft(d0,fftax,flags=FFTW.PATIENT), # planned in-place iFFT operator 𝓕⁻¹!
			ε⁻¹,
			[ 3. * inv(sum(diag(ε⁻¹[:,:,I]))) for I in eachindex(gr)], # [ 3. * inv(sum(diag(einv))) for einv in ε⁻¹],
			[ inv(mm) for mm in mag ], # inverse |k⃗+g⃗| magnitudes for precond. ops
			shift,
		)
end

function HelmholtzMap(kz::T, ε⁻¹, gr::Grid; shift=0.) where {T<:Real}
	HelmholtzMap(SVector{3,T}(0.,0.,kz), ε⁻¹, gr::Grid; shift)
end

# function ModeSolver(k⃗::SVector{3,T}, geom::Geometry, grid::Grid{ND}; nev=1, tol=1e-8, maxiter=3000, ω₀=1/1.55, constraint=nothing,) where {ND,T<:Real}
function ModeSolver(k⃗::SVector{3,T}, ε⁻¹, grid::Grid{2}; nev=1, tol=1e-8, maxiter=3000, ω₀=1/1.55, constraint=nothing,) where {T<:Real}
	# run inital smoothing sub-processes
	# ε⁻¹ = εₛ⁻¹( (1. / ω₀), geom, grid)

	# sinds,sinds_proc,Srvol,mats,minds,ε⁻¹ = _εₛ⁻¹_init( (1. / ω₀), geom.shapes, grid)
	M̂ = HelmholtzMap(k⃗, ε⁻¹, grid)
	P̂ = HelmholtzPreconditioner(M̂)
	# eigs_itr = LOBPCGIterator(M̂,false,randn(eltype(M̂),(size(M̂)[1],nev)),P̂,constraint)
	# λ⃗ = randn(Complex{T},2*M̂.N)
	# b⃗ = similar(λ⃗)
	# adj_itr = bicgstabl_iterator!(λ⃗, M̂ - ( 1. * I ), b⃗, 2;		# last entry is `l`::Int = # of GMRES iterations
    #                          Pl = Identity(),
    #                          max_mv_products = size(M̂, 2),
    #                          abstol = zero(T),
    #                          reltol = sqrt(eps(T)),
    #                          initial_zero = false)
	ModeSolver{2,T,3,4}(
		# geom,
		# mats,
		grid,
		M̂,
		P̂,
		# eigs_itr,
		randn(Complex{T},2*N(grid),nev), #eigs_itr.XBlocks.block,
		zeros(Complex{T},nev),
		zeros(T,nev),
		# λ⃗,
		# b⃗,
		# similar(M̂.d),							# λ⃗d
		# similar(M̂.e),							# λ⃗ẽ
		# similar(M̂.ε⁻¹),						# ε⁻¹_bar
		# similar(M̂.m⃗),							 # kx̄_m⃗
		# similar(M̂.n⃗),							# kx̄_n⃗
		# similar(M̂.mag),						# māg
		# zero(SVector{3,Float64}),				# k̄_kx
		# 0.,										# ω̄
		# adj_itr,
		# sinds,
		# sinds_proc,
		# Srvol,
		# minds,
	)
end

# function ModeSolver(kz::T, ε⁻¹, grid::Grid{2}; nev=1, tol=1e-8, maxiter=3000,constraint=nothing,) where {T<:Real}
# 	ModeSolver(SVector{3,T}(0.,0.,kz), ε⁻¹, grid; nev, tol, maxiter, constraint)
# end
#
# function ModeSolver(ω, p::AbstractVector, geom_fn::F, grid::Grid{2}; nev=1, tol=1e-8, maxiter=3000,constraint=nothing,) where {F<:Function}
# 	ε⁻¹ = smooth(ω,p,:fεs,true,geom_fn,grid);
# 	return ModeSolver(k_guess(ω,ε⁻¹), ε⁻¹, grid; nev, tol, maxiter, constraint)
# end

function ModeSolver(k⃗::SVector{3,T}, ε⁻¹, grid::Grid{3}; nev=1, tol=1e-8, maxiter=3000, ω₀=1/1.55, constraint=nothing,) where {T<:Real}
	# run inital smoothing sub-processes
	# ε⁻¹ = εₛ⁻¹( (1. / ω₀), geom, grid)

	# sinds,sinds_proc,Srvol,mats,minds,ε⁻¹ = _εₛ⁻¹_init( (1. / ω₀), geom.shapes, grid)
	M̂ = HelmholtzMap(k⃗, ε⁻¹, grid)
	P̂ = HelmholtzPreconditioner(M̂)
	# eigs_itr = LOBPCGIterator(M̂,false,randn(eltype(M̂),(size(M̂)[1],nev)),P̂,constraint)
	# λ⃗ = randn(Complex{T},2*M̂.N)
	# b⃗ = similar(λ⃗)
	# adj_itr = bicgstabl_iterator!(λ⃗, M̂ - ( 1. * I ), b⃗, 2;		# last entry is `l`::Int = # of GMRES iterations
    #                          Pl = Identity(),
    #                          max_mv_products = size(M̂, 2),
    #                          abstol = zero(T),
    #                          reltol = sqrt(eps(T)),
    #                          initial_zero = false)
	ModeSolver{3,T,4,5}(
		# geom,
		# mats,
		grid,
		M̂,
		P̂,
		# eigs_itr,
		randn(Complex{T},2*N(grid),nev), #eigs_itr.XBlocks.block,
		zeros(Complex{T},nev),
		zeros(T,nev),
		# λ⃗,
		# b⃗,
		# similar(M̂.d),							# λ⃗d
		# similar(M̂.e),							# λ⃗ẽ
		# similar(M̂.ε⁻¹),						# ε⁻¹_bar
		# similar(M̂.m⃗),							 # kx̄_m⃗
		# similar(M̂.n⃗),							# kx̄_n⃗
		# similar(M̂.mag),						# māg
		# zero(SVector{3,Float64}),				# k̄_kx
		# 0.,										# ω̄
		# adj_itr,
		# sinds,
		# sinds_proc,
		# Srvol,
		# minds,
	)
end

function ModeSolver(kz::T, ε⁻¹, grid::Grid{ND}; nev=1, tol=1e-8, maxiter=3000,constraint=nothing,) where {ND,T<:Real}
	ModeSolver(SVector{3,T}(0.,0.,kz), ε⁻¹, grid; nev, tol, maxiter, constraint)
end

function ModeSolver(ω, p::AbstractVector, geom_fn::F, grid::Grid{ND}; nev=1, tol=1e-8, maxiter=3000,constraint=nothing,) where {ND,F<:Function}
	ε⁻¹ = smooth(ω,p,:fεs,true,geom_fn,grid);
	return ModeSolver(k_guess(ω,ε⁻¹), ε⁻¹, grid; nev, tol, maxiter, constraint)
end

"""
################################################################################
#																			   #
#							  	Struct Methods 								   #
#																			   #
################################################################################
"""

function (M̂::HelmholtzMap{2,T,3,4})(Hout::AbstractArray{Complex{T},3}, Hin::AbstractArray{Complex{T},3}) where T<:Real
	_M!(Hout,Hin,M̂.e,M̂.d,M̂.ε⁻¹,M̂.mn,M̂.mag,M̂.𝓕!,M̂.𝓕⁻¹!,M̂.Ninv)
end

function (M̂::HelmholtzMap{3,T,4,5})(Hout::AbstractArray{Complex{T},4}, Hin::AbstractArray{Complex{T},4}) where T<:Real
	_M!(Hout,Hin,M̂.e,M̂.d,M̂.ε⁻¹,M̂.mn,M̂.mag,M̂.𝓕!,M̂.𝓕⁻¹!,M̂.Ninv)
end

function (M̂::HelmholtzMap{2,T,3,4})(Hout::AbstractVector{Complex{T}}, Hin::AbstractVector{Complex{T}}) where T<:Real
	@inbounds Hin_arr = reshape(Hin,(2,M̂.Nx,M̂.Ny))
	@inbounds Hout_arr = reshape(Hout,(2,M̂.Nx,M̂.Ny))
	vec( _M!(Hout_arr,Hin_arr,M̂.e,M̂.d,M̂.ε⁻¹,M̂.mn,M̂.mag,M̂.𝓕!,M̂.𝓕⁻¹!,M̂.Ninv) )
end

function (M̂::HelmholtzMap{3,T,4,5})(Hout::AbstractVector{Complex{T}}, Hin::AbstractVector{Complex{T}}) where T<:Real
	@inbounds Hin_arr = reshape(Hin,(2,M̂.Nx,M̂.Ny,M̂.Nz))
	@inbounds Hout_arr = reshape(Hout,(2,M̂.Nx,M̂.Ny,M̂.Nz))
	vec( _M!(Hout_arr,Hin_arr,M̂.e,M̂.d,M̂.ε⁻¹,M̂.mn,M̂.mag,M̂.𝓕!,M̂.𝓕⁻¹!,M̂.Ninv) )
end

function (P̂::HelmholtzPreconditioner)(Hout::AbstractArray{T,3}, Hin::AbstractArray{T,3}) where T<:Union{Real, Complex}
	_P!(Hout,Hin,P̂.M̂.e,P̂.M̂.d,P̂.M̂.ε_ave,P̂.M̂.mn,P̂.M̂.inv_mag,P̂.M̂.𝓕!,P̂.M̂.𝓕⁻¹!,P̂.M̂.Ninv)
end

function (P̂::HelmholtzPreconditioner)(Hout::AbstractArray{T,4}, Hin::AbstractArray{T,4}) where T<:Union{Real, Complex}
	_P!(Hout,Hin,P̂.M̂.e,P̂.M̂.d,P̂.M̂.ε_ave,P̂.M̂.mn,P̂.M̂.inv_mag,P̂.M̂.𝓕!,P̂.M̂.𝓕⁻¹!,P̂.M̂.Ninv)
end

function (P̂::HelmholtzPreconditioner{2})(Hout::AbstractVector{T}, Hin::AbstractVector{T}) where T<:Union{Real, Complex}
	@inbounds Hin_arr = reshape(Hin,(2,P̂.M̂.Nx,P̂.M̂.Ny))
	@inbounds Hout_arr = reshape(Hout,(2,P̂.M̂.Nx,P̂.M̂.Ny))
	vec( _P!(Hout_arr,Hin_arr,P̂.M̂.e,P̂.M̂.d,P̂.M̂.ε_ave,P̂.M̂.mn,P̂.M̂.inv_mag,P̂.M̂.𝓕!,P̂.M̂.𝓕⁻¹!,P̂.M̂.Ninv) )
end

function (P̂::HelmholtzPreconditioner{3})(Hout::AbstractVector{T}, Hin::AbstractVector{T}) where T<:Union{Real, Complex}
	@inbounds Hin_arr = reshape(Hin,(2,P̂.M̂.Nx,P̂.M̂.Ny,P̂.M̂.Nz))
	@inbounds Hout_arr = reshape(Hout,(2,P̂.M̂.Nx,P̂.M̂.Ny,P̂.M̂.Nz))
	vec( _P!(Hout_arr,Hin_arr,P̂.M̂.e,P̂.M̂.d,P̂.M̂.ε_ave,P̂.M̂.mn,P̂.M̂.inv_mag,P̂.M̂.𝓕!,P̂.M̂.𝓕⁻¹!,P̂.M̂.Ninv) )
end

function Base.:(*)(M̂::HelmholtzMap, x::AbstractVector)
    #length(x) == A.N || throw(DimensionMismatch())
    y = similar(x, promote_type(eltype(M̂), eltype(x)), 2*M̂.N)
    M̂(y, x)
end

function _unsafe_mul!(y::AbstractVecOrMat, M̂::HelmholtzMap, x::AbstractVector)
    M̂(y, x)
end

function Base.:(*)(P̂::HelmholtzPreconditioner, x::AbstractVector)
    #length(x) == A.N || throw(DimensionMismatch())
    y = similar(x, promote_type(eltype(P̂.M̂), eltype(x)), 2*P̂.M̂.N)
    P̂(y, x)
end

function _unsafe_mul!(y::AbstractVecOrMat, P̂::HelmholtzPreconditioner, x::AbstractVector)
    P̂(y, x)
end

# Update k methods

function update_k(M̂::HelmholtzMap{ND,T,NDp1,NDp2},k⃗::SVector{3,T2}) where {ND,T<:Real,NDp1,NDp2,T2<:Real}
	mag, mn = mag_mn(k⃗,M̂.g⃗)
	M̂.mag = mag
	M̂.inv_mag = [inv(mm) for mm in mag]
	M̂.mn = mn
	M̂.k⃗ = k⃗
end

function update_k(M̂::HelmholtzMap{ND,T1,NDp1,NDp2},kz::T2) where {ND,T1<:Real,NDp1,NDp2,T2<:Real}
	update_k(M̂,SVector{3,T}(0.,0.,kz))
end

update_k(ms::ModeSolver,k) = update_k(ms.M̂,k)

function update_k!(M̂::HelmholtzMap{ND,T1,NDp1,NDp2},k⃗::SVector{3,T2}) where {ND,T1<:Real,NDp1,NDp2,T2<:Real}
	# mag_m_n!(M̂.mag,M̂.m⃗,M̂.n⃗,k⃗,M̂.g⃗)
	M̂.k⃗ = k⃗
	mag_mn!(M̂.mag,M̂.mn,M̂.k⃗,M̂.g⃗)
	M̂.inv_mag = [inv(mm) for mm in M̂.mag]

end

function update_k!(M̂::HelmholtzMap{ND,T1,NDp1,NDp2},kz::T2) where {ND,T1<:Real,NDp1,NDp2,T2<:Real}
	update_k!(M̂,SVector{3,T2}(0.,0.,kz))
end

update_k!(ms::ModeSolver,k) = update_k!(ms.M̂,k)

# Update ε⁻¹ methods

function update_ε⁻¹(M̂::HelmholtzMap{ND,T,NDp1,NDp2},ε⁻¹) where {ND,T<:Real,NDp1,NDp2}
	@assert size(M̂.ε⁻¹) == size(ε⁻¹)
	M̂.ε⁻¹ = ε⁻¹
	M̂.ε_ave = [ 3. * inv(sum(diag(ε⁻¹[:,:,ix,iy]))) for ix in axes(ε⁻¹,3), iy in axes(ε⁻¹,3)]
end

function update_ε⁻¹(ms::ModeSolver{ND,T},ε⁻¹) where {ND,T<:Real}
	@assert size(ms.M̂.ε⁻¹) == size(ε⁻¹)
	ms.M̂.ε⁻¹ = ε⁻¹
	ms.M̂.ε_ave = [ 3. * inv(sum(diag(ε⁻¹[:,:,ix,iy]))) for ix=1:ms.grid.Nx, iy=1:ms.grid.Ny]
end

# property methods
Base.size(A::HelmholtzMap) = (2*A.N, 2*A.N)
Base.size(A::HelmholtzMap,d::Int) = 2*A.N
Base.eltype(A::HelmholtzMap{ND,T,NDp1,NDp2}) where {ND,T<:Real,NDp1,NDp2}  = Complex{T}
LinearAlgebra.issymmetric(A::HelmholtzMap) = false # A._issymmetric
LinearAlgebra.ishermitian(A::HelmholtzMap) = true # A._ishermitian
LinearAlgebra.isposdef(A::HelmholtzMap)    = true # A._isposdef
ismutating(A::HelmholtzMap) = true # A._ismutating
#_ismutating(f) = first(methods(f)).nargs == 3
Base.size(A::HelmholtzPreconditioner) = (2*A.M̂.N, 2*A.M̂.N)
Base.size(A::HelmholtzPreconditioner,d::Int) = 2*A.M̂.N
Base.eltype(A::HelmholtzPreconditioner) = eltype(A.M̂)
LinearAlgebra.issymmetric(A::HelmholtzPreconditioner) = true # A._issymmetric
LinearAlgebra.ishermitian(A::HelmholtzPreconditioner) = true # A._ishermitian
LinearAlgebra.isposdef(A::HelmholtzPreconditioner)    = true # A._isposdef
ismutating(A::HelmholtzPreconditioner) = true # A._ismutating

import Base: *, transpose, adjoint
function Base.:(*)(M::HelmholtzMap,X::Matrix)
	#if isequal(size(M),size(X)) # size check?
	ncolsX = size(X)[2]
	# @assert ncolsX == size(M)[1]
	Y = similar(X)
	for i in 1:ncolsX
		@views Y[:,i] = M * X[:,i]
	end
	return Y
end

function LinearAlgebra.mul!(y::AbstractVecOrMat, M̂::HelmholtzMap, x::AbstractVector)
    LinearMaps.check_dim_mul(y, M̂, x)
	M̂(y, x)
end

function LinearAlgebra.mul!(y::AbstractVecOrMat, P̂::HelmholtzPreconditioner, x::AbstractVector)
    LinearMaps.check_dim_mul(y, P̂, x)
	P̂(y, x)
end

Base.adjoint(A::HelmholtzMap) = A
Base.transpose(P̂::HelmholtzPreconditioner) = P̂
LinearAlgebra.ldiv!(c,P̂::HelmholtzPreconditioner,b) = mul!(c,P̂,b) # P̂(c, b) #
LinearAlgebra.ldiv!(P̂::HelmholtzPreconditioner,b) = mul!(b,P̂,b)

mag_m_n!(M̂::HelmholtzMap,k) = mag_m_n!(M̂.mag,M̂.m⃗,M̂.n⃗,M̂.g⃗,k)
mag_m_n!(ms::ModeSolver,k) = mag_m_n!(ms.M̂.mag,ms.M̂.m⃗,ms.M̂.n⃗,ms.M̂.g⃗,k)


text:
"""
    scan_rundir(pth::String)

Scan a MITgcm run directory and standard output text file
("output.txt" or "STDOUT.0000") and return a NamedTuple of
collected information (various formats)

Initially, the output looked like `(grid=gr,packages=pac,params_time=par1,params_grid=par2,completed=co)`
"""
function scan_rundir(pth::String)
    #1 standard output
    filout=joinpath(pth,"output.txt")
    !isfile(filout) ? filout=joinpath(pth,"STDOUT.0000") : nothing
    if isfile(filout)
        stdout=scan_stdout(filout)
    else
        stdout=missing
    end
    return stdout
end

"""
    scan_stdout(filout::String)

Scan a MITgcm standard output text file ("output.txt" or "STDOUT.0000") and return a NamedTuple of collected information (various formats).

- packages : report of packages being compiled and used
- params_time : initial time, model duation, output frequency, etc
- params_grid : type of grid (Curvilinear, Cartesian, ...) and array sizes
- params_files : type of output (use_mdsio, use_mnc) and array size (ioSize)
- completed : true / false depending on the end of standard output file (filout)
"""
function scan_stdout(filout::String)
    pth=dirname(filout)

    tmp=readlines(filout)
    #1.0 run completed
    co = tmp[end]=="PROGRAM MAIN: Execution ended Normally"

    ll = findall(occursin.("======================================================",tmp))
    #1.1 grid
    l0 = findall(occursin.("Computational Grid Specification",tmp))[1]
    l1 = minimum(ll[findall(ll.>l0+2)])
    gr = tmp[l0+2:l1]
    #1.2 packages
    l0 = findall(occursin.("PACKAGES_BOOT: On/Off package Summary",tmp))[1]
    l1 = findall(occursin.("PACKAGES_BOOT: End of package Summary",tmp))[1]
    pac = tmp[l0:l1]
    #1.3 parameters (time, grid)
    l0 = findall(occursin.("Time stepping paramters",tmp))[1]
    l1 = findall(occursin.("Gridding paramters",tmp))[1]
    l2 = findall(occursin.("End of Model config. summary",tmp))[1]

    par1=Dict()
    list1=["deltaTClock","nIter0","nTimeSteps","nEndIter",
            "pChkPtFreq","dumpFreq","monitorFreq"]
    for i in list1
        lx = findall(occursin.(i,tmp))
        lx = lx[findall((lx.>l0).*(lx.<l1))[1]]
        par1[Symbol(i)] = parse(Float64,tmp[lx+1][20:end])
    end
    par1=(; zip(Symbol.(keys(par1)), values(par1))...)

    par2=Dict()
    list1=["usingCartesianGrid","usingCylindricalGrid","usingSphericalPolarGrid","usingCurvilinearGrid"]
    for i in list1
        lx = findall(occursin.(i,tmp))
        lx = lx[findall((lx.>l1).*(lx.<l2))[1]]
        par2[Symbol(i)] = strip(tmp[lx+1][20:end])=="T"
    end

	b=["nPx","nPy","nSx","nSy","sNx","sNy","Nx","Ny","Nr"]
	for bb in b
		i1=findall(occursin.(bb,gr))[1]
		tmp=parse(Int,split(split(gr[i1][20:end],";")[1],"=")[2])
		par2[Symbol(bb)] = tmp
	end

    par2=(; zip(Symbol.(keys(par2)), values(par2))...)

    #1.4 monitors

    #2 output files

    #2.1 type (mdsio or mnc?) and overall size (ioSize) of output
    tst_mdsio = !isempty(filter(x -> occursin("XC",x), readdir(pth)))
    pth_mnc=joinpath(pth,"mnc_test_0001")
    tst_mnc = isdir(pth_mnc)
    ioSize=Array{Any}(undef,1)
    if tst_mdsio
        tmp=read_mdsio(pth,"XC")
        ioSize[1]=size(tmp)
    end
    if tst_mnc
        tmp=MITgcmTools.read_mnc(pth_mnc,"grid","XC")
        ioSize[1]=size(tmp)
    end

    par3=(use_mdsio=tst_mdsio,use_mnc=tst_mnc,ioSize=ioSize[1])

    #2.1 pickups
    #2.2 diags
    #2.3 snapshots
    #2.4 tave
    #2.5 mnc

    return (packages=pac,params_time=par1,params_grid=par2,params_files=par3,completed=co)
end

"""
    read_nctiles(fileName,fldName,mygrid)

Read model output from NCTiles file and convert to MeshArray instance.
```
mygrid=GridSpec("LatLonCap")
fileName="nctiles_grid/GRID"
Depth=read_nctiles(fileName,"Depth",mygrid)
hFacC=read_nctiles(fileName,"hFacC",mygrid)
hFacC=read_nctiles(fileName,"hFacC",mygrid,I=(:,:,1))
```
"""
function read_nctiles(fileName::String,fldName::String,mygrid::gcmgrid;
    I::Union{Missing,Tuple{Colon,Colon,Vararg{Union{Colon,Integer}}}}=missing)

    if (mygrid.class!="LatLonCap")||(mygrid.ioSize!=[90 1170])
        error("non-llc90 cases have not yet been tested with read_nctiles")
    end

    pth0=dirname(fileName)

    nam=split(fileName,"/")[end]
    isempty(nam) ? nam=split(fileName,"/")[end-1] : nothing
    occursin(".nctiles",nam) ? nam=nam[1:end-8] : nothing

    isdir(fileName) ? pth1=fileName : pth1=pth0
    lst=readdir(pth1)
    lst=lst[findall(occursin.(nam,lst).*occursin.(".nc",lst))]

    fileIn=joinpath(pth1,lst[1])
    fileRoot=fileIn[1:end-8]
    ntile=MITgcmTools.ncgetatt(fileIn,"Global","ntile")

    x = MITgcmTools.ncread(fileIn,fldName)
    s = [size(x,i) for i in 1:ndims(x)]
    n=length(size(x))
    start=ones(Int,n)
    count=-ones(Int,n)

    ~ismissing(I) && length(I)!=n ? error("ncdims v I inconsistency") : nothing
    if ~ismissing(I)
        k=findall([!isa(I[i],Colon) for i=1:length(I)])
        j=[I[i] for i in k]
        start[k]=j
        count[k].=1
        s[k].=1
    end

    f=Array{Float64, n}[]
    m0=[0]
    for ff in 1:mygrid.nFaces
        (ni,nj)=Int.(mygrid.fSize[ff]./s[1:2])
        nn=ni*nj
        i0=(mod1.(1:nn,ni).-1)*s[1]
        j0=div.(0:nn-1,ni)*s[2]

        #f0=Array{Float64}(undef,mygrid.fSize[ff]...,s[3:end]...)
        f0=fill(NaN,mygrid.fSize[ff]...,s[3:end]...)
        n0=m0[1]
        for n in 1:nn
            fileIn=@sprintf("%s.%04d.nc",fileRoot,n+n0)
            if isfile(fileIn) #skip if no file / blank tile
                x = ncread(fileIn,fldName,start,count)
                i=collect(1:s[1]) .+ i0[n]
                j=collect(1:s[2]) .+ j0[n]
                f0[i,j,:,:]=x[:,:,:,:]
            end
            m0[1]+=1
        end
        #f0[findall(isnan.(f0))].=0.0
        push!(f,f0)
    end

    fld=MeshArray(mygrid,f)
    return fld
end

"""
    findtiles(ni::Int,nj::Int,mygrid::gcmgrid)
    findtiles(ni::Int,nj::Int,grid::String="LatLonCap",GridParentDir="../inputs/GRID_LLC90/")

Return a `MeshArray` map of tile indices, `mytiles["tileNo"]`, for tile
size `ni,nj` and extract grid variables accordingly.
"""
function findtiles(ni::Int,nj::Int,mygrid::gcmgrid)
    mytiles = OrderedDict()

    GridVariables=GridLoad(mygrid)

    mytiles["nFaces"]=mygrid.nFaces;
    mytiles["ioSize"]=mygrid.ioSize;

    XC=GridVariables.XC;
    YC=GridVariables.YC;
    XC11=similar(XC); YC11=similar(XC);
    XCNINJ=similar(XC); YCNINJ=similar(XC);
    iTile=similar(XC); jTile=similar(XC); tileNo=similar(XC);
    tileCount=0;
    for iF=1:XC11.grid.nFaces
        face_XC=XC.f[iF]; face_YC=YC.f[iF];
        for jj=Int.(1:size(face_XC,2)/nj);
            for ii=Int.(1:size(face_XC,1)/ni);
                tileCount=tileCount+1;
                tmp_i=(1:ni).+ni*(ii-1)
                tmp_j=(1:nj).+nj*(jj-1)
                tmp_XC=face_XC[tmp_i,tmp_j]
                tmp_YC=face_YC[tmp_i,tmp_j]
                XC11.f[iF][tmp_i,tmp_j].=tmp_XC[1,1]
                YC11.f[iF][tmp_i,tmp_j].=tmp_YC[1,1]
                XCNINJ.f[iF][tmp_i,tmp_j].=tmp_XC[end,end]
                YCNINJ.f[iF][tmp_i,tmp_j].=tmp_YC[end,end]
                iTile.f[iF][tmp_i,tmp_j]=collect(1:ni)*ones(Int,1,nj)
                jTile.f[iF][tmp_i,tmp_j]=ones(Int,ni,1)*collect(1:nj)'
                tileNo.f[iF][tmp_i,tmp_j]=tileCount*ones(Int,ni,nj)
            end
        end
    end

    mytiles["tileNo"] = tileNo;
    mytiles["XC"] = XC;
    mytiles["YC"] = YC;
    mytiles["XC11"] = XC11;
    mytiles["YC11"] = YC11;
    mytiles["XCNINJ"] = XCNINJ;
    mytiles["YCNINJ"] = YCNINJ;
    mytiles["iTile"] = iTile;
    mytiles["jTile"] = jTile;

    return mytiles

end

findtiles(ni::Int,nj::Int,GridName::String="LatLonCap",GridParentDir="../inputs/GRID_LLC90/") = findtiles(ni,nj,GridSpec(GridName,GridParentDir))

## read_bin function with full list of argument

"""
    read_bin(fil::String,kt::Union{Int,Missing},kk::Union{Int,Missing},prec::DataType,mygrid::gcmgrid)

Read model output from binary file and convert to MeshArray. Other methods:

```
read_bin(fil::String,prec::DataType,mygrid::gcmgrid)
read_bin(fil::String,mygrid::gcmgrid)
```
"""
function read_bin(fil::String,kt::Union{Int,Missing},kk::Union{Int,Missing},prec::DataType,mygrid::gcmgrid)

  if ~ismissing(kt);
    error("non-empty kt option not implemented yet");
  end;

  if ~ismissing(kk);
    error("non-empty kk option not implemented yet");
  end;

  (n1,n2)=mygrid.ioSize

  if prec==Float64;
    reclen=8;
  else;
    reclen=4;
  end;
  tmp1=stat(fil);
  n3=Int64(tmp1.size/n1/n2/reclen);

  fid = open(fil);
  fld = Array{prec,1}(undef,(n1*n2*n3));
  read!(fid,fld);
  fld = hton.(fld);
  close(fid)

  n3>1 ? s=(n1,n2,n3) : s=(n1,n2)
  v0=reshape(fld,s);

  convert2gcmfaces(v0,mygrid)

end

## read_bin with reduced list of argument

# read_bin(fil::String,prec::DataType,mygrid::gcmgrid)
function read_bin(fil::String,prec::DataType,mygrid::gcmgrid)
  read_bin(fil,missing,missing,prec,mygrid::gcmgrid)
end

# read_bin(fil::String,mygrid::gcmgrid)
function read_bin(fil::String,mygrid::gcmgrid)
  read_bin(fil,missing,missing,mygrid.ioPrec,mygrid::gcmgrid)
end

## read_bin with alternative arguments

# read_bin(fil::String,x::MeshArray)
function read_bin(fil::String,x::MeshArray)
  read_bin(fil,missing,missing,eltype(x),x.grid::gcmgrid)
end

# read_bin(tmp::Array,mygrid::gcmgrid)
function read_bin(tmp::Array,mygrid::gcmgrid)
  convert2gcmfaces(tmp,mygrid)
end

# read_bin(tmp::Array,x::MeshArray)
function read_bin(tmp::Array,x::MeshArray)
  convert2gcmfaces(tmp,x.grid)
end

"""
    read_meta(metafile)

Read a `MITgcm` metadata file, parse it, and return as a NamedTuple

```
p="./hs94.cs-32x32x5/run/"
meta=read_meta(p*"surfDiag.0000000020.002.001.meta")
pairs(meta)
meta.dimList
```
"""
function read_meta(metafile)

    metafile[end-4:end]==".data" ? fil=metafile[1:end-5]*".meta" : fil=metafile[:]
    fil[end-4:end]!==".meta" ? error("inconsistent file name") : nothing

    meta = read(fil,String)
    meta = split(meta,";\n")
    meta = meta[isempty.(meta).==false]
    meta = replace.(meta,Ref(",\n"=>";"))
    meta = replace.(meta,Ref("\n"=>""))
    meta = replace.(meta,Ref("}"=>"]"))
    meta = replace.(meta,Ref("{"=>"["))
    #meta = replace.(meta,Ref(" "=>""))
    meta = replace.(meta,Ref("'"=>"\""))
    meta = replace.(meta,Ref(";]"=>"]"))
    meta = replace.(meta,Ref(","=>" "))

    meta = split.(meta,"=")
    meta = [[replace(x[1]," "=>"") x[2]] for x in meta]

    metaDict = OrderedDict{String,Any}(m[1] => m[2] for m in meta)

    for k in keys(metaDict)
        val = eval(Meta.parse(metaDict[k]))
        if isa(val[1],String)
            val = replace.(val,Ref(" "=>""))
        end
        if length(val) == 1
            val = val[1]
        end
        metaDict[k] = val
    end
    metaDict["dataprec"] = titlecase(metaDict["dataprec"])

    meta = (; zip(Symbol.(keys(metaDict)),values(metaDict))...)
    return meta

end

"""
    read_namelist(fil)

Read a `MITgcm` namelist file, parse it, and return as a NamedTuple

```
using MITgcmTools
testreport("advect_xy")
fil=joinpath(MITgcm_path[1],"verification","advect_xy","run","data")
namelist=read_namelist(fil)
```
"""
function read_namelist(fil)

    meta = read(fil,String)
    meta = split(meta,"\n")
    meta = meta[findall((!isempty).(meta))]
    meta = meta[findall(first.(meta).!=='#')]
    groups = meta[findall(occursin.('&',meta))]
	groups = [Symbol(groups[1+2*(i-1)][3:end]) for i in 1:Int(length(groups)/2)]
	params = fill(OrderedDict(),length(groups))

	for i in 1:length(groups)
		ii=1+findall(occursin.(String(groups[i]),meta))[1]
		i1=ii
		tmp0=OrderedDict()
        k0=[:unknown]
		while !occursin('&',meta[ii])
			if occursin('=',meta[ii])
				tmp1=split(meta[ii],'=')
                k0[1]=Symbol(strip(tmp1[1]))
				tmp2=split(tmp1[2],',')
                if length(tmp2)==2
                    tmp0[k0[1]]=strip(tmp2[1])
                else
                    tmp0[k0[1]]=strip(tmp1[2])
                end
            else
                try
                    tmp0[k0[1]]=tmp0[k0[1]]*","*strip(meta[ii])
                catch
                    println("ignoring line -- unclear why ...")
                end
			end
			ii += 1
		end
        for ii in keys(tmp0)
            tmp0[ii]=parse_param(tmp0[ii])
        end
		params[i]=tmp0
	end

    return MITgcm_namelist(Symbol.(groups),params)
end

"""
    parse_param(p1)

Parse namelist parameter and return in corresponding type
"""
function parse_param(p1)
	p2=missing
	if p1==".TRUE."||p1==".true."
		p2=true
	elseif p1==".FALSE."||p1==".false."
		p2=false
	else
        if first(p1)=='\''&&!occursin(',',p1)
			p2=p1[2:end-1]
        elseif occursin('.',p1)
			try
				p2=parse(Float64,p1)
			catch
				p2=p1
			end
		else
			try
				p2=parse(Int,p1)
			catch
				p2=p1
			end
		end
	end
    if isa(p2,AbstractString)
        if occursin(',',p2)
            p2=strip.(split(p2,','))
            p2=p2[findall( (!isempty).(p2) )]
			p2=[parse_param(p3) for p3 in p2]
        end
    end
	return p2
end

"""
    write_namelist(fil)

Save a `MITgcm` namelist file. In the example below, one is read from file, modified, and then saved to a new file using write_namelist.

```
using MITgcmTools
fil=joinpath(MITgcm_path[1],"verification","advect_xy","run","data")
nml=read_namelist(fil)
write_namelist(fil*"_new",namelist)
```

or

```
nml=read(fil,MITgcm_namelist())
write(fil*"_new",nml)
```
"""
function write_namelist(fil,namelist)
	fid = open(fil, "w")
	for jj in 1:length(namelist.groups)
        ii=namelist.groups[jj]
		tmpA=namelist.params[jj]
		params=(; zip(keys(tmpA),values(tmpA))...)

        txt=fill("",length(params))
        for i in 1:length(params)
            x=params[i]
            y=missing
            isa(x,Bool)&&x==true ? y=".TRUE." : nothing
            isa(x,Bool)&&x==false ? y=".FALSE." : nothing
            if isa(x,Array)&&(eltype(x)<:AbstractString)
                tmpy=[""]
                [tmpy[1]*="'"*x[ii]*"', \n " for ii in 1:length(x)]
                y=tmpy[1][1:end-4]
            elseif isa(x,Array)
                tmpy=[""]
                [tmpy[1]*=string(x[ii])*", \n " for ii in 1:length(x)]
                y=tmpy[1][1:end-4]
            end
            ismissing(y)&&isa(x,AbstractString)&&(!occursin('*',x)) ? y="'$x'" : nothing
            ismissing(y) ? y="$x" : nothing
            y[end]==',' ? y=y[1:end-1] : nothing
            txt[i]=y
        end

		txtparams=[" $(keys(params)[i]) = $(txt[i]),\n" for i in 1:length(params)]

		write(fid," &$(ii)\n")
		[write(fid,txtparams[i]) for i in 1:length(txtparams)]
		write(fid," &\n")
		write(fid," \n")
	end
	close(fid)
end

"""
    read_meta(pth::String,fil::String)

Read a `MITgcm` metadata files, parse them, and return as an array of NamedTuple

```
p="./hs94.cs-32x32x5/run/"
meta=read_meta(p,"surfDiag.0000000020")
pairs(meta[end])
[meta[i].dimList for i in 1:length(meta)]
```
"""
function read_meta(pth::String,fil::String)
    f=readdir(pth)
    kk=findall(occursin.(fil,f).*occursin.(".meta",f))
    [read_meta(pth*f[k]) for k in kk]
end

"""
    read_mdsio(fil::String)

Read a single `MITgcm` MDSIO-type file (".data" binary + ".meta" text pair), and return as an Array

```
p="./hs94.cs-32x32x5/run/"
x=read_mdsio(p*"surfDiag.0000000020.002.001.data")
y=read_mdsio(p*"pickup.ckptA.002.001.data")
z=read_mdsio(p*"T.0000000000.002.001.data")
```
"""
function read_mdsio(fil::String)
    m=read_meta(fil)
    T=eval(:($(Symbol(m.dataprec))))

    s=m.dimList[:,3]-m.dimList[:,2].+1
    isempty(findall(keys(m).==:nrecords)) ? n=1 : n=m.nrecords

    fid = open(fil)
    x = Array{T,1}(undef,(n*prod(s)))
    read!(fid,x)
    x = hton.(x)
    close(fid)

    n>1 ? s=Tuple([s;n]) : s=Tuple(s)
    x=reshape(x,s)

    return x
end

"""
    read_mdsio(pth::String,fil::String)

Read a set of `MITgcm` MDSIO-type files (".data" binary + ".meta" text pair), combine, and return as an Array.
Unlike `read_mdsio(fil::String)` where `fil` is one complete file name, this method will search within `pth`
for files that start with `fil`.

```
p="./hs94.cs-32x32x5/run/"
x=read_mdsio(p,"surfDiag.0000000020")
y=read_mdsio(p,"pickup.ckptA")
z=read_mdsio(p,"T.0000000000")
```
"""
function read_mdsio(pth::String,fil::String)
    f=readdir(pth)
    kk=findall(occursin.(fil,f).*occursin.(".data",f))
    f=f[kk]
    kk=findall([ (f[i][1:length(fil)]==fil)*(f[i][1+length(fil)]=='.') for i in 1:length(f)])

    m=[read_meta(joinpath(pth,f[k])) for k in kk]
    T=eval(:($(Symbol(m[1].dataprec))))

    m[1].nrecords>1 ? s=Tuple([m[1].dimList[:,1];m[1].nrecords]) : s=Tuple(m[1].dimList[:,1])
    x = Array{T,length(s)}(undef,s)

    for k=1:length(m)
        ii=m[k].dimList[1,2]:m[k].dimList[1,3]
        jj=m[k].dimList[2,2]:m[k].dimList[2,3]
        x[ii,jj,:,:]=read_mdsio(joinpath(pth,f[kk[k]]))
    end

    return x
end

#for the gcmgrid interface in MeshArrays.jl
function read_mdsio(fil::String,x::MeshArray)
    bas=split(basename(fil),'.')[1]
    xx=read_mdsio(dirname(fil),string(bas))
    return x.grid.read(xx,x)
end

#for the gcmgrid interface in MeshArrays.jl
read_mdsio(xx::Array,x::MeshArray) = MeshArrays.read(xx::Array,x::MeshArray)

"""
    read_mnc(pth::String,fil::String,var::String)

Read variable `var` from a set of `MITgcm` MNC-type files (netcdf files), combine, and
return as an Array. This method will search within `pth` for files that start with `fil`.
"""
function read_mnc(pth::String,fil::String,var::String)
    lst=readdir(pth)
    lst=lst[findall(occursin.(fil,lst).*occursin.(".nc",lst))]

    fil=joinpath(pth,lst[1])
    ncfile=MITgcmTools.NetCDF.open(fil)
    ncatts=ncfile.gatts

    v = MITgcmTools.NetCDF.open(fil, var)
    if haskey(v.atts,"coordinates")
        has_RC=occursin("RC",v.atts["coordinates"])
        has_iter=occursin("iter",v.atts["coordinates"])
    else
        has_RC=false
        has_iter=false
    end

    if has_RC*has_iter
        s=(ncatts["Nx"],ncatts["Ny"],Int64(v.dim[3].dimlen),Int64(v.dim[4].dimlen))
    elseif has_RC|has_iter
        s=(ncatts["Nx"],ncatts["Ny"],Int64(v.dim[3].dimlen))
    else
        s=(ncatts["Nx"],ncatts["Ny"])
    end

    T = Float64
    x = Array{T,length(s)}(undef,s)

    for f in lst
        fil=joinpath(pth,f)
        ncfile=MITgcmTools.NetCDF.open(fil)
        ncatts=ncfile.gatts
        b=(ncatts["bi"],ncatts["bj"])
        s=(ncatts["sNx"],ncatts["sNy"])
        ii=(b[1]-1)*s[1] .+ collect(1:s[1])
        jj=(b[2]-1)*s[2] .+ collect(1:s[2])
        v = MITgcmTools.NetCDF.open(fil, var)
        if has_RC*has_iter
            x[ii,jj,:,:]=v[:,:,:,:]
        elseif has_RC|has_iter
            x[ii,jj,:]=v[:,:,:]
        elseif has_RC|has_iter
            x[ii,jj]=v[:,:]
        end
    end

    x
end

"""
    GridLoad_mnc(γ::gcmgrid)

Load grid variabes (XC, YC, Depth) model run directory (`joinpath(rundir,"mnc_test_0001")`).
"""
function GridLoad_mnc(γ::gcmgrid)
    pth=joinpath(γ.path,"mnc_test_0001")
    XC=read_mnc(pth,"grid","XC")
    YC=read_mnc(pth,"grid","YC")
    Depth=read_mnc(pth,"grid","Depth")
    tmp=MeshArray(γ,γ.ioPrec)
    (XC=γ.read(XC,tmp), YC=γ.read(YC,tmp), Depth=γ.read(Depth,tmp))
end

"""
    GridLoad_mnc(myexp::MITgcm_config)

Load grid variables (XC, YC, Depth) from model run directory (`joinpath(rundir,"mnc_test_0001")`).
"""
function GridLoad_mnc(myexp::MITgcm_config)
    if isdir(joinpath(myexp.folder,string(myexp.ID),"run"))
        rundir=joinpath(myexp.folder,string(myexp.ID),"run")
    else
        pth=joinpath(MITgcm_path[1],"verification")
        rundir=joinpath(pth,myexp.configuration,"run")
    end
    GridLoad_mnc(rundir)
end

"""
    GridLoad_mnc(rundir::String)

Load grid variables (XC, YC, Depth) from model run directory (`joinpath(rundir,"mnc_test_0001")`).
"""
function GridLoad_mnc(rundir::String)

    data_mnc=joinpath(rundir,"data.mnc")
    if isfile(data_mnc)
        nml=read_namelist(data_mnc)
        pth=joinpath(rundir,nml.params[1][:mnc_outdir_str]*"0001")
    else
        pth=joinpath(rundir,"mnc_test_0001")
    end

	tmp=MITgcmTools.read_mnc(pth,"grid","XC")
    exps_ioSize=size(tmp)
    elty=eltype(tmp)
    sc=MITgcmTools.scan_rundir(rundir)
    #
    if sc.params_grid.usingCurvilinearGrid&&(exps_ioSize==(192,32))
        γ=gcmgrid(rundir,"CubeSphere",6,fill((32,32),6),[32 32*6],elty, read, write)
        c=cube2compact
    elseif sc.params_grid.usingCurvilinearGrid&&(exps_ioSize==(384,16))
        γ=gcmgrid(rundir,"CubeSphere",6,fill((32,32),6),[32 32*6],elty, read, write)
        function c(tmp)
            tmp2=reshape(tmp,(32,12,16))
            tmp3=Array{eltype(tmp2)}(undef,32*6,32)
            for i in 1:6
                ii=collect((i-1)*32 .+(1:32))
                tmp3[ii,1:16]=tmp2[:,2*i-1,:]
                tmp3[ii,17:32]=tmp2[:,2*i,:]
            end
            cube2compact(tmp3)
        end
    else
        s1=exps_ioSize
        s2=[s1[1] s1[2]]
        γ=gcmgrid(rundir,"PeriodicDomain",1,fill(s1,1),s2,elty, read, write)
        c=(x->x)
    end
    #
    pth=joinpath(γ.path,"mnc_test_0001")
    XC=c(read_mnc(pth,"grid","XC"))
    YC=c(read_mnc(pth,"grid","YC"))
    Depth=c(read_mnc(pth,"grid","Depth"))
    tmp=MeshArray(γ,γ.ioPrec)
    (XC=γ.read(XC,tmp), YC=γ.read(YC,tmp), Depth=γ.read(Depth,tmp))
end

"""
    GridLoad_mdsio(myexp::MITgcm_config)

Load grid variables (XC, YC, Depth, etc) from model run directory (`rundir`).
"""
function GridLoad_mdsio(myexp::MITgcm_config)
    if isdir(joinpath(myexp.folder,string(myexp.ID),"run"))
        rundir=joinpath(myexp.folder,string(myexp.ID),"run")
    else
        pth=joinpath(MITgcm_path[1],"verification")
        rundir=joinpath(pth,myexp.configuration,"run")
    end
    GridLoad_mdsio(rundir)
end

"""
    GridLoad_mdsio(rundir::String)

Load grid variables (XC, YC, Depth, etc) from model run directory (`rundir`).
"""
function GridLoad_mdsio(rundir::String)
    tmp=read_mdsio(rundir,"XC")
    exps_ioSize=size(tmp)
    elty=eltype(tmp)
    sc=MITgcmTools.scan_rundir(rundir)
    #
    if sc.params_grid.usingCurvilinearGrid
        readcube(xx::Array,x::MeshArray) = read_mdsio(cube2compact(xx),x)
        readcube(fil::String,x::MeshArray) = read_mdsio(fil::String,x::MeshArray)
        writecube(x::MeshArray) = compact2cube(write(x))
        writecube(fil::String,x::MeshArray) = write(fil::String,x::MeshArray)
        γ=gcmgrid(rundir,"CubeSphere",6,fill((32,32),6),[32 32*6],elty, readcube, writecube)
    else
        s1=exps_ioSize
        s2=[s1[1] s1[2]]
        γ=gcmgrid(rundir,"PeriodicDomain",1,fill(s1,1),s2,elty, read_mdsio, write)
    end
    Γ=GridLoad(γ;option="full")
end

"""
    read_available_diagnostics(fldname::String; filename="available_diagnostics.log")

Get the information for a particular variable `fldname` from the
`available_diagnostics.log` text file generated by `MITgcm`.
"""
function read_available_diagnostics(fldname::String; filename="available_diagnostics.log")
    availdiags = readlines(filename)
    line = availdiags[findall(occursin.(@sprintf("%-8s",fldname),availdiags))[1]]

    line = split(line,'|')
    line = lstrip.(rstrip.(line))

    diagInfo = OrderedDict([
    "diagNum" => parse(Int,line[1]),
    "fldname" => line[2],
    "levs" => parse(Int,line[3]),
    "mate" => line[4],
    "code" => line[5],
    "units" => line[6],
    "title" => line[7]
    ])

end

"""
    read_flt(dirIn::String,prec::DataType)

Read displacements from MITgcm/pkg/flt output file into a DataFrame.
"""
function read_flt(dirIn::String,prec::DataType)

   #load the data into one array
   prec==Float64 ? reclen=8 : reclen=4
   n1=13

   filIn="float_trajectories"
   tmp1=readdir(dirIn)
   tmp1=filter(x -> occursin(filIn,x),tmp1)
   filList=filter(x -> occursin(".data",x),tmp1)
   #hack:
   #filList=filter(x -> occursin("002.002.data",x),tmp1)
   nf=length(filList)

   n2=Array{Int,1}(undef,nf)
   for ff=1:nf
      fil=dirIn*filList[ff]
      #println(fil)
      tmp=stat(fil)
      n2[ff]=Int64(tmp.size/n1/reclen)-1
   end

   arr = Array{prec,2}(undef,(n1+1,sum(n2)));
   ii=0;
   #@softscope for ff=1:nf
   for ff=1:nf
      fil=dirIn*filList[ff]
      fid = open(fil)
      tmp = Array{prec,2}(undef,(n1,n2[ff]+1))
      read!(fid,tmp)
      arr[1:n1,ii+1:ii+n2[ff]] = hton.(tmp[:,2:n2[ff]+1])
      arr[n1+1,ii+1:ii+n2[ff]] .= ff
      ii=ii+n2[ff]
   end

   #sort the whole dataset by time
   jj = sort!([1:ii;], by=i->arr[2,i]); arr=arr[:,jj];
   #arr = sort!(arr, dims=2, by=i->arr[2,i]);

   #nfloats=Int(maximum(arr[1,:]))
   #npoints=counts(Int.(arr[1,:]))

   #reformat data as a DataFrame
   df=DataFrame()
   df.ID=Int.(arr[1,:])
   df.time=Int.(arr[2,:])
   df.lon=arr[3,:]
   df.lat=arr[4,:]
   df.dep=arr[5,:]
   if true
      df.i=arr[6,:]
      df.j=arr[7,:]
      df.k=arr[8,:]
      df.etaN=arr[9,:]
      df.uVel=arr[10,:]
      df.vVel=arr[11,:]
      df.theta=arr[12,:]
      df.salt=arr[13,:]
      df.tile=Int.(arr[14,:])
   end

#   nfloats=maximum(df.ID);
#   nsteps=maximum(counts(df.ID));
#   println("# floats=$nfloats")
#   println("# steps=$nsteps")

   return df
end

ext: jl


text:
using sb

@test sum(abs.(3*sb.σy - [[0 -3im]; [3im 0]])) < eps(1.0)
@test sum(abs.(sb.zm)) < eps(1.0)
@test sum(abs.(3*sb.σz - [[3 0]; [0 -3]])) < eps(1.0)
@test sum(abs.(3*sb.eye - [[3 0]; [0 3]])) < eps(1.0)
@test sum(abs.(3*sb.σx - [[0 3]; [3 0]])) < eps(1.0)


text:
@reexport module Utils
include("base.jl")
include("tiling.jl")
include("sum_tree.jl")
include("readers_writer_lock.jl")
include("parameter_server.jl")
include("devices.jl")
end


text:
"""
    GitLabCI(; config_file::Union{AbstractString, Nothing}="", coverage::Bool=true) -> GitLabCI

Add `GitLabCI` to a template's plugins to add a `.gitlab-ci.yml` configuration file to
generated repositories, and appropriate badge(s) to the README.

# Keyword Arguments:
* `config_file::Union{AbstractString, Nothing}=""`: Path to a custom `.gitlab-ci.yml`.
  If `nothing` is supplied, no file will be generated.
* `coverage::Bool=true`: Whether or not GitLab CI's built-in code coverage analysis should
  be enabled.
"""
struct GitLabCI <: GenericPlugin
    gitignore::Vector{String}
    src::Union{String, Nothing}
    dest::String
    badges::Vector{Badge}
    view::Dict{String, Any}

    function GitLabCI(; config_file::Union{AbstractString, Nothing}="", coverage::Bool=true)
        if config_file != nothing
            config_file = if isempty(config_file)
                config_file = joinpath(DEFAULTS_DIR, "gitlab-ci.yml")
            elseif isfile(config_file)
                abspath(config_file)
            else
                throw(ArgumentError("File $(abspath(config_file)) does not exist"))
            end
        end
        badges = [
            Badge(
                "Build Status",
                "https://gitlab.com/{{USER}}/{{PKGNAME}}.jl/badges/master/build.svg",
                "https://gitlab.com/{{USER}}/{{PKGNAME}}.jl/pipelines",
            ),
        ]
        if coverage
            push!(
                badges,
                Badge(
                    "Coverage",
                    "https://gitlab.com/{{USER}}/{{PKGNAME}}.jl/badges/master/coverage.svg",
                    "https://gitlab.com/{{USER}}/{{PKGNAME}}.jl/commits/master",
                ),
            )
        end

        new(
            coverage ? ["*.jl.cov", "*.jl.*.cov", "*.jl.mem"] : [],
            config_file,
            ".gitlab-ci.yml",
            badges,
            Dict("GITLABCOVERAGE" => coverage),
        )
    end
end

function interactive(::Type{GitLabCI})
    name = "GitLabCI"
    kwargs = Dict{Symbol, Any}()
    default_config_file = joinpath(DEFAULTS_DIR, "gitlab-ci.yml")

    print("$name: Enter the config template filename (\"None\" for no file) ")
    print("[", replace(default_config_file, homedir() => "~"), "]: ")
    config_file = readline()
    kwargs[:config_file] = if uppercase(config_file) == "NONE"
        nothing
    elseif isempty(config_file)
        default_config_file
    else
        config_file
    end

    print("$name: Enable test coverage analysis? [yes]: ")
    coverage = readline()
    kwargs[:coverage] = if isempty(coverage)
        true
    else
        !in(uppercase(coverage), ["N", "NO", "FALSE", "NONE"])
    end

    return GitLabCI(; kwargs...)
end


text:
# # Example 3: Pulse Parametrization

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`state_to_state_parametrizations.ipynb`](@__NBVIEWER_ROOT_URL__/examples/state_to_state_parametrizations.ipynb).

#md # ``\gdef\op#1{\hat{#1}}``
#md # ``\gdef\init{\text{init}}``
#md # ``\gdef\tgt{\text{tgt}}``

#nb # $
#nb # \newcommand{tr}[0]{\operatorname{tr}}
#nb # \newcommand{diag}[0]{\operatorname{diag}}
#nb # \newcommand{abs}[0]{\operatorname{abs}}
#nb # \newcommand{pop}[0]{\operatorname{pop}}
#nb # \newcommand{aux}[0]{\text{aux}}
#nb # \newcommand{opt}[0]{\text{opt}}
#nb # \newcommand{tgt}[0]{\text{tgt}}
#nb # \newcommand{init}[0]{\text{init}}
#nb # \newcommand{lab}[0]{\text{lab}}
#nb # \newcommand{rwa}[0]{\text{rwa}}
#nb # \newcommand{bra}[1]{\langle#1\vert}
#nb # \newcommand{ket}[1]{\vert#1\rangle}
#nb # \newcommand{Bra}[1]{\left\langle#1\right\vert}
#nb # \newcommand{Ket}[1]{\left\vert#1\right\rangle}
#nb # \newcommand{Braket}[2]{\left\langle #1\vphantom{#2}\mid{#2}\vphantom{#1}\right\rangle}
#nb # \newcommand{op}[1]{\hat{#1}}
#nb # \newcommand{Op}[1]{\hat{#1}}
#nb # \newcommand{dd}[0]{\,\text{d}}
#nb # \newcommand{Liouville}[0]{\mathcal{L}}
#nb # \newcommand{DynMap}[0]{\mathcal{E}}
#nb # \newcommand{identity}[0]{\mathbf{1}}
#nb # \newcommand{Norm}[1]{\lVert#1\rVert}
#nb # \newcommand{Abs}[1]{\left\vert#1\right\vert}
#nb # \newcommand{avg}[1]{\langle#1\rangle}
#nb # \newcommand{Avg}[1]{\left\langle#1\right\rangle}
#nb # \newcommand{AbsSq}[1]{\left\vert#1\right\vert^2}
#nb # \newcommand{Re}[0]{\operatorname{Re}}
#nb # \newcommand{Im}[0]{\operatorname{Im}}
#nb # $

# This example illustrates the parametrization of control pulses as a
# form of amplitude constraint.

using DrWatson
@quickactivate "KrotovTests"
#-
using QuantumControl
using QuantumControl.Shapes: flattop
using Krotov:
    SquareParametrization,
    TanhParametrization,
    TanhSqParametrization,
    LogisticParametrization,
    LogisticSqParametrization
using LinearAlgebra

#-
#jl using Test; println("")
using Plots
Plots.default(
    linewidth               = 3,
    size                    = (550, 300),
    legend                  = :top,
    foreground_color_legend = nothing,
    background_color_legend = RGBA(1, 1, 1, 0.8),
)
#-

# ## Parametrizations

# ## Symmetric Bounded Controls

#-
include(joinpath(@__DIR__, "plots", "symmetric_parametrization_comparison.jl"))  # hide
fig = plot_symmetric_parametrization_comparison()  # hide
#jl display(fig)

# ## Positive (Bounded) Controls

#-
include(joinpath(@__DIR__, "plots", "positive_parametrization_comparison.jl"))  # hide
fig = plot_positive_parametrization_comparison()  # hide
#jl display(fig)

# ## Two-level Hamiltonian

# We consider the Hamiltonian $\op{H}_{0} = - \frac{\omega}{2} \op{\sigma}_{z}$, representing
# a simple qubit with energy level splitting $\omega$ in the basis
# $\{\ket{0},\ket{1}\}$. The control field $\epsilon(t)$ is assumed to couple via
# the Hamiltonian $\op{H}_{1}(t) = \epsilon(t) \op{\sigma}_{x}$ to the qubit,
# i.e., the control field effectively drives transitions between both qubit
# states.
#
# We we will use

ϵ(t) = 0.2 * flattop(t, T=5, t_rise=0.3, func=:blackman);

#-
"""Two-level-system Hamiltonian."""
function hamiltonian(Ω=1.0, ϵ=ϵ)
    σ̂_z = ComplexF64[
        1  0
        0 -1
    ]
    σ̂_x = ComplexF64[
        0  1
        1  0
    ]
    Ĥ₀ = -0.5 * Ω * σ̂_z
    Ĥ₁ = σ̂_x
    return (Ĥ₀, (Ĥ₁, ϵ))
end;
#-

H = hamiltonian();
#jl @test length(H) == 2

# The control field here switches on from zero at $t=0$ to it's maximum amplitude
# 0.2 within the time period 0.3 (the switch-on shape is half a [Blackman pulse](https://en.wikipedia.org/wiki/Window_function#Blackman_window)).
# It switches off again in the time period 0.3 before the
# final time $T=5$). We use a time grid with 500 time steps between 0 and $T$:

tlist = collect(range(0, 5, length=500));

#-
function plot_control(pulse::Vector, tlist)
    plot(tlist, pulse, xlabel="time", ylabel="amplitude", legend=false)
end

plot_control(ϵ::T, tlist) where {T<:Function} = plot_control([ϵ(t) for t in tlist], tlist)

plot_control(H[2][2], tlist)

# ## Optimization target

# The `krotov` package requires the goal of the optimization to be described by a
# list of `Objective` instances. In this example, there is only a single
# objective: the state-to-state transfer from initial state $\ket{\Psi_{\init}} =
# \ket{0}$ to the target state $\ket{\Psi_{\tgt}} = \ket{1}$, under the dynamics
# of the Hamiltonian $\op{H}(t)$:

function ket(label)
    result = Dict("0" => Vector{ComplexF64}([1, 0]), "1" => Vector{ComplexF64}([0, 1]),)
    return result[string(label)]
end;

#-
#jl @test dot(ket(0), ket(1)) ≈ 0
#-

objectives = [Objective(initial_state=ket(0), generator=H, target_state=ket(1))]

#-
#jl @test length(objectives) == 1
#-

# ## Square-parametrization for positive pulses

problem = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ => Dict(
            :lambda_a => 5,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => SquareParametrization(),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence=res -> begin
        ((res.J_T < 1e-3) && (res.converged = true) && (res.message = "J_T < 10⁻³"))
    end
);
#-
opt_result_positive, file = @optimize_or_load(
    datadir(),
    problem;
    method=:krotov,
    filename="parametrization#opt_result_positive.jld2"
);
#-
opt_result_positive

# We can plot the optimized field:

#-
#!jl plot_control(opt_result_positive.optimized_controls[1], tlist)
#-

#-
#jl @test minimum(opt_result_positive.optimized_controls[1]) ≥ 0.0
#jl @test minimum(opt_result_positive.optimized_controls[1]) < 1e-16
#jl @test maximum(opt_result_positive.optimized_controls[1]) > 0.0
#-

# ## Tanh-Square-Parametrization for positive amplitude-constrained pulses

problem_tanhsq = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ => Dict(
            :lambda_a => 10,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => TanhSqParametrization(3),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence=res -> begin
        ((res.J_T < 1e-3) && (res.converged = true) && (res.message = "J_T < 10⁻³"))
    end
);
#-
opt_result_tanhsq, file = @optimize_or_load(
    datadir(),
    problem_tanhsq;
    method=:krotov,
    filename="parametrization#opt_result_tanhsq.jld2"
);
#-
opt_result_tanhsq

# We can plot the optimized field:

#-
#!jl plot_control(opt_result_tanhsq.optimized_controls[1], tlist)
#-

#-
#jl @test minimum(opt_result_tanhsq.optimized_controls[1]) ≥ 0.0
#jl @test minimum(opt_result_tanhsq.optimized_controls[1]) < 1e-16
#jl @test maximum(opt_result_tanhsq.optimized_controls[1]) > 0.0
#jl @test maximum(opt_result_tanhsq.optimized_controls[1]) < 3.0
#-

# ## Logistic-Square-Parametrization for positive amplitude-constrained pulses

problem_logisticsq = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ => Dict(
            :lambda_a => 1,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => LogisticSqParametrization(3, k=1.0),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence=res -> begin
        ((res.J_T < 1e-3) && (res.converged = true) && (res.message = "J_T < 10⁻³"))
    end
);
#-
opt_result_logisticsq, file = @optimize_or_load(
    datadir(),
    problem_logisticsq;
    method=:krotov,
    filename="parametrization#opt_result_logisticsq.jld2"
);
# We can plot the optimized field:

#-
#!jl plot_control(opt_result_logisticsq.optimized_controls[1], tlist)
#-

#-
#jl @test minimum(opt_result_logisticsq.optimized_controls[1]) ≥ 0.0
#jl @test minimum(opt_result_logisticsq.optimized_controls[1]) < 1e-16
#jl @test maximum(opt_result_logisticsq.optimized_controls[1]) > 0.0
#jl @test maximum(opt_result_logisticsq.optimized_controls[1]) < 3.0
#-
# ## Tanh-parametrization for amplitude-constrained pulses

problem_tanh = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        ϵ => Dict(
            :lambda_a => 1,
            :update_shape => t -> flattop(t, T=5, t_rise=0.3, func=:blackman),
            :parametrization => TanhParametrization(-0.5, 0.5),
        )
    ),
    tlist=tlist,
    iter_stop=50,
    J_T=QuantumControl.Functionals.J_T_ss,
    check_convergence=res -> begin
        ((res.J_T < 1e-3) && (res.converged = true) && (res.message = "J_T < 10⁻³"))
    end
);
#-
opt_result_tanh, file = @optimize_or_load(
    datadir(),
    problem_tanh;
    method=:krotov,
    filename="parametrization#opt_result_tanh.jld2"
);
#-
#!jl plot_control(opt_result_tanh.optimized_controls[1], tlist)
#-
#jl @test minimum(opt_result_tanh.optimized_controls[1]) > -0.5
#jl @test maximum(opt_result_tanh.optimized_controls[1]) < 0.5
#-


text:
module DecisionProgramming

include("influence_diagram.jl")
include("decision_model.jl")
include("random.jl")
include("analysis.jl")
include("printing.jl")

export Node,
    AbstractNode,
    ChanceNode,
    DecisionNode,
    ValueNode,
    State,
    States,
    Path,
    paths,
    Probabilities,
    Consequences,
    AbstractPathProbability,
    DefaultPathProbability,
    AbstractPathUtility,
    DefaultPathUtility,
    LocalDecisionStrategy,
    DecisionStrategy,
    validate_influence_diagram

export DecisionVariables,
    PathCompatibilityVariables,
    ForbiddenPath,
    lazy_probability_cut,
    PositivePathUtility,
    NegativePathUtility,
    expected_value,
    value_at_risk,
    conditional_value_at_risk

export random_diagram

export CompatiblePaths,
    UtilityDistribution,
    StateProbabilities

export print_decision_strategy,
    print_utility_distribution,
    print_state_probabilities,
    print_statistics,
    print_risk_measures

# For API docs
export AbstractRNG, Model, VariableRef

end # module


text:
using CUDA.APIUtils

@testset "@enum_without_prefix" begin
    mod = @eval module $(gensym())
        using CUDA.APIUtils
        @enum MY_ENUM MY_ENUM_VALUE
        @enum_without_prefix MY_ENUM MY_
    end

    @test mod.ENUM_VALUE == mod.MY_ENUM_VALUE
end
@testset "@checked" begin
    mod = @eval module $(gensym())
        using CUDA.APIUtils

        const checks = Ref(0)
        macro check(ex)
            esc(quote
                $checks[] += 1
                $ex
            end)
        end

        @checked function foo()
            ccall(:jl_getpid, Cint, ())
        end
    end

    @test mod.checks[] == 0
    @test mod.foo() == getpid()
    @test mod.checks[] == 1
    @test mod.unsafe_foo() == getpid()
    @test mod.checks[] == 1
end

@testset "@argout" begin
  f() = 1
  f(a) = 2
  f(a,b) = 3

  @test CUDA.@argout(f()) == nothing
  @test CUDA.@argout(f(out(4))) == 4
  @test CUDA.@argout(f(out(5), out(6))) == (5,6)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Objconv_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Objconv")
JLLWrappers.@generate_main_file("Objconv", UUID("dfe201cf-525d-5101-8f66-6dc5e2f3b224"))
end  # module Objconv_jll


text:
# Copyright 2013 Peter Chronz
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

export
    # Typealiases
    Bytes,

    # Types
    OfpStruct,
    OfpHeader,
    OfpMessage,
    OfpQueueGetConfigRequest,
    OfpQueueProp,
    OfpQueuePropHeader,
    OfpQueuePropNone,
    UnknownQueueProperty,
    OfpQueuePropMinRate,
    OfpPacketQueue,
    OfpQueueGetConfigReply,
    OfpError,
    UnrecognizedErrorType,
    UnrecognizedErrorCode,
    OfpPhyPort,
    OfpSwitchFeatures,
    OfpSwitchConfig,
    OfpMatch,
    OfpActionHeader,
    OfpActionEmpty,
    OfpActionOutput,
    OfpActionEnqueue,
    OfpActionVlanVid,
    OfpActionVlanPcp,
    WrongMessageType,
    MessageInvarianceException,
    OfpActionDlAddress,
    OfpActionNwAddress,
    OfpActionNwTos,
    OfpActionTpPort,
    OfpActionVendor,
    OfpPacketIn,
    OfpPortStatus,
    OfpFlowMod,
    OfpPortMod,
    OfpEmptyMessage,
    UnrecognizedMessageError,
    OfpStatsRequestBody,
    OfpStatsReplyBody,
    OfpDescStats,
    OfpFlowStatsRequest,
    OfpFlowStats,
    OfpAggregateStatsRequest,
    OfpAggregateStatsReply,
    OfpTableStats,
    OfpPortStatsRequest,
    OfpPortStats,
    OfpQueueStatsRequest,
    OfpQueueStats,
    OfpVendorStatsRequest,
    OfpVendorStatsReply,
    OfpStatsRequest,
    OfpStatsReply,
    OfpPacketOut,
    OfpFlowRemoved,
    OfpVendorHeader,

    # Constants.
    OFPT_HELLO,
    OFPT_ERROR,
    OFPT_ECHO_REQUEST,
    OFPT_ECHO_REPLY,
    OFPT_VENDOR,
    OFPT_FEATURES_REQUEST,
    OFPT_FEATURES_REPLY,
    OFPT_GET_CONFIG_REQUEST,
    OFPT_GET_CONFIG_REPLY,
    OFPT_SET_CONFIG,
    OFPT_PACKET_IN,
    OFPT_FLOW_REMOVED,
    OFPT_PORT_STATUS,
    OFPT_PACKET_OUT,
    OFPT_FLOW_MOD,
    OFPT_PORT_MOD,
    OFPT_STATS_REQUEST,
    OFPT_STATS_REPLY,
    OFPT_BARRIER_REQUEST,
    OFPT_BARRIER_REPLY,
    OFPT_QUEUE_GET_CONFIG_REQUEST,
    OFPT_QUEUE_GET_CONFIG_REPLY,
    OFPPC_PORT_DOWN,
    OFPPC_NO_STP,
    OFPPC_NO_RECV,
    OFPPC_NO_RECV_STP,
    OFPPC_NO_FLOOD,
    OFPPC_NO_FWD,
    OFPPC_NO_PACKET_IN,
    OFPPS_LINK_DOWN,
    OFPPS_STP_LISTEN,
    OFPPS_STP_LEARN,
    OFPPS_STP_FORWARD,
    OFPPS_STP_BLOCK,
    OFPPS_STP_MASK,
    OFPP_MAX,
    OFPP_IN_PORT,
    OFPP_TABLE,
    OFPP_NORMAL,
    OFPP_FLOOD,
    OFPP_ALL,
    OFPP_CONTROLLER,
    OFPP_LOCAL,
    OFPP_NONE,
    OFPPF_10MB_HD,
    OFPPF_10MB_FD,
    OFPPF_100MB_HD,
    OFPPF_100MB_FD,
    OFPPF_1GB_HD,
    OFPPF_1GB_FD,
    OFPPF_10GB_FD,
    OFPPF_COPPER,
    OFPPF_FIBER,
    OFPPF_AUTONEG,
    OFPPF_PAUSE,
    OFPPF_PAUSE_ASYM,
    OFPC_FLOW_STATS,
    OFPC_TABLE_STATS,
    OFPC_PORT_STATS,
    OFPC_STP,
    OFPC_RESERVED,
    OFPC_IP_REASM,
    OFPC_QUEUE_STATS,
    OFPC_ARP_MATCH_IP,
    OFPC_FRAG_NORMAL,
    OFPC_FRAG_DROP,
    OFPC_FRAG_REASM,
    OFPC_FRAG_MASK,
    OFPR_NO_MATCH,
    OFPR_ACTION,
    OFPPR_ADD,
    OFPPR_DELETE,
    OFPPR_MODIFY,
    OFPFW_IN_PORT,
    OFPFW_DL_VLAN,
    OFPFW_DL_SRC,
    OFPFW_DL_DST,
    OFPFW_DL_TYPE,
    OFPFW_NW_PROTO,
    OFPFW_TP_SRC,
    OFPFW_TP_DST,
    OFPFW_NW_SRC_SHIFT,
    OFPFW_NW_SRC_BITS,
    OFPFW_NW_SRC_MASK,
    OFPFW_NW_SRC_ALL,
    OFPFW_NW_DST_SHIFT,
    OFPFW_NW_DST_BITS,
    OFPFW_NW_DST_MASK,
    OFPFW_NW_DST_ALL,
    OFPFW_DL_VLAN_PCP,
    OFPFW_NW_TOS,
    OFPFW_ALL,
    OFPAT_OUTPUT,
    OFPAT_SET_VLAN_VID,
    OFPAT_SET_VLAN_PCP,
    OFPAT_STRIP_VLAN,
    OFPAT_SET_DL_SRC,
    OFPAT_SET_DL_DST,
    OFPAT_SET_NW_SRC,
    OFPAT_SET_NW_DST,
    OFPAT_SET_NW_TOS,
    OFPAT_SET_TP_SRC,
    OFPAT_SET_TP_DST,
    OFPAT_ENQUEUE,
    OFPAT_VENDOR,
    OFPFC_ADD,
    OFPFC_MODIFY,
    OFPFC_MODIFY_STRICT,
    OFPFC_DELETE,
    OFPFC_DELETE_STRICT,
    OFPFF_SEND_FLOW_REM,
    OFPFF_CHECK_OVERLAP,
    OFPFF_EMERG,
    OFPET_HELLO_FAILED,
    OFPET_BAD_REQUEST,
    OFPET_BAD_ACTION,
    OFPET_FLOW_MOD_FAILED,
    OFPET_PORT_MOD_FAILED,
    OFPET_QUEUE_OP_FAILED,
    OFPBRC_BAD_VERSION,
    OFPBRC_BAD_TYPE,
    OFPBRC_BAD_STAT,
    OFPBRC_BAD_VENDOR,
    OFPBRC_BAD_SUBTYPE,
    OFPBRC_EPERM,
    OFPBRC_BAD_LEN,
    OFPBRC_BUFFER_EMPTY,
    OFPBRC_BUFFER_UNKNOWN,
    OFPBAC_BAD_TYPE,
    OFPBAC_BAD_LEN,
    OFPBAC_BAD_VENDOR,
    OFPBAC_BAD_VENDOR_TYPE,
    OFPBAC_BAD_OUT_PORT,
    OFPBAC_BAD_ARGUMENT,
    OFPBAC_EPERM,
    OFPBAC_TOO_MANY,
    OFPBAC_BAD_QUEUE,
    OFPFMFC_ALL_TABLES_FULL,
    OFPFMFC_OVERLAP,
    OFPFMFC_EPERM,
    OFPFMFC_BAD_EMERG_TIMEOUT,
    OFPFMFC_BAD_COMMAND,
    OFPFMFC_UNSUPPORTED,
    OFPMFC_BAD_PORT,
    OFPMFC_BAD_HW_ADDR,
    OFPQOFC_BAD_PORT,
    OFPQOFC_BAD_QUEUE,
    OFPQOFC_EPERM,
    OFPQT_NONE,
    OFPQT_MIN_RATE,
    OFPST_DESC,
    OFPST_FLOW,
    OFPST_AGGREGATE,
    OFPST_TABLE,
    OFPST_PORT,
    OFPST_QUEUE,
    OFPST_VENDOR,
    OFPRR_IDLE_TIMEOUT,
    OFPRR_HARD_TIMEOUT,
    OFPRR_DELTE,
    OFP_MAX_ETH_ALEN,
    OFP_MAX_PORT_NAME_LEN,
    DESC_STR_LEN,
    SERIAL_NUM_LEN,
    OFP_MAX_TABLE_NAME_LEN,

    # functions
    give_length,
    string,
    bytes,
    btoui,

    # API
    start_server

max_stars_repo_name: pchronz/juliaflow

max_issues_repo_name: pchronz/juliaflow

max_forks_repo_name: pchronz/juliaflow


text:
Sonum8(x::UInt8) = reinterpret(Sonum8,x)
Sonum16(x::UInt16) = reinterpret(Sonum16,x)

UInt8(x::Sonum8) = reinterpret(UInt8,x)
UInt16(x::Sonum16) = reinterpret(UInt16,x)

Int(x::Sonum8) = Int(UInt8(x))
Int(x::Sonum16) = Int(UInt16(x))


text:
using QuadGK
using StaticArrays

abstract type AbstractIntegrateFixedOrder end

struct IntegrateFixedOrder{FuncT,Order,NumGW,RetT} <: AbstractIntegrateFixedOrder  where { FuncT, Order<:Integer, RetT<:Real }
    x::SVector{Order,RetT}  # Technically, order = Order-1
    w::SVector{Order,RetT}
    gw::SVector{NumGW,RetT}
    f::FuncT
end

function IntegrateFixedOrder(f::FuncT,order::Integer)  where { FuncT } #, RetT }
    @assert 1<=order <=10
    #ret_type = RetT
    ret_type = typeof(f(0.0))
    x,w,gw = QuadGK.cachedrule(ret_type,order)
    s_x = SVector{order+1,ret_type}(x)
    s_w = SVector{order+1,ret_type}(w)
    num_gw = floor(Int,(order+1)//2)
    s_gw = SVector{num_gw,ret_type}(gw)
    #s_gw = @SVector zeros(floor(Int,(order+1)//2))
    #s_gw[1:length(gw)] .= gw   # since s_gw is generally too big

    IntegrateFixedOrder{FuncT,order+1,num_gw,ret_type}(s_x,s_w,s_gw,f)
end

function integrate_scalar_fixed_order_gauss(ifo::T, a::Real, b::Real)  where { T<:AbstractIntegrateFixedOrder }
    s = convert(eltype(ifo.x), 0.5) * (b-a)
    n1 = 1 - (length(ifo.x) & 1) # 0 if even order, 1 if odd order
    # unroll first iterationof loop to get correct type of Ik and Ig
    fg = ifo.f(a + (1+ifo.x[2])*s) + ifo.f(a + (1-ifo.x[2])*s)
    Ig = fg * ifo.gw[1]
    for i = 2:length(ifo.gw)-n1
        fg = ifo.f(a + (1+ifo.x[2i])*s) + ifo.f(a + (1-ifo.x[2i])*s)
        Ig += fg * ifo.gw[i]
    end
    if n1 != 0 # even: Gauss rule does not include x == 0
        f0 = ifo.f(a + s)
        Ig += f0 * ifo.gw[end]
    end
    Ig_s = Ig * s # new variable since this may change the type
    return Ig_s
end

function integrate_scalar_fixed_order_kronrod(ifo::T, a::Real, b::Real)  where { T<:AbstractIntegrateFixedOrder }
    s = convert(eltype(ifo.x), 0.5) * (b-a)
    n1 = 1 - (length(ifo.x) & 1) # 0 if even order, 1 if odd order
    # unroll first iterationof loop to get correct type of Ik and Ig
    fg = ifo.f(a + (1+ifo.x[2])*s) + ifo.f(a + (1-ifo.x[2])*s)
    fk = ifo.f(a + (1+ifo.x[1])*s) + ifo.f(a + (1-ifo.x[1])*s)
    Ik = fg * ifo.w[2] + fk * ifo.w[1]
    for i = 2:length(ifo.gw)-n1
        fg = ifo.f(a + (1+ifo.x[2i])*s) + ifo.f(a + (1-ifo.x[2i])*s)
        fk = ifo.f(a + (1+ifo.x[2i-1])*s) + ifo.f(a + (1-ifo.x[2i-1])*s)
        Ik += fg * ifo.w[2i] + fk * ifo.w[2i-1]
    end
    if n1 == 0 # even: Gauss rule does not include x == 0
        Ik += f(a + s) * ifo.w[end]
    else # odd: don't count x==0 twice in Gauss rule
        f0 = f(a + s)
        Ik += f0 * ifo.w[end] +
            (f(a + (1+ifo.x[end-1])*s) + f(a + (1-ifo.x[end-1])*s)) * ifo.w[end-1]
    end
    Ik_s = Ik * s # new variable since this may change the type
    return Ik_s
end

function (ifo::IntegrateFixedOrder)(a::Real, b::Real)
    integrate_scalar_fixed_order_gauss(ifo,a,b)
end

function error_estimate(ifo::T, a::Real, b::Real)  where { T<:AbstractIntegrateFixedOrder }
    Ig = integrate_scalar_fixed_order_gauss(ifo,a,b)
    Ik = integrate_scalar_fixed_order_kronrod(ifo,a,b)
    return abs(Ik - Ig)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
"""
# Description

Rewrite an expression to remove all use of backticks.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `e::Any`: An expression in which backticks have been removed.

# Examples

```
julia> remove_backticks(:(`mean(a)`))
:(mean(a))
```
"""
function remove_backticks(@nospecialize(e::Any))
    if isa(e, Expr) && e.head == :macrocall && isa(e.args[1], GlobalRef) && e.args[1].name == Symbol("@cmd")
        Meta.parse(e.args[3])
    elseif isa(e, Expr)
        Expr(
            e.head,
            [remove_backticks(a) for a in e.args]...,
        )
    else
        e
    end
end