# Vectors defining a point source: the direction of the dipole current.

export PointSrc

# When put a point source in space discretized by a finite-difference grid, the point source
# must specify current dipole I∆r rather current I or current density J.
#
# - If it assigns current density J, it is assumed to flow uniformly through the area
# element ∆a where J is centered, so the total current is I = J ∆a, which decreases with
# increasing resolution.  If J is electric current density and if it decreases into a half,
# the induced magnetic field that circuits around I decreases into a half at the same
# distance from I.  Therefore, in order to keep the simulation result independent of grid
# resolution for sufficiently high resolution (i.e., in order to make the simulation result
# converges), current I instead of current density J should be assigned.
#
# - Now, if it assigns current I, the current is assumed to flow uniformly in the edge
# element ∆l along the direction of the current.  (The fact that this assumption is made in
# finite difference is evident if we consider a loop of electric current formed along the
# perimeter of a pixel: the magnetic current dipole induced by this electric current loop is
# proportional to the line integral of the current, and on the finite-difference grid, this
# line integral is calculated by assuming each edge has a uniform current over its length.)
# For the same current I, if ∆l decreases, the current dipole I ∆l decreases.  When the
# current dipole decreases, the E-field induced by the dipole decreases.
#
# In conclusion, a point source must specify the dipole current I∆r, and from that it must
# recover J = I∆r / (∆a ∆l) to assign to the field point.
#
# Wait, if we put a z-directional electric current component in a system that is periodic
# and one-cell-thick in the z-direction, shouldn't Jz be independent of the thickness ∆z of
# the system?  If we put a point source described above, then changing ∆z will change Jz to
# assign.  This is not what we want in general in simulating a 2D system.
#
# This means that a point source is not the correct source type to simulate a 2D system.
# This is obvious if we think about the actual physical systems represented by discretized
# systems.  The physical source distribution represented by the above-described situation is
# a periodic array of point sources separated by ∆z.  Here, by reducing ∆z while keeping the
# strength of each point source, we actually change the physical situation.  We have more
# charges per unit z oscillating along the z-direction; in other words, the charge density
# per unit z along the source line increases.  This change of charge density with varying ∆z
# should be captured in the discretized system, and the above-described point source
# specification in terms of the current dipole I∆r does so: as ∆z decreases, Iz = Jz ∆x ∆y
# = I∆r / ∆z increases indeed.
#
# On the other hand, the similar 2D situation we usually want to simulate is a 3D system
# with a z-directional line source with a constant current Iz.  Here, choosing ∆z is just to
# pick a certain ∆z-long section of the current, in which the current Iz does not change
# with the choice of ∆z.  Therefore, we must have a line source type that is specified with
# a current I instead of a current dipole I∆r.  Using such a line source would not change
# Jz = I / (∆x ∆y) to assign on each edge element even if ∆z changes.

mutable struct PointSrc{K,Kf} <: Source{K,Kf}
    # Specify geometry first, current next (because the value of current is optional).
    c::SFloat{K}  # location in shape dimension (c menas center)
    p::SFloat{Kf}  # polarization direction in field dimension
    I∆r::ComplexF  # source strength (= dipole current I∆r); complex in order to represent phase difference between point 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

    PointSrc{K,Kf}(c::AbsVecReal, p::AbsVecReal, I∆r::Number, isfield˔shp::Bool) where {K,Kf} = new(c, normalize(p), I∆r, isfield˔shp)
end

PointSrc(c::AbsVecReal, p::AbsVecReal, I∆r::Number=1.0; isfield˔shp::Bool=isfield_ortho_shape(length(p),length(c))) =
    (K = length(c); Kf = length(p); PointSrc{K,Kf}(c, p, I∆r, isfield˔shp))

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::PointSrc{K,Kf}  # point source to add
                  ) where {K,K₊₁,Kf}
    @assert K₊₁==K+1

    for nw = 1:Kf  # set Jw (w-component of J)
        gt_cmp = src.isfield˔shp ? gt₀ : gt_w(nw, gt₀)  # grid type of location of Jw

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

        ind_cmp = ntuple(k->@SVector(zeros(Int,2)), Val(K))  # ind_cmp[nu]: indices of voxel points in nu-direction
        wt_vxl = SArray{Tuple{ntuple(k->2,Val(K))...}}(ntuple(k->1.0, Val(2^K)))  # 2(-by-2)^(K-1) array of weight factors
        for nu = 1:K  # u-coordinate of location
            u, ∆u = lcmp[nu], ∆lcmp[nu]
            ind_u, wt_u = distweights(src.c[nu], gt_cmp[nu], t_ind(bounds,nu,nu), u, ∆u, isbloch[nu])

            # Using ntuple(), replace nu-entry of ind_vxl to create a new tuple
            wt_vxl = wt_vxl .* SArray{Tuple{ntuple(k->(k==nu ? 2 : 1),Val(K))...}}(wt_u)
            ind_cmp = ntuple(k->(k==nu ? ind_u : ind_cmp[k]), Val(K))
        end

        # Calculate the w-component (I∆r)_w of the dipole current I∆r.  Once the w-component
        # of the dipole current is obtained, we obtain the w-component Iw of the current
        # by dividing (I∆r)_w by ∆w.  Therefore, (I∆r)_w / ∆w = Iw, or Iw∆w = (I∆r)_w.
        Iw∆w = src.I∆r * src.p[nw]
        Jw = Iw∆w .* wt_vxl  # Jw = Iw / (∆u ∆v) = (Iw ∆w / ∆w) / (∆u ∆v) = (Iw ∆w) / (∆u ∆v ∆w)
        for ci = CartesianIndices(ntuple(k->1:2, Val(K)))
            ind_cnr = t_ind(ind_cmp, ci)  # SInt{K}: index of voxel corner
            jKd[ind_cnr.data..., nw] += Jw[ci]
        end
    end

    return nothing
end


# # The above function for an arbitrary K is the generalization of the following function for
# # K = 3.
# function add_src!(jKd::AbsArrNumber{4},  # (3+1)-dimensional array of Je (electric) or Jm (magnetic); first 3 dimensions specify location; last 1 dimension specify field component
#                   gt₀::SVec{3,GridType},  # grid type of voxel corners; generated by ft2gt.(ft, boundft)
#                   bounds::Tuple2{SFloat{3}},  # bounds[NEG][k] = boundary of domain at negative end in k-direction
#                   l::Tuple23{VecFloat},  # l[PRIM][k] = primal vertex locations in k-direction
#                   ∆l::Tuple23{VecFloat},  # ∆l[PRIM][k] = (∆l at primal vertices in w) == diff(l[DUAL][k] including ghost point)
#                   isbloch::SBool{3},  # Bloch boundary conditions
#                   src::PointSrc{3,3})  # point source to add
#     for nw = 1:3  # set Jw (w-component of J)
#         Iw∆w = src.I∆r * src.p[nw]  # w-component of I∆r
#
#         gt_cmp = gt_w(nw, gt₀)  # grid type of Jw
#
#         lcmp = t_ind(l, gt_cmp)
#         ∆lcmp = t_ind(∆l, gt_cmp)
#
#         x, ∆x = lcmp[1], ∆lcmp[1]
#         indx, wtx = distweights(src.c[1], gt_cmp[1], t_ind(bounds,1,1), x, ∆x, isbloch[1])
#
#         y, ∆y = lcmp[2], ∆lcmp[2]
#         indy, wty = distweights(src.c[2], gt_cmp[2], t_ind(bounds,2,2), y, ∆y, isbloch[2])
#
#         z, ∆z = lcmp[3], ∆lcmp[3]
#         indz, wtz = distweights(src.c[3], gt_cmp[3], t_ind(bounds,3,3), z, ∆z, isbloch[3])
#
#         wt_vxl = SArray{Tuple{2,1,1}}(wtx) .* SArray{Tuple{1,2,1}}(wty) .* SArray{Tuple{1,1,2}}(wtz)  # 2×2×2 SArray
#         Jw = Iw∆w .* wt_vxl
#
#         for k = 1:2, j = 1:2, i = 1:2
#             j3d[indx[i],indy[j],indz[k],nw] += Jw[i,j,k]
#         end
#     end
#
#     return nothing
# end
