# Vectors defining a plane source: the plane normal and the direction of dipole current.

# PlaneSrc supports the polarization normal to the plane.  In that case, the source strength
# is interpreted as J∆n rather than sheet current density K.

# For the TE equations, the electric plane source has K = 2 and Kf = 2, whereas the magnetic
# plane source has K = 2 and Kf = 1.  This might be nonintuitive to the users, but it
# becomes obvious when the equations are written down: in the TE equations on the xy-plane,
# the electric current source has both the x- and y-components (so Kf = 2), whereas the
# magnetic current source has only the z-component (so Kf = 1).

export PlaneSrc

mutable struct PlaneSrc{K,Kf} <: Source{K,Kf}
    # Specify geometry first, current next (because the value of current is optional).
    n::SFloat{K}  # direction normal to plane; currently supports only Cartesian directions
    c::Float  # intercept (c means center)
    p::SFloat{Kf}  # polarization direction in field dimension
    J∆n::ComplexF # sheet current density (current per unit in-plane length normal to flow direction); complex in order to represent phase difference between plane sources
    isfield˔shp::Bool  # true if spaces where field and shapes are defined are orthogonal complement of each other; false if they are the same

    PlaneSrc{K,Kf}(n::AbsVecReal, c::Real, p::AbsVecReal, J∆n::Number, isfield˔shp::Bool) where {K,Kf} = new(normalize(n), c, normalize(p), J∆n, isfield˔shp)
end

function PlaneSrc(n::AbsVecReal, c::Real, p::AbsVecReal, J∆n::Number=1.0, isfield˔shp::Bool=isfield_ortho_shape(length(p),length(n)))
    count(n.!=0) == 1 || throw(ArgumentError("n = $n must be along Cartesian direction."))

    K = length(n)
    Kf = length(p)

    return PlaneSrc{K,Kf}(n, c, p, J∆n, isfield˔shp)
end

function add_src!(jKd::AbsArrNumber{K₊₁},  # (K+1)-dimensional array of Je (electric) or Jm (magnetic); first K dimensions specify location; last 1 dimension specify field component
                  gt₀::SVec{K,GridType},  # grid type of voxel corners; generated by ft2gt.(ft, boundft)
                  bounds::Tuple2{SFloat{K}},  # bounds[NEG][k] = boundary of domain at negative end in k-direction
                  l::Tuple2{NTuple{K,VecFloat}},  # l[PRIM][k] = primal vertex locations in k-direction
                  ∆l::Tuple2{NTuple{K,VecFloat}},  # ∆l[PRIM][k] = (∆l at primal vertices in w) == diff(l[DUAL][k] including ghost point)
                  isbloch::SBool{K},  # Bloch boundary conditions
                  src::PlaneSrc{K,Kf}  # plane source to add
                  ) where {K,K₊₁,Kf}
    nn = findfirst(src.n.==1)  # index of Cartesian direction corresponding to plane normal
    for nw = 1:Kf  # set Jw (w-component of J)
        gt_cmp = src.isfield˔shp ? gt₀ : gt_w(nw, gt₀)  # grid type of Jw

        lcmp = t_ind(l, gt_cmp)
        ∆lcmp = t_ind(∆l, gt_cmp)

        gtn = gt_cmp[nn]
        ln = lcmp[nn]
        ∆ln = ∆lcmp[nn]
        # Above, ln is the field locations in the plane normal direction.
        # - If src.isfield˔shp = true, then gt_cmp = gt₀ for all nw, so ln does not change
        # with nw.
        # - If src.isfield˔shp = false, then gt_cmp changes with nw, but ln changes only
        # when nw = nn, because lcmp[nn], which is determined by gt_cmp[nn], remains gt₀[nn]
        # for nw ≠ nn and changes to alter(gt₀[nn]) only for nw = nn.
        # Using the above observation, it would be possible to simplify the above code, but
        # I leave it the way it is in order to maintain the similarity of the implementation
        # between concrete Source types.

        ind, wt = distweights(src.c, gtn, t_ind(bounds,nn,nn), ln, ∆ln, isbloch[nn])

        Jw∆n = src.J∆n * src.p[nw]  # nw is @inbounds, because src.p is length-Kf and 1 ≤ nw ≤ Kf
        Jw = Jw∆n * wt  # Jw has two elements like wt

        # The code below mimics the implementation of slicedim.
        # For example, V[Base.setindex(axes(V), ir:ir, nr)...] means V[:,ir:ir,:] for r = y.
        # (axes(V) = (1:Nx,1:Ny,1:Nz) and setindex(t,v,i) does t[i] = v.  Therefore,
        # Base.setindex(axes(V), ir:ir, nr) creates (1:Nx, ir:ir, 1:Nz) for r = y.)
        # The use of ir:ir instead of ir is to support 1D arrays.  If V is 1D, then V[ir] is
        # a scalar and the dot equal on V[ir] fails.
        for k = 1:2  # 2 == length(ind)
            # Below, [1:K] is used to take the first K entries of the length-K₊₁ tuple.
            jKd[Base.setindex(axes(jKd), ind[k]:ind[k], nn)[1:K]..., nw] .+= Jw[k]
        end
    end

    return nothing
end
