text:
"""
The vacuum speed of light is a global constant (299792458 meters per second).
"""
const c = 299792458;

## The first set of functions in this file are baseline functions that
##  modify (or are required by those that modify) a matrix representing
##  the Green's function, operating at complex frequency

"""

    GFVACGG!(GvacGG::AbstractArray{Complex{FT}, 2},
             startidx1::Integer, startidx2::Integer,
             freq::Union{FT, Complex{FT}}, blochk::Array{FT, 1},
             posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
             latticevecs::Array{FT, 2}, reciprocalvecs::Array{FT, 2},
             alphafield::FT, alphasource::FT,
             numdims::Integer) where FT<:AbstractFloat

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq` and wavevector `blochk` (a
3-element vector), corresponding to source position `possource` and field position
`posfield` (both 3-element vectors) described by Gaussian basis functions parameterized by
real atomic polarizabilities `alphasource` and `alphafield` respectively, using the
appropriate function for `numdims` periodic dimensions (with lattice parameters
`latticevecs` for the real lattice vectors and `reciprocalvecs` for the reciprocal lattice
vectors as appropriate). `numdims` must be 0, 1, or 2.

"""
function GFVACGG!(GvacGG::AbstractArray{Complex{FT}, 2},
                  startidx1::Integer, startidx2::Integer,
                  freq::Union{FT, Complex{FT}}, blochk::Array{FT, 1},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  latticevecs::Array{FT, 2}, reciprocalvecs::Array{FT, 2},
                  alphafield::FT, alphasource::FT,
                  numdims::Integer) where FT<:AbstractFloat

    if (numdims == 0)
        GFVACGG0!(GvacGG, startidx1, startidx2, freq,
                  posfield, possource, alphafield, alphasource);
    elseif (numdims == 1)
        GFVACGG1Ewald!(GvacGG, startidx1, startidx2, freq, blochk,
                       posfield, possource, latticevecs, reciprocalvecs,
                       alphafield, alphasource);
    elseif (numdims == 2)
        GFVACGG2Ewald!(GvacGG, startidx1, startidx2, freq, blochk,
                       posfield, possource, latticevecs, reciprocalvecs,
                       alphafield, alphasource);
    else
        error("GFVACGG!(): Periodic dimensionality numdims (currently ", numdims,
              ") must be 0 (no periodicity), 1, or 2 for GFVACGG");
    end

end

"""

    GFVACGG_sca!(GvacGG::AbstractArray{Complex{FT}, 2},
                 startidx1::Integer, startidx2::Integer,
                 freq::Union{FT, Complex{FT}}, blochk::Array{FT, 1},
                 posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                 latticevecs::Array{FT, 2}, reciprocalvecs::Array{FT, 2},
                 alphafield::FT, alphasource::FT,
                 numdims::Integer) where FT <: AbstractFloat

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with a 3-by-3 block of
zeros, which is the scattering Green's function in vacuum by definition, irrespective of the
other function arguments.

"""
function GFVACGG_sca!(GvacGG::AbstractArray{Complex{FT}, 2},
                      startidx1::Integer, startidx2::Integer,
                      freq::Union{FT, Complex{FT}}, blochk::Array{FT, 1},
                      posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                      latticevecs::Array{FT, 2}, reciprocalvecs::Array{FT, 2},
                      alphafield::FT, alphasource::FT,
                      numdims::Integer) where FT <: AbstractFloat

    if (length(size(GvacGG)) != 2 || size(GvacGG, 1) < 3 || size(GvacGG, 2) < 3)
        error("GFVACGG_sca!(): GvacGG (current size: ",
              size(GvacGG), ") must be 3-by-3");
    end
    if (startidx1 < 1 || startidx1 > size(GvacGG, 1) - 2)
        error("GFVACGG_sca!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ", size(GvacGG, 1) - 2,
              " as GvacGG has ", size(GvacGG, 1), " rows");
    end
    if (startidx2 < 1 || startidx2 > size(GvacGG, 2) - 2)
        error("GFVACGG_sca!(): startidx2 (currently ", startidx2,
              ") must be between 1 and ", size(GvacGG, 2) - 2,
              " as GvacGG has ", size(GvacGG, 2), " columns");
    end
    GvacGGtemp = view(GvacGG, startidx1:startidx1+2, startidx2:startidx2+2);
    fill!(GvacGGtemp, zero(eltype(GvacGGtemp)));

end

"""

    GFVACGG0!(GvacGG::AbstractArray{Complex{FT}, 2},
              startidx1::Integer, startidx2::Integer,
              freq::Union{FT, Complex{FT}},
              posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
              alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Fill `GvacGG[startidx1:startidx1+2,startidx2: startidx2+2]` in-place with the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq` without spatial periodicity,
corresponding to source position `possource` and field position `posfield` (both 3-element
vectors) described by Gaussian basis functions parameterized by real atomic polarizabilities
`alphasource` and `alphafield` respectively.

"""
function GFVACGG0!(GvacGG::AbstractArray{Complex{FT}, 2},
                   startidx1::Integer, startidx2::Integer,
                   freq::Union{FT, Complex{FT}},
                   posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                   alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    if (length(size(GvacGG)) != 2 || size(GvacGG, 1) < 3 || size(GvacGG, 2) < 3)
        error("GFVACGG0!(): GvacGG (current size: ", size(GvacGG), ") must be 3-by-3");
    end
    if (startidx1 < 1 || startidx1 > size(GvacGG, 1) - 2)
        error("GFVACGG0!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ", size(GvacGG, 1) - 2,
              " as GvacGG has ", size(GvacGG, 1), " rows");
    end
    if (startidx2 < 1 || startidx2 > size(GvacGG, 2) - 2)
        error("GFVACGG0!(): startidx2 (currently ", startidx2,
              ") must be between 1 and ", size(GvacGG, 2) - 2,
              " as GvacGG has ", size(GvacGG, 2), " columns");
    end

  # sigma_p = (alpha_p / 3)^(1/3) / (2*sqrt(pi)), sigmatot = sqrt(2(sigma_p^2 + sigma_q^2))
    sigmatot = hypot(cbrt(alphafield), cbrt(alphasource)) / (cbrt(3) * sqrt(2*pi));

    myrhat = posfield .- possource;
    rho = norm(myrhat)/sigmatot;
    GvacGGtemp = view(GvacGG, startidx1:startidx1+2, startidx2:startidx2+2);
    fill!(GvacGGtemp, zero(eltype(GvacGGtemp)));
    if (rho > eps())

        q = sigmatot*freq/c;
        normalize!(myrhat);

        term1 = exp(-1*(q^2 / 4) + 1im*rho*q) * erfc(-1im*(q/2) - rho);
        if (isnan(term1)) # just in case
            term1 = zero(eltype(GvacGGtemp));
        end
# exp(-0.25*q^2 - 1im*q*rho)*erfc(-0.5im*q + rho), computed with better numerical stability given rho >= 0
        term2 = exp(-1*(q^2 / 4) - 1im*q*rho - (-1im*(q/2) +
                                                rho)^2 + log(erfcx(-1im*(q/2) + rho)));
        if (isnan(term2)) # just in case
            term2 = zero(eltype(GvacGGtemp));
        end

        gvacGG = (term1 - term2)/(8*pi*sigmatot*rho);
        dgvacGGdrho = ((4/sqrt(pi))*rho*exp(-1*rho^2) +
                       (-1 + 1im*q*rho)*term1 +
                       (1 + 1im*q*rho)*term2)/(8*pi*sigmatot*rho^2);
        d2gvacGGdrho2 = ((-8/sqrt(pi))*(rho + rho^3)*exp(-1*rho^2) -
                         (-2 + 2im*q*rho + (q*rho)^2)*term1 +
                         (-2 - 2im*q*rho + (q*rho)^2)*term2)/(8*pi*sigmatot*rho^3);

        GvacGGtemp .+= (d2gvacGGdrho2/(sigmatot^2) -
                        dgvacGGdrho/(sigmatot^2 * rho)).*complex(myrhat*transpose(myrhat)) .+
        Matrix{Complex{FT}}((dgvacGGdrho/(sigmatot^2 * rho) + (freq/c)^2 * gvacGG)*I, 3, 3);

        if (sum(isnan.(GvacGGtemp)) > 0)
            fill!(GvacGGtemp, zero(eltype(GvacGGtemp)));
        end

    end

end

"""

    GFVACGG0_sca!(GvacGG::AbstractArray{Complex{FT}, 2},
                  startidx1::Integer, startidx2::Integer,
                  freq::Union{FT, Complex{FT}},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  alphafield::FT, alphasource::FT) where FT <: AbstractFloat

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with a 3-by-3 block of
zeros, which is the scattering Green's function in vacuum by definition, irrespective of the
other function arguments.

"""
function GFVACGG0_sca!(GvacGG::AbstractArray{Complex{FT}, 2},
                       startidx1::Integer, startidx2::Integer,
                       freq::Union{FT, Complex{FT}},
                       posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                       alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    if (length(size(GvacGG)) != 2 || size(GvacGG, 1) < 3 || size(GvacGG, 2) < 3)
        error("GFVACGG0_sca!(): GvacGG (current size: ", size(GvacGG), ") must be 3-by-3");
    end
    if (startidx1 < 1 || startidx1 > size(GvacGG, 1) - 2)
        error("GFVACGG0_sca!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ", size(GvacGG, 1) - 2,
              " as GvacGG has ", size(GvacGG, 1), " rows");
    end
    if (startidx2 < 1 || startidx2 > size(GvacGG, 2) - 2)
        error("GFVACGG: startidx2 (currently ", startidx2,
              ") must be between 1 and ", size(GvacGG, 2) - 2,
              " as GvacGG has ", size(GvacGG, 2), " columns");
    end
    GvacGGtemp = view(GvacGG, startidx1:startidx1+2, startidx2:startidx2+2);
    fill!(GvacGGtemp, zero(eltype(GvacGGtemp)));

end

"""

    GFVACGGcoincident!(GvacGG::AbstractArray{Complex{FT}, 2},
                       startidx1::Integer, startidx2::Integer,
                       freq::Union{FT, Complex{FT}}, alphafield::FT,
                       alphasource::FT) where FT<:AbstractFloat

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 tensor
block corresponding to vacuum Green's function interaction between two Gaussian basis
functions, defined by respective atomic polarizabilities `alphasource` and `alphafield`,
when their centers coincide; this quantity always has a finite imaginary part, and has a
finite real part when `alphafield` and `alphasource` are not both zero (meaning at least one
basis function has a finite width).

"""
function GFVACGGcoincident!(GvacGG::AbstractArray{Complex{FT}, 2},
                            startidx1::Integer, startidx2::Integer,
                            freq::Union{FT, Complex{FT}}, alphafield::FT,
                            alphasource::FT) where FT<:AbstractFloat

    if (length(size(GvacGG)) != 2 || size(GvacGG, 1) < 3 || size(GvacGG, 2) < 3)
        error("GFVACGGcoincident!(): GvacGG (current size: ",
              size(GvacGG), ") must be 3-by-3");
    end
    if (startidx1 < 1 || startidx1 > size(GvacGG, 1) - 2)
        error("GFVACGGcoincident!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ", size(GvacGG, 1) - 2,
              " as GvacGG has ", size(GvacGG, 1), " rows");
    end
    if (startidx2 < 1 || startidx2 > size(GvacGG, 2) - 2)
        error("GFVACGGcoincident!(): startidx2 (currently ", startidx2,
              ") must be between 1 and ", size(GvacGG, 2) - 2,
              " as GvacGG has ", size(GvacGG, 2), " columns");
    end

    setindex!(GvacGG,
              Matrix{eltype(GvacGG)}(GFVACGGcoincident(freq, alphafield, alphasource)*I,
                                     3, 3), startidx1:startidx1+2, startidx2:startidx2+2);

end

"""

    GFVACGG1Ewald!(GvacGG::AbstractArray{Complex{FT}, 2},
                   startidx1::Integer, startidx2::Integer,
                   freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                   posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                   latticevecs::AbstractArray{FT, 2},
                   reciprocalvecs::AbstractArray{FT, 2},
                   alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq` and wavevector `blochk` (a
3-element vector), corresponding to source position `possource` and field position
`posfield` (both 3-element vectors) described by Gaussian basis functions parameterized by
real atomic polarizabilities `alphasource` and `alphafield` respectively, in 1 periodic
dimension (with lattice parameters `latticevecs` for the real lattice vectors and
`reciprocalvecs` for the reciprocal lattice vectors as appropriate). This is done via Ewald
summation over Gaussian basis functions in 1 dimension.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFVACGG1Ewald!(GvacGG::AbstractArray{Complex{FT}, 2},
                        startidx1::Integer, startidx2::Integer,
                        freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                        posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                        latticevecs::AbstractArray{FT, 2},
                        reciprocalvecs::AbstractArray{FT, 2},
                        alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    if (length(size(GvacGG)) != 2 || size(GvacGG, 1) < 3 || size(GvacGG, 2) < 3)
        error("GFVACGG1Ewald!(): GvacGG (current size: ", size(GvacGG), ") must be 3-by-3");
    end
    if (startidx1 < 1 || startidx1 > size(GvacGG, 1) - 2)
        error("GFVACGG1Ewald!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ", size(GvacGG, 1) - 2,
              " as GvacGG has ", size(GvacGG, 1), " rows");
    end
    if (startidx2 < 1 || startidx2 > size(GvacGG, 2) - 2)
        error("GFVACGG1Ewald!(): startidx2 (currently ", startidx2,
              ") must be between 1 and ", size(GvacGG, 2) - 2,
              " as GvacGG has ", size(GvacGG, 2), " columns");
    end

    sigmatot = hypot(cbrt(alphafield), cbrt(alphasource)) / (cbrt(3) * sqrt(2*pi));

    myrvec = posfield .- possource;
    mylatticeunitvec = normalize(latticevecs[:,1]);
    myrperp = norm(myrvec .- mylatticeunitvec .* dot(mylatticeunitvec, myrvec));
    kappa = min(sqrt(pi)/norm(latticevecs[:,1]), 1/myrperp);
    N1max = round(Int, max(10, 10*myrperp/sigmatot));
    N2max = 0;
    Mmax = 10;
    Smax = Mmax;
    GvacGGtemp = view(GvacGG, startidx1:startidx1+2, startidx2:startidx2+2);
    fill!(GvacGGtemp, zero(Complex{FT}));

    if (norm(myrvec)/sigmatot <= eps())
        GvacGGtemp .+= GFVACGGEwaldSR(freq, blochk, posfield, possource,
                                      latticevecs, alphafield, alphasource,
                                      kappa, N1max, N2max) .+
        GFVACGGEwaldLR1D(freq, blochk, posfield, possource, reciprocalvecs[:,1],
                         alphafield, alphasource, kappa, Mmax, Smax) -
        GFVACGGcoincident(freq, alphafield, alphasource)*I;
    else
        GvacGGtemp .+= GFVACGGEwaldSR(freq, blochk, posfield, possource,
                                      latticevecs, alphafield, alphasource,
                                      kappa, N1max, N2max) .+
        GFVACGGEwaldLR1D(freq, blochk, posfield, possource, reciprocalvecs[:,1],
                         alphafield, alphasource, kappa, Mmax, Smax);
    end

end

"""

    GFVACGG2Ewald!(GvacGG::AbstractArray{Complex{FT}, 2},
                   startidx1::Integer, startidx2::Integer,
                   freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                   posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                   latticevecs::AbstractArray{FT, 2},
                   reciprocalvecs::AbstractArray{FT, 2},
                   alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq` and wavevector `blochk` (a
3-element vector), corresponding to source position `possource` and field position
`posfield` (both 3-element vectors) described by Gaussian basis functions parameterized by
real atomic polarizabilities `alphasource` and `alphafield` respectively, in 2 periodic
dimensions (with lattice parameters `latticevecs` for the real lattice vectors and
`reciprocalvecs` for the reciprocal lattice vectors as appropriate). This is done via Ewald
summation over Gaussian basis functions in 2 dimensions.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFVACGG2Ewald!(GvacGG::AbstractArray{Complex{FT}, 2},
                        startidx1::Integer, startidx2::Integer,
                        freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                        posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                        latticevecs::AbstractArray{FT, 2},
                        reciprocalvecs::AbstractArray{FT, 2},
                        alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    if (length(size(GvacGG)) != 2 || size(GvacGG, 1) < 3 || size(GvacGG, 2) < 3)
        error("GFVACGG2Ewald!(): GvacGG (current size: ", size(GvacGG), ") must be 3-by-3");
    end
    if (startidx1 < 1 || startidx1 > size(GvacGG, 1) - 2)
        error("GFVACGG2Ewald!(): startidx1 (currently ", startidx1,
              ") must be between 1 and ", size(GvacGG, 1) - 2,
              " as GvacGG has ", size(GvacGG, 1), " rows");
    end
    if (startidx2 < 1 || startidx2 > size(GvacGG, 2) - 2)
        error("GFVACGG2Ewald!(): startidx2 (currently ", startidx2,
              ") must be between 1 and ", size(GvacGG, 2) - 2,
              " as GvacGG has ", size(GvacGG, 2), " columns");
    end

    Auc = norm(cross(latticevecs[:,1], latticevecs[:,2]));
    M1max = 5;
    M2max = 5;
    kappa = 1e16;
    GvacGGtemp = view(GvacGG, startidx1:startidx1+2, startidx2:startidx2+2);
    fill!(GvacGGtemp, zero(Complex{FT}));
    GvacGGSR = zeros(Complex{FT}, 3, 3);
    sigmatot = hypot(cbrt(alphafield), cbrt(alphasource)) / (cbrt(3) * sqrt(2*pi));
    if (sigmatot < sqrt(Auc))
        kappa = sqrt(pi/Auc);
        N1max = 5;
        N2max = 5;
        setindex!(GvacGGSR, GFVACGGEwaldSR(freq, blochk, posfield, possource,
                                           latticevecs, alphafield, alphasource,
                                           kappa, N1max, N2max), 1:3, 1:3);
    end
    if (norm(posfield .- possource)/sigmatot <= eps())
        GvacGGtemp .+= GvacGGSR .+
        GFVACGGEwaldLR2D(freq, blochk, posfield, possource, reciprocalvecs, alphafield,
                         alphasource, kappa, M1max, M2max) -
        GFVACGGcoincident(freq, alphafield, alphasource)*I;
    else
        GvacGGtemp .+= GvacGGSR .+
        GFVACGGEwaldLR2D(freq, blochk, posfield, possource, reciprocalvecs, alphafield,
                         alphasource, kappa, M1max, M2max);
    end

end

"""

    GFVACGGEwaldSR(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                   posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                   latticevecs::AbstractArray{FT, 2},
                   alphafield::FT, alphasource::FT, kappa::FT,
                   N1max::Integer, N2max::Integer=0) where FT<:AbstractFloat

Return the short-range Ewald summation contribution to the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq` and wavevector `blochk` (a
3-element vector), corresponding to source position `possource` and field position
`posfield` (both 3-element vectors) described by Gaussian basis functions parameterized by
real atomic polarizabilities `alphasource` and `alphafield` respectively, in 1 or 2 periodic
dimensions (with lattice parameters `latticevecs` for the real lattice vectors as
appropriate). The sums over the real lattice run over `-1*N1max:N1max` and `-1*N2max:N2max`
along each lattice vector, and `kappa` is the Ewald cutoff parameter.

"""
function GFVACGGEwaldSR(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                        posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                        latticevecs::AbstractArray{FT, 2},
                        alphafield::FT, alphasource::FT, kappa::FT,
                        N1max::Integer, N2max::Integer=0) where FT<:AbstractFloat

    sigmatot = hypot(cbrt(alphafield), cbrt(alphasource)) / (cbrt(3) * sqrt(2*pi));

    q = sigmatot*freq/c;
    mynu = 1/hypot(1, 1/(sigmatot*kappa));

    n1list = -1*N1max:N1max;
    n2list = -1*N2max:N2max;

    GvacGGSR = zeros(Complex{FT}, 3, 3);
    currGvacGGSR = zeros(Complex{FT}, 3, 3);

    for n2 in n2list, n1 in n1list

        currlatticevec = n1*latticevecs[:,1] .+ n2*latticevecs[:,2];
        kdotR = dot(blochk, currlatticevec);
        myrhat = posfield .+ currlatticevec .- possource;
        rho = norm(myrhat)/sigmatot;
        nurho = mynu*rho;

        # simple expression for atoms coinciding in the same unit cell
        if (rho <= eps())

            setindex!(currGvacGGSR,
                      Matrix{eltype(currGvacGGSR)}((((2/sqrt(pi))*(q^2 - 1 - mynu*(q^2 - mynu^2)*exp((q^2 / 4) * (1/mynu^2 - 1))) + 1im*q^3 * (erf(1im*q/2) - erf(1im*q/(2*mynu))))/(6*pi*sigmatot^3))*I, 3, 3),
                      1:3, 1:3);

    # more complicated otherwise
        else

            fill!(currGvacGGSR, zero(eltype(currGvacGGSR)));
            normalize!(myrhat);
            term1 = exp(-1*(q^2 / 4) + 1im*q*rho)*erfc(-1im*(q/2) - rho);#exp(-0.25*q^2 + 1im*q*rho - (-0.5im*q - rho)^2 + log(erfcx(-0.5im*q - rho)));
            if (isnan(term1)) # just in case
                term1 = zero(eltype(currGvacGGSR));
            end
            term2 = exp(-1*(q^2 / 4) + 1im*q*rho)*erfc(-1im*q/(2*mynu) - nurho);#exp(-0.25*q^2 + 1im*q*rho - (-0.5im*q/mynu - nurho)^2 + log(erfcx(-0.5im*q/mynu - nurho)));
            if (isnan(term2)) # just in case
                term2 = zero(eltype(currGvacGGSR));
            end
# exp(-0.25*q^2 - 1im*q*rho)*erfc(-0.5im*q + rho), computed with better numerical stability given rho >= 0
            term3 = exp(-1*(q^2 / 4) - 1im*q*rho - (-1im*(q/2) + rho)^2 +
                        log(erfcx(-1im*(q/2) + rho)));
            if (isnan(term3)) # just in case
                term3 = zero(eltype(currGvacGGSR));
            end
# exp(-0.25*q^2 - 1im*q*rho)*erfc(-0.5im*q/mynu + nurho), computed with better numerical stability given rho >= 0
            term4 = exp(-1*(q^2 / 4) - 1im*q*rho - (-1im*q/(2*mynu) + nurho)^2 +
                        log(erfcx(-1im*q/(2*mynu) + nurho)));
            if (isnan(term4)) # just in case
                term4 = zero(eltype(currGvacGGSR));
            end

            gvacGGSR = ((term1 - term2) - (term3 - term4))/(8*pi*sigmatot*rho);
            if (isnan(gvacGGSR)) # just in case
              gvacGGSR = zero(gvacGGSR);
            end
            dgvacGGSRdrho = ((4/sqrt(pi))*(exp(-1*rho^2) * rho -
                                           exp((q^2 / 4) * (1/mynu^2 - 1) - (nurho)^2) *
                                           nurho) +
                             (-1 + 1im*q*rho)*(term1 - term2) +
                             (1 + 1im*q*rho)*(term3 - term4))/(8*pi*sigmatot*rho^2);
            if (isnan(dgvacGGSRdrho)) # just in case
                dgvacGGSRdrho = zero(dgvacGGSRdrho);
            end
            d2gvacGGSRdrho2 = ((8/sqrt(pi))*(exp((q^2 / 4) * (1/mynu^2 - 1) -
                                                 nurho^2) * (nurho + nurho^3) -
                                             exp(-1*rho^2) * (rho + rho^3)) -
                               (-2 + 2im*q*rho + (q*rho)^2)*(term1 - term2) +
                               (-2 - 2im*q*rho +
                                (q*rho)^2)*(term3 - term4))/(8*pi*sigmatot*rho^3);
            if (isnan(d2gvacGGSRdrho2)) # just in case
                d2gvacGGSRdrho2 = zero(d2gvacGGSRdrho2);
            end

            currGvacGGSR .+= exp(-1im*kdotR) .*
            (((d2gvacGGSRdrho2/(sigmatot^2) - dgvacGGSRdrho/(sigmatot^2 * rho))) .*
             complex(myrhat * transpose(myrhat)) .+
             Matrix{eltype(currGvacGGSR)}(((dgvacGGSRdrho/(sigmatot^2 * rho) +
                                            (freq/c)^2 * gvacGGSR))*I, 3, 3));

            if (sum(isnan.(currGvacGGSR)) > 0) # just in case
                fill!(currGvacGGSR, zero(eltype(currGvacGGSR)));
            end
        end

        GvacGGSR .+= currGvacGGSR;

    end

    return GvacGGSR;

end

"""

    GFVACGGEwaldLR1D(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                     posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                     reciprocalvec::AbstractArray{FT, 1},
                     alphafield::FT, alphasource::FT, kappa::FT,
                     Mmax::Integer, Smax::Integer) where FT<:AbstractFloat

Return the long-range Ewald summation contribution to the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq` and wavevector `blochk` (a
3-element vector), corresponding to source position `possource` and field position
`posfield` (both 3-element vectors) described by Gaussian basis functions parameterized by
real atomic polarizabilities `alphasource` and `alphafield` respectively, in 1 periodic
dimension (with lattice parameters `reciprocalvec` for the reciprocal lattice vector as
appropriate). The sum over the reciprocal lattice runs over `-1*Mmax:Mmax`, the sum
computing the series representation of the reciprocal lattice expression runs over
`0:Smax`, and `kappa` is the Ewald cutoff parameter.

"""
function GFVACGGEwaldLR1D(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                          posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                          reciprocalvec::AbstractArray{FT, 1},
                          alphafield::FT, alphasource::FT, kappa::FT,
                          Mmax::Integer, Smax::Integer) where FT<:AbstractFloat

    if (real(freq) != 0.0 && imag(freq) != 0.0)
        error(string("GFVACGGEwaldLR1D(): freq is", freq,
                     "; GFVACGG for 1D-periodic systems can only handle freq purely",
                     " along the real or imaginary axes"));
    end

    sigmatot = hypot(cbrt(alphafield), cbrt(alphasource)) / (cbrt(3) * sqrt(2*pi));

    q = sigmatot*freq/c;

    mynu2 = 1/(1 + 1/(sigmatot*kappa)^2);
    GvacGGLR = zeros(Complex{FT}, 3, 3);

# extract component of atomic displacement radial (perpendicular) with respect to array axis
    myreciprocalunitvec = normalize(reciprocalvec);
    myrvec = posfield .- possource;
    projperpreciprocal = myreciprocalunitvec * transpose(myreciprocalunitvec);
    subtractfromdiag!(projperpreciprocal, ones(FT, 3));
    myrhoperpvec = projperpreciprocal * myrvec/sigmatot;
    myrhoperp = norm(myrhoperpvec);
    nu2rhoperp2 = mynu2*myrhoperp^2;

    for mm=-1*Mmax:Mmax

        kgvec = blochk .+ reciprocalvec.*mm;
        h2over4nu2 = real(sigmatot^2 * dot(kgvec, kgvec) - q^2)/(4*mynu2);

        freqckg = Matrix{Complex{FT}}((freq/c)^2 * I, 3, 3) .-
        complex(kgvec * transpose(kgvec));

# if the source and field points lie along the same line, evaluate a simpler expression
        if (myrhoperp <= eps())

            currE1 = E1func(h2over4nu2);
            currE2 = exp(-1*h2over4nu2) - h2over4nu2*currE1;
            GvacGGLR .+= exp(-1*(q^2 / 4) + 1im*dot(kgvec, myrvec)) .*
            (freqckg .* currE1 .- (2*mynu2/sigmatot^2).*projperpreciprocal.*currE2);

        else # do the full sum over exponential integral orders

            sum1 = zero(FT);
            sum2 = zero(FT);
            sum3 = zero(FT);
            prevEint = E1func(h2over4nu2);
            currEint = copy(prevEint);
            coeff1 = zero(FT);
            coeff2 = zero(FT);
            coeff3 = zero(FT);

            for ss=0:Smax

                if (ss == 0)
                    coeff1 = one(coeff1);
                else

                    coeff1 = (-1/ss)*nu2rhoperp2 * coeff1;
                    currEint = (1/ss)*(exp(-1*h2over4nu2) - h2over4nu2*prevEint);
                    prevEint = copy(currEint);

                    if (ss == 1)
                        coeff2 = -2.0*mynu2;
                    else

                        coeff2 = (-1/(ss-1))*nu2rhoperp2 * coeff2;
                        if (ss == 2)
                            coeff3 = 4*mynu2^2;
                        else
                            coeff3 = (-1/(ss-2))*nu2rhoperp2 * coeff3;
                        end

                    end

                end

                sum1 += coeff1*currEint;
                sum2 += coeff2*currEint;
                sum3 += coeff3*currEint;

            end

            GvacGGLR .+= exp(-1*(q^2 / 4) + 1im*dot(kgvec, myrvec)) .*
            (sum1.*freqckg .+ sum2.*(projperpreciprocal./(sigmatot^2) .+
                                     (2im/sigmatot) .*
                                     sympart(complex(kgvec*transpose(myrhoperpvec)))) .+
             sum3.*complex(myrhoperpvec * transpose(myrhoperpvec))./(sigmatot^2));

        end

    end

    lmul!(norm(reciprocalvec)/(8*pi^2), GvacGGLR);
    return GvacGGLR;

end

"""

    GFVACGGEwaldLR2D(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                     posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                     reciprocalvecs::AbstractArray{FT, 2},
                     alphafield::FT, alphasource::FT, kappa::FT,
                     M1max::Integer, M2max::Integer) where FT<:AbstractFloat

Return the long-range Ewald summation contribution to the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq` and wavevector `blochk` (a
3-element vector), corresponding to source position `possource` and field position
`posfield` (both 3-element vectors) described by Gaussian basis functions parameterized by
real atomic polarizabilities `alphasource` and `alphafield` respectively, in 2 periodic
dimensions (with lattice parameters `reciprocalvecs` for the reciprocal lattice vectors as
appropriate). The sum over the reciprocal lattice runs over `-1*M1max:M1max` and
`-1*M2max:M2max` along each reciprocal lattice vector, and `kappa` is the Ewald cutoff
parameter.

"""
function GFVACGGEwaldLR2D(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                          posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                          reciprocalvecs::AbstractArray{FT, 2},
                          alphafield::FT, alphasource::FT, kappa::FT,
                          M1max::Integer, M2max::Integer) where FT<:AbstractFloat

    sigmatot = hypot(cbrt(alphafield), cbrt(alphasource)) / (cbrt(3) * sqrt(2*pi));

    q = sigmatot*freq/c;

    mymu = sigmatot*kappa;
    mynu = mymu/hypot(1, mymu);
    mynu2 = 1/(1 + 1/mymu^2);
    GvacGGLR = zeros(Complex{FT}, 3, 3);

    ## unit vector normal to plane, and prefactor involving Brillouin zone area
    perpunitvec = cross(reciprocalvecs[:,1], reciprocalvecs[:,2]);
    ABZover16pi2 = norm(perpunitvec) / (16*pi^2);
    normalize!(perpunitvec);

    ## component of atomic displacement normal to the plane
    myrvec = posfield .- possource;
    myrhoperpsigned = dot(perpunitvec, myrvec)/sigmatot;

    for m2=-1*M2max:M2max, m1=-1*M1max:M1max

        kgvec = blochk .+ m1.*reciprocalvecs[:,1] .+ m2.*reciprocalvecs[:,2];
        kgnorm2 = dot(kgvec, kgvec);
        kgdotr = dot(kgvec, myrvec);

        freqckg = Matrix{Complex{FT}}((freq/c)^2 * I, 3, 3) .-
        complex(kgvec * transpose(kgvec));

        h = sqrt(complex(sigmatot^2 * kgnorm2 - q^2));

        ## if the source and field points do lie parallel to the lattice plane
        f = 2*exp(-1*(q^2 / 4) + 1im*kgdotr)*erfc(h/(2*mynu));
        dfdrhoperp = zero(f);
        d2fdrhoperp2 = h^2 * f -
        (4*mynu*h/sqrt(pi))*exp((q^2 / (4*mymu^2)) +
                                1im*kgdotr - sigmatot^2 * kgnorm2 / (4*mynu^2));

# if the source and field points do not lie parallel to the lattice plane
        if (abs(myrhoperpsigned) > eps())

# exp(-0.25*q^2 + 1im*kgdotr + h*myrhoperpsigned) * erfc(h/(2*mynu) + mynu*myrhoperpsigned), computed with possibly better numerical stability
            term1 = exp(-1*(q^2 / 4) +
                        1im*kgdotr + h*myrhoperpsigned -
                        (h/(2*mynu) + mynu*myrhoperpsigned)^2 +
                        log(erfcx(h/(2*mynu) + mynu*myrhoperpsigned)));
# no need for that in exponential decay
            term2 = exp(-1*(q^2 / 4) + 1im*kgdotr - h*myrhoperpsigned) *
            erfc(h/(2*mynu) - mynu*myrhoperpsigned);
#      f = exp(-0.25*q^2 + 1im*kgdotr + h*myrhoperpsigned) * erfc(h/(2*mynu) + mynu*myrhoperpsigned) + exp(-0.25*q^2 + 1im*kgdotr - h*myrhoperpsigned) * erfc(h/(2*mynu) - mynu*myrhoperpsigned);
            f = term1 + term2;
#      dfdrhoperp = h*(exp(-0.25*q^2 + 1im*kgdotr + h*myrhoperpsigned) * erfc(h/(2*mynu) - mynu*myrhoperpsigned) + exp(-0.25*q^2 + 1im*kgdotr - h*myrhoperpsigned) * erfc(h/(2*mynu) - mynu*myrhoperpsigned));
            dfdrhoperp = h*(term1 - term2);

            d2fdrhoperp2 = h^2 * f -
            (4*mynu*h/sqrt(pi))*
            exp((q^2 / (4*mymu^2)) + 1im*kgdotr -
                sigmatot^2 * kgnorm2/(4*mynu^2) - (mynu*myrhoperpsigned)^2);

        end

        GvacGGLR .+= (1/h) .* (sigmatot .* freqckg .* f .+
                               2im .* dfdrhoperp .*
                               sympart(complex(kgvec * transpose(perpunitvec))) .+
                               (1/sigmatot) .* d2fdrhoperp2 .*
                               complex(perpunitvec * transpose(perpunitvec)));

    end

    lmul!(ABZover16pi2, GvacGGLR);
    return GvacGGLR;

end

## The second set of functions in this file consists of aliases for
##  functions with fewer arguments: in particular, GFVACGG!(GvacGG,
##  startidx1, startidx2, freq, posfield, possource, alphafield,
##  alphasource) is the same as GFVACGG0!(GvacGG, startidx1,
##  startidx2, freq, posfield, possource, alphafield, alphasource),
##  which is convenient for compact molecules

"""

    GFVACGG!(GvacGG::AbstractArray{Complex{FT}, 2}, startidx1::Integer, startidx2::Integer,
             freq::Union{FT, Complex{FT}},
             posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
             alphafield::FT, alphasource::FT) where {FT<:AbstractFloat}

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 Green's
function interaction tensor in vacuum at frequency `freq`, corresponding to source position
`possource` and field position `posfield` (both 3-element vectors) described by Gaussian
basis functions parameterized by real atomic polarizabilities `alphasource` and `alphafield`
respectively, for nonperiodic geometries.

"""
GFVACGG!(GvacGG::AbstractArray{Complex{FT}, 2}, startidx1::Integer, startidx2::Integer,
         freq::Union{FT, Complex{FT}},
         posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
         alphafield::FT, alphasource::FT) where {FT<:AbstractFloat} =
GFVACGG0!(GvacGG, startidx1, startidx2, freq, posfield, possource, alphafield, alphasource);

"""

    GFVACGG_sca!(GvacGG::AbstractArray{Complex{FT}, 2}, startidx1::Integer, startidx2::Integer,
                 freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
                 possource::AbstractArray{FT, 1}, alphafield::FT,
                 alphasource::FT) where {FT<:AbstractFloat}

Fill `GvacGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with a 3-by-3 block of
zeros, which is the scattering Green's function in vacuum by definition, irrespective of the
other function arguments.

"""
GFVACGG_sca!(GvacGG::AbstractArray{Complex{FT}, 2}, startidx1::Integer, startidx2::Integer,
             freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
             possource::AbstractArray{FT, 1}, alphafield::FT,
             alphasource::FT) where {FT<:AbstractFloat} =
GFVACGG0_sca!(GvacGG, startidx1, startidx2, freq, posfield, possource, alphafield,
              alphasource);

## The third set of functions in this file consists of aliases for not
##  in-place function versions of the first and second sets

"""

    GFVACGG(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
            posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
            latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
            alphafield::FT, alphasource::FT, numdims::Integer) where FT<:AbstractFloat

Return the 3-by-3 Green's function interaction tensor in vacuum at frequency `freq` and
wavevector `blochk` (a 3-element vector), corresponding to source position `possource` and
field position `posfield` (both 3-element vectors) described by Gaussian basis functions
parameterized by real atomic polarizabilities `alphasource` and `alphafield` respectively,
using the appropriate function for `numdims` periodic dimensions (with lattice parameters
`latticevecs` for the real lattice vectors and `reciprocalvecs` for the reciprocal lattice
vectors as appropriate). `numdims` must be 0, 1, or 2.

"""
function GFVACGG(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                 posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                 latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
                 alphafield::FT, alphasource::FT, numdims::Integer) where FT<:AbstractFloat

    GvacGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFVACGG!(GvacGG, 1, 1, freq, blochk, posfield, possource,
             latticevecs, reciprocalvecs, alphafield, alphasource, numdims);
    return GvacGG;

end

"""

    GFVACGG_sca(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
                alphafield::FT, alphasource::FT,
                numdims::Integer) where {FT<:AbstractFloat}

Return a 3-by-3 block of zeros, which is the scattering Green's function in vacuum by
definition, irrespective of the other function arguments.

"""
GFVACGG_sca(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
            posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
            latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
            alphafield::FT, alphasource::FT,
            numdims::Integer) where {FT<:AbstractFloat} = zeros(Complex{FT}, 3, 3);

"""

    GFVACGG0(freq::Union{FT, Complex{FT}},
             posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
             alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Return the 3-by-3 Green's function interaction tensor in vacuum at frequency `freq` without
spatial periodicity, corresponding to source position `possource` and field position
`posfield` (both 3-element vectors) described by Gaussian basis functions parameterized by
real atomic polarizabilities `alphasource` and `alphafield` respectively.

"""
function GFVACGG0(freq::Union{FT, Complex{FT}},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    GvacGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFVACGG0!(GvacGG, 1, 1, freq, posfield, possource, alphafield, alphasource);
    return GvacGG;

end

"""

    GFVACGG0_sca(freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
                 possource::AbstractArray{FT, 1}, alphafield::FT,
                 alphasource::FT) where {FT<:AbstractFloat}

Return a 3-by-3 block of zeros, which is the scattering Green's function in vacuum by
definition, irrespective of the other function arguments.

"""
GFVACGG0_sca(freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
             possource::AbstractArray{FT, 1}, alphafield::FT,
             alphasource::FT) where {FT<:AbstractFloat} = zeros(Complex{FT}, 3, 3);

"""

    GFVACGGcoincident(freq::Union{FT, Complex{FT}},
                      alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Return the scalar (as the tensor is the scalar multiplied by the 3-by-3 identity)
corresponding to vacuum Green's function interaction between two Gaussian basis functions,
defined by respective atomic polarizabilities `alphasource` and `alphafield`, when their
centers coincide; this quantity always has a finite imaginary part, and has a finite real
part when `alphafield` and `alphasource` are not both zero (meaning at least one basis
function has a finite width).

"""
function GFVACGGcoincident(freq::Union{FT, Complex{FT}},
                           alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    sigmatot = hypot(cbrt(alphafield), cbrt(alphasource)) / (cbrt(3) * sqrt(2*pi));
    q = sigmatot*freq/c;
    return (((2/sqrt(pi))*(q^2 - 1) + 1im*q^3 * erfc(-1im*q/2))/(6*pi*sigmatot^3));

end

"""

    GFVACGG1Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  latticevecs::AbstractArray{FT, 2},
                  reciprocalvecs::AbstractArray{FT, 2},
                  alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Return the 3-by-3 Green's function interaction tensor in vacuum at frequency `freq` and
wavevector `blochk` (a 3-element vector), corresponding to source position `possource` and
field position `posfield` (both 3-element vectors) described by Gaussian basis functions
parameterized by real atomic polarizabilities `alphasource` and `alphafield` respectively,
in 1 periodic dimension (with lattice parameters `latticevecs` for the real lattice vectors
and `reciprocalvecs` for the reciprocal lattice vectors as appropriate). This is done via
Ewald summation over Gaussian basis functions in 1 dimension.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFVACGG1Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                       posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                       latticevecs::AbstractArray{FT, 2},
                       reciprocalvecs::AbstractArray{FT, 2},
                       alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    GvacGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFVACGG1Ewald!(GvacGG, 1, 1, freq, blochk, posfield, possource,
                   latticevecs, reciprocalvecs, alphafield, alphasource);
    return GvacGG;

end

"""

    GFVACGG2Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  latticevecs::AbstractArray{FT, 2},
                  reciprocalvecs::AbstractArray{FT, 2}, alphafield::FT,
                  alphasource::FT) where FT<:AbstractFloat

Return the 3-by-3 Green's function interaction tensor in vacuum at frequency `freq` and
wavevector `blochk` (a 3-element vector), corresponding to source position `possource` and
field position `posfield` (both 3-element vectors) described by Gaussian basis functions
parameterized by real atomic polarizabilities `alphasource` and `alphafield` respectively,
in 2 periodic dimensions (with lattice parameters `latticevecs` for the real lattice vectors
and `reciprocalvecs` for the reciprocal lattice vectors as appropriate). This is done via
Ewald summation over Gaussian basis functions in 2 dimensions.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFVACGG2Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                       posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                       latticevecs::AbstractArray{FT, 2},
                       reciprocalvecs::AbstractArray{FT, 2}, alphafield::FT,
                       alphasource::FT) where FT<:AbstractFloat

    GvacGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFVACGG2Ewald!(GvacGG, 1, 1, freq, blochk, posfield, possource,
                   latticevecs, reciprocalvecs, alphafield, alphasource);
    return GvacGG;

end

"""

    GFVACGG(freq::Union{FT, Complex{FT}},
            posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
            alphafield::FT, alphasource::FT) where {FT<:AbstractFloat}

Return the 3-by-3 Green's function interaction tensor in vacuum at frequency `freq`,
corresponding to source position `possource` and field position `posfield` (both 3-element
vectors) described by Gaussian basis functions parameterized by real atomic polarizabilities
`alphasource` and `alphafield` respectively, for nonperiodic geometries.

"""
GFVACGG(freq::Union{FT, Complex{FT}},
        posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
        alphafield::FT, alphasource::FT) where {FT<:AbstractFloat} =
GFVACGG0(freq, posfield, possource, alphafield, alphasource);


text:
# Print the benchmark results as a markdown table

times = Dict{String, Dict}()

open("benchmarks.csv") do f
    for line in eachline(f)
        if !startswith(line, "Package")
            software, benchmark, runtime = split(line, ",")
            if haskey(times, benchmark)
                times[benchmark][software] = parse(Float64, runtime)
            else
                times[benchmark] = Dict(software=> parse(Float64, runtime))
            end
        end
    end
end

for (benchmark, label, millisecond) in (
                            ("Parse PDB 1CRN"  , "Parse PDB 1CRN / ms"  , true ),
                            ("Parse PDB 1HTQ"  , "Parse PDB 1HTQ / s"   , false),
                            ("Parse mmCIF 1CRN", "Parse mmCIF 1CRN / ms", true ),
                            ("Parse mmCIF 1HTQ", "Parse mmCIF 1HTQ / s" , false),
                            ("Parse MMTF 1CRN" , "Parse MMTF 1CRN / ms" , true ),
                            ("Parse MMTF 1HTQ" , "Parse MMTF 1HTQ / s"  , false),
                            ("Count"           , "Count / ms"           , true ),
                            ("Distance"        , "Distance / ms"        , true ),
                            ("Ramachandran"    , "Ramachandran / ms"    , true ))
    print("| $(rpad(label, 21)) |")
    for software in ("BioStructures", "MIToS", "Biopython", "ProDy", "MDAnalysis", "biotite",
                        "atomium", "Bio3D", "Rpdb", "BioJava", "BioPerl", "BioRuby", "GEMMI",
                        "Victor", "ESBTL", "chemfiles-python", "chemfiles-cxx")
        if haskey(times[benchmark], software)
            if millisecond
                val = string(round(1000 * times[benchmark][software], sigdigits=2))
            else
                val = string(round(times[benchmark][software], sigdigits=2))
            end
        else
            val = "-"
        end
        print(" $(rpad(val, 16)) |")
    end
    println()
end

ext: jl

max_stars_repo_name: wojdyr/pdb-benchmarks

max_issues_repo_name: wojdyr/pdb-benchmarks

max_forks_repo_name: wojdyr/pdb-benchmarks


text:
module Types

export BucketService,
    DaemonService,
    DatasourceService,
    QueueService,
    CacheService

abstract BucketService
abstract DatasourceService
abstract QueueService
abstract CacheService

export DaemonTaskDetails

"""
    DaemonTaskDetails

This is the base composite abstract class used to compose Details and Payload
i.e. compose a task with
```julia
type YourDaemonTaskDetails <: DaemonTaskDetails
    basic_info::BasicTask.Info
    payload_info::YourTask.Info
end
```
"""
abstract DaemonTaskDetails

end # module Types

max_forks_repo_name: seung-lab/SimpleTasks


text:
# ---
# title: 57. Insert Interval
# id: problem57
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sort
# link: <https://leetcode.com/problems/insert-interval/description/>
# hidden: true
# ---
#
# Given a set of _non-overlapping_ intervals, insert a new interval into the
# intervals (merge if necessary).
#
# You may assume that the intervals were initially sorted according to their
# start times.
#
#
#
# **Example 1:**
#
#
#
#     Input: intervals = [[1,3],[6,9]], newInterval = [2,5]
#     Output: [[1,5],[6,9]]
#
#
# **Example 2:**
#
#
#
#     Input: intervals = [[1,2],[3,5],[6,7],[8,10],[12,16]], newInterval = [4,8]
#     Output: [[1,2],[3,10],[12,16]]
#     Explanation: Because the new interval [4,8] overlaps with [3,5],[6,7],[8,10].
#
# **Example 3:**
#
#
#
#     Input: intervals = [], newInterval = [5,7]
#     Output: [[5,7]]
#
#
# **Example 4:**
#
#
#
#     Input: intervals = [[1,5]], newInterval = [2,3]
#     Output: [[1,5]]
#
#
# **Example 5:**
#
#
#
#     Input: intervals = [[1,5]], newInterval = [2,7]
#     Output: [[1,7]]
#
#
#
#
# **Constraints:**
#
#   * `0 <= intervals.length <= 104`
#   * `intervals[i].length == 2`
#   * `0 <= intervals[i][0] <= intervals[i][1] <= 105`
#   * `intervals` is sorted by `intervals[i][0]` in **ascending**  order.
#   * `newInterval.length == 2`
#   * `0 <= newInterval[0] <= newInterval[1] <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "@add_predicates" begin
    qc = Quickcheck("Set for @add_predicates tests")

    @test isa(@add_predicate(qc, "Identity", (x::Float64 -> x == x)),
              Quickcheck)
    @test isa(@add_predicate(qc, "Is odd", (n::Int -> isodd(n))),
              Quickcheck)
    @test isa(@add_predicate(qc, "Sum commute",
                             ((n::Int, x::Float64) -> n + x == x + n)),
              Quickcheck)

    @test_throws(ErrorException,
                 @add_predicate(qc, "Identity", x == x))

    @test_throws(ErrorException,
                 @add_predicate(qc, "Identity", x -> x == x))
    @test_throws(ErrorException,
                 @add_predicate(qc, "sum commute",
                                (n, x::Float64) -> n + x == x + n))
    @test_throws(ErrorException,
                 @add_predicate(qc, "sum commute",
                                (n::Int, x) -> n + x == x + n))
    @test_throws(ErrorException,
                 @add_predicate(qc, "sum commute",
                                (n, x) -> n + x == x + n))

    @test_throws(ErrorException,
                 @add_predicate(qc, "sum commute",
                                (n, x -> n + x == x + n)))
    @test_throws(ErrorException,
                 @add_predicate(qc, "sum commute",
                                [n, x] -> n + x == x + n))

    ## Test that every type for which a `generate` method exists can
    ## be added to a `Quickcheck`.
    for type ∈ types_with_generate
        @eval begin
            qc = Quickcheck("dummy")
            @test @add_predicate(qc,
                                 string($type) * " true",
                                 x::$type -> true) isa Quickcheck
            @test isnothing(@quickcheck qc)
        end
    end
end


text:
export GaussianWeightedMeanPrecision

"""
Description:

    A Gaussian with weighted-mean-precision parameterization:

    f(out,xi,w) = 𝒩(out|xi,w) = (2π)^{-D/2} |w|^{1/2} exp(-1/2 xi' w^{-1} xi) exp(-1/2 xi' w xi + out' xi)

Interfaces:

    1. out
    2. xi (weighted mean, w*m)
    3. w (precision)

Construction:

    GaussianWeightedMeanPrecision(out, xi, w, id=:some_id)
"""
mutable struct GaussianWeightedMeanPrecision <: Gaussian
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function GaussianWeightedMeanPrecision(out, xi, w; id=generateId(GaussianWeightedMeanPrecision))
        @ensureVariables(out, xi, w)
        self = new(id, Array{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:xi] = self.interfaces[2] = associate!(Interface(self), xi)
        self.i[:w] = self.interfaces[3] = associate!(Interface(self), w)

        return self
    end
end

slug(::Type{GaussianWeightedMeanPrecision}) = "𝒩"

format(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = "$(slug(GaussianWeightedMeanPrecision))(xi=$(format(dist.params[:xi])), w=$(format(dist.params[:w])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{GaussianWeightedMeanPrecision}; xi=0.0, w=1.0) = ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}(Dict(:xi=>xi, :w=>w))
ProbabilityDistribution(::Type{GaussianWeightedMeanPrecision}; xi::Number=0.0, w::Number=1.0) = ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}(Dict(:xi=>xi, :w=>w))
ProbabilityDistribution(::Type{Multivariate}, ::Type{GaussianWeightedMeanPrecision}; xi=[0.0], w=transpose([1.0])) = ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}(Dict(:xi=>xi, :w=>w))

dims(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = length(dist.params[:xi])

vague(::Type{GaussianWeightedMeanPrecision}) = ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=tiny)
vague(::Type{GaussianWeightedMeanPrecision}, dims::Int64) = ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(dims), w=tiny*diageye(dims))
vague(::Type{GaussianWeightedMeanPrecision}, dims::Tuple{Int64}) = ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=zeros(dims), w=tiny*diageye(dims[1]))

unsafeMean(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = cholinv(dist.params[:w])*dist.params[:xi]

unsafeMode(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = cholinv(dist.params[:w])*dist.params[:xi]

unsafeVar(dist::ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}) = 1.0/dist.params[:w] # unsafe variance
unsafeVar(dist::ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}) = diag(cholinv(dist.params[:w]))

unsafeCov(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = cholinv(dist.params[:w])

function unsafeMeanCov(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType
    v = cholinv(dist.params[:w])
    return (v*dist.params[:xi], v)
end

unsafeWeightedMean(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = deepcopy(dist.params[:xi]) # unsafe weighted mean

unsafePrecision(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = deepcopy(dist.params[:w]) # unsafe precision

unsafeWeightedMeanPrecision(dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType = (deepcopy(dist.params[:xi]), deepcopy(dist.params[:w]))

function logPdf(dist::ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}, x)
    m = dist.params[:xi]/dist.params[:w]
    return -0.5*(log(2pi) - log(dist.params[:w]) + (x-m)^2*dist.params[:w])
end

function logPdf(dist::ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}, x)
    m = cholinv(dist.params[:w])*dist.params[:xi]
    return -0.5*(dims(dist)*log(2pi) - logdet(dist.params[:w]) + transpose(x-m)*dist.params[:w]*(x-m))
end

isProper(dist::ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}) = (floatmin(Float64) < dist.params[:w] < floatmax(Float64))
isProper(dist::ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}) = isRoundedPosDef(dist.params[:w])

function ==(t::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, u::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType
    (t === u) && return true
    isApproxEqual(t.params[:xi], u.params[:xi]) && isApproxEqual(t.params[:w], u.params[:w])
end

function ==(t::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, u::ProbabilityDistribution{V, F}) where {V<:VariateType, F<:Gaussian}
    (t === u) && return true
    isApproxEqual(t.params[:xi], unsafeWeightedMean(u)) && isApproxEqual(t.params[:w], unsafePrecision(u))
end


text:
"""
    solve_adbc(equation, u, tlist, n, dt = (tlist[2] - tlist[1]) / 200)

Solve filtered equation using approximate deconvolution boundary conditions (ADBC). This
approach was proposed in [Borggaard 2006].

J. Borggaard, T. Iliescu, Approximate deconvolution boundary conditions for large eddy
simulation, Applied Mathematics Letters 19 (8) (2006) 735–740.
doi: https://doi.org/10.1016/j.aml.2005.08.022.
URL: https://www.sciencedirect.com/science/article/pii/S0893965905003319
"""
function solve_adbc(
    equation::DiffusionEquation{ClosedIntervalDomain{T},TopHatFilter},
    u,
    tlist,
    n,
    Δt = (tlist[2] - tlist[1]) / 1000,
) where {T}
    @unpack domain, filter, f, g_a, g_b = equation
    x = discretize(domain, n)
    Δx = x[2] - x[1]
    h = filter.width
    all(h.(x) .≈ h(x[1])) || error("ADBC requires constant filter width")
    h₀ = h(x[1])
    a, b = domain.left, domain.right

    # Get matrices
    D = diffusion_matrix(domain, n)

    # Filter matrix
    inds = [-1, 0, 1]
    stencil = [1 / 24, 11 / 12, 1 / 24]
    diags = [i => fill(s, n + 1 - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    W = spdiagm(diags...)
    W[1, [1, 2]] = [23 / 24, 1 / 24]
    W[end, [end, end - 1]] = [23 / 24, 1 / 24]
    w₀ = W[1, :]
    wₙ = W[end, :]

    # Reconstruction matrix
    inds = [-1, 0, 1]
    stencil = [-1 / 24, 13 / 12, -1 / 24]
    diags = [i => fill(s, n + 1 - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    R = spdiagm(diags...)
    R[1, [1, 2]] = [25 / 24, -1 / 24]
    R[end, [end, end - 1]] = [25 / 24, -1 / 24]

    function perform_step!(ūᵏ, tᵏ, Δt, p)
        uᵏ, uᵏ⁺¹, ūᵏ⁺¹ = p

        # Current approximate unfiltered solution
        mul!(uᵏ, R, ūᵏ)

        # Next approximate unfiltered solution
        uᵏ⁺¹ .= uᵏ
        uᵏ⁺¹[1] = g_a(tᵏ + Δt)
        uᵏ⁺¹[end] = g_b(tᵏ + Δt)

        # Filtered boundary conditions
        ūᵏ⁺¹[1] = w₀'uᵏ⁺¹
        ūᵏ⁺¹[end] = wₙ'uᵏ⁺¹

        # Next inner points for filtered solution
        ūᵏ⁺¹[2:end-1] .+=
            Δt .* D[2:end-1, :]ūᵏ .+ Δt .* W[2:end-1, :]f.(x, tᵏ) .+
            Δt .* (abs.(x[2:end-1] .- b) .≤ h₀) ./ 2h₀ .* (g_b(tᵏ) - uᵏ[end-1]) / Δx .-
            Δt .* (abs.(x[2:end-1] .- a) .≤ h₀) ./ 2h₀ .* (uᵏ[2] - g_a(tᵏ)) / Δx

        # Advance by Δt
        ūᵏ .= ūᵏ⁺¹
    end

    ūᵏ = W * u.(x)
    p = (copy(ūᵏ), copy(ūᵏ), copy(ūᵏ))
    tᵏ = tlist[1]
    while tᵏ + Δt < tlist[2]
        perform_step!(ūᵏ, tᵏ, Δt, p)
        tᵏ += Δt
    end

    # Perform last time step (with adapted step size)
    Δt_last = tlist[2] - tᵏ
    perform_step!(ūᵏ, tᵏ, Δt_last, p)
    tᵏ += Δt_last

    ūᵏ
end

function solve_adbc(
    equation::DiffusionEquation{ClosedIntervalDomain{T},ConvolutionalFilter},
    u,
    tlist,
    n,
    Δt = (tlist[2] - tlist[1]) / 1000;
    solver = :euler,
    δ = √(6 / π) / equation.filter.kernel(0),
    ū_ext = nothing,
    uₓ = nothing,
) where {T}
    @unpack domain, filter, f, g_a, g_b = equation
    x = discretize(domain, n)
    Δx = x[2] - x[1]
    h₀ = filter.width(0)
    G = filter.kernel
    a, b = domain.left, domain.right

    # Get matrices
    D = diffusion_matrix(domain, n)

    # Filter stencil
    normalize(w) = w / sum(w)
    nw = floor(Int, h₀ / 2Δx)
    inds = -nw:nw
    w = normalize(G.(Δx .* inds))

    # Filter matrix
    diags = [i => fill(w, n + 1 - abs(i)) for (i, w) ∈ zip(inds, w)]
    W = spdiagm(diags...)

    # Move weights of outside points to endpoint
    for i = 1:nw
        W[i, 1] += sum(w[1:(nw+1-i)])
        W[end+1-i, end] += sum(w[(end-nw+i):end])
    end

    # Boundary stencils
    w₀ = W[1, :]
    wₙ = W[end, :]

    G_a = G.(x[2:end-1] .- a)
    G_b = G.(b .- x[2:end-1])

    f̄ᵏ_constant = apply_filter_extend(x -> f(x, tlist[1]), filter, domain).(x[2:end-1])

    function du!(du, ūᵏ, tᵏ, Δt, p)
        @unpack uᵏ, uᵏ⁺¹, Dūᵏ, fᵏ, f̄ᵏ, δ = p

        # Diffusion term
        mul!(Dūᵏ, D, ūᵏ)

        # Source term
        # fᵏ .= f.(x, tᵏ)
        # mul!(f̄ᵏ, W[2:end-1, :], fᵏ)
        # f̄ᵏ .= apply_filter_extend(x -> f(x, tᵏ), filter, domain).(x[2:end-1])
        f̄ᵏ = f̄ᵏ_constant

        # Current approximate unfiltered solution
        @. uᵏ = ūᵏ - δ^2 / 24 * Dūᵏ

        # Next approximate unfiltered solution
        uᵏ⁺¹ .= uᵏ
        uᵏ⁺¹[1] = g_a(tᵏ + Δt)
        uᵏ⁺¹[end] = g_b(tᵏ + Δt)

        # Filtered boundary conditions
        if isnothing(ū_ext)
            # Approximate deconvolution boundary conditions
            du[1] = (w₀'uᵏ⁺¹ - ūᵏ[1]) / Δt
            du[end] = (wₙ'uᵏ⁺¹ - ūᵏ[end]) / Δt
        else
            # Exact filtered boundary conditions
            du[1] = (ū_ext(tᵏ + Δt)(a) - ūᵏ[1]) / Δt
            du[end] = (ū_ext(tᵏ + Δt)(b) - ūᵏ[end]) / Δt
        end

        if isnothing(uₓ)
            # Approximate boundary terms
            uₓ_a = (uᵏ[2] - g_a(tᵏ)) / Δx
            uₓ_b = (g_b(tᵏ) - uᵏ[end-1]) / Δx
        else
            # Exact boundary terms
            uₓ_a = uₓ(a, tᵏ)
            uₓ_b = uₓ(b, tᵏ)
        end

        # Next inner points for filtered solution
        du_in = @view du[2:end-1]
        Dūᵏ_in = @view Dūᵏ[2:end-1]
        @. du_in = Dūᵏ_in + f̄ᵏ + G_b * uₓ_b - G_a * uₓ_a
    end

    function perform_step_euler!(ūᵏ, tᵏ, Δt, p)
        @unpack du = p
        du!(du, ūᵏ, tᵏ, Δt, p)
        ūᵏ .+= Δt .* du
    end

    function perform_step_RK4!(ūᵏ, tᵏ, Δt, p)
        @unpack k₁, k₂, k₃, k₄ = p

        # Compute change
        du!(k₁, ūᵏ, tᵏ, Δt / 2, p)
        du!(k₂, ūᵏ + Δt / 2 * k₁, tᵏ + Δt / 2, Δt / 2, p)
        du!(k₃, ūᵏ + Δt / 2 * k₂, tᵏ + Δt / 2, Δt / 2, p)
        du!(k₄, ūᵏ + Δt * k₃, tᵏ + Δt, Δt / 2, p)

        # Advance by Δt
        @. ūᵏ += Δt * (k₁ / 6 + k₂ / 3 + k₃ / 3 + k₄ / 6)
    end

    if solver == :euler
        perform_step! = perform_step_euler!
    elseif solver == :RK4
        perform_step! = perform_step_RK4!
    else
        error("Solver is not supported")
    end

    # ūᵏ = W * u.(x)
    ūᵏ = apply_filter_extend(u, filter, domain).(x)
    p = (;
        du = copy(ūᵏ),
        k₁ = copy(ūᵏ),
        k₂ = copy(ūᵏ),
        k₃ = copy(ūᵏ),
        k₄ = copy(ūᵏ),
        uᵏ = copy(ūᵏ),
        uᵏ⁺¹ = copy(ūᵏ),
        Dūᵏ = copy(ūᵏ),
        fᵏ = copy(ūᵏ),
        f̄ᵏ = copy(ūᵏ[2:end-1]),
        δ,
    )
    tᵏ = tlist[1]
    while tᵏ + Δt < tlist[2]
        # @show tᵏ
        perform_step!(ūᵏ, tᵏ, Δt, p)
        tᵏ += Δt
    end

    # Perform last time step (with adapted step size)
    Δt_last = tlist[2] - tᵏ
    perform_step!(ūᵏ, tᵏ, Δt_last, p)
    tᵏ += Δt_last

    # nT = round(Int, (tlist[2] - tlist[1]) / Δt)
    # uu = zeros(n + 1, nT + 1)
    # uu[:, 1] = ūᵏ
    # for i = 1:nT
    #     perform_step!(ūᵏ, tᵏ, Δt, p)
    #     uu[:, i+1] = ūᵏ
    #     tᵏ += Δt
    # end

    ūᵏ
    # uu
end


text:
using Dynstack
using Test

@testset "Dynstack.jl" begin
    # Write your tests here.
end


text:
function MakieLayout.Legend(fg::FigureGrid)
    colorbar = _Colorbar_(fg)
    legend = _Legend_(fg)
    if !isnothing(colorbar)
        Colorbar(fg.figure[:, end + 1]; colorbar...)
    end
    if !isnothing(legend)
        Legend(fg.figure[:, end + 1], legend...)
    end
end

function has_zcolor(entry::Entry)
    return entry.plottype <: Union{Heatmap, Contour, Contourf, Surface} &&
        !haskey(entry.primary, :color) &&
        !haskey(entry.named, :color) &&
        !haskey(entry.attributes, :color)
end

function getlabeledcolorbar(grid)
    scales, labels = first(grid).scales, first(grid).labels
    key = any(has_zcolor, entries(grid)) ? 3 : :color
    label, scale = get(labels, key, nothing), get(scales, key, nothing)
    return scale isa ContinuousScale ? (label, scale) : nothing
end

function _Colorbar_(fg::FigureGrid)
    grid = fg.grid
    labeledcolorbar = getlabeledcolorbar(grid)
    isnothing(labeledcolorbar) && return
    label, colorscale = labeledcolorbar
    colormap = current_default_theme().colormap[]
    for entry in entries(grid)
        colormap = to_value(get(entry.attributes, :colormap, colormap))
    end
    limits = colorscale.extrema
    return (; label, limits, colormap)
end

"""
    plottypes_attributes(entries)

Return plottypes and relative attributes, as two vectors of the same length,
for the given `entries`.
"""
function plottypes_attributes(entries)
    plottypes = PlotFunc[]
    attributes = Vector{Symbol}[]
    for entry in entries
        # FIXME: this should probably use the rescaled values
        plottype = Makie.plottype(entry.plottype, entry.positional...)
        n = findfirst(==(plottype), plottypes)
        attrs = (keys(entry.primary)..., keys(entry.named)...)
        if isnothing(n)
            push!(plottypes, plottype)
            push!(attributes, collect(Symbol, attrs))
        else
            union!(attributes[n], attrs)
        end
    end
    return plottypes, attributes
end

hassymbolkey((k, v)::Pair) = k isa Symbol

function _Legend_(fg::FigureGrid)
    grid = fg.grid

    # assume all subplots have same scales, to be changed to support free scales
    scales = filter(hassymbolkey, first(grid).scales)
    labels = filter(hassymbolkey, first(grid).labels)

    # remove keywords that don't support legends
    for key in [:row, :col, :layout, :stack, :dodge, :group]
        pop!(labels, key, nothing)
    end
    for (key, val) in scales
        val isa ContinuousScale && pop!(labels, key, nothing)
    end

    # if no legend-worthy keyword remains return nothing
    isempty(labels) && return nothing

    titles = unique!(collect(String, values(labels)))
    # empty strings create difficulties with the layout
    nonemptytitles = map(t -> isempty(t) ? " " : t, titles)

    plottypes, attributes = plottypes_attributes(entries(grid))

    labels_list = Vector{String}[]
    elements_list = Vector{Vector{LegendElement}}[]

    for title in titles
        label_attrs = [key for (key, val) in labels if val == title]
        uniquevalues = mapreduce(k -> scales[k].data, assert_equal, label_attrs)
        elements = map(eachindex(uniquevalues)) do idx
            local elements = LegendElement[]
            for (P, attrs) in zip(plottypes, attributes)
                shared_attrs = attrs ∩ label_attrs
                isempty(shared_attrs) && continue
                options = [attr => scales[attr].plot[idx] for attr in shared_attrs]
                append!(elements, legend_elements(P; options...))
            end
            return elements
        end
        push!(labels_list, map(string, uniquevalues))
        push!(elements_list, elements)
    end
    return elements_list, labels_list, nonemptytitles
end

#Notes

# TODO: correctly handle composite plot types (now fall back to poly)
# TODO: check that all scales for the same label agree on the data
# TODO: make legend updateable?
# TODO: allow custom attributes in legend elements?


text:
# Packages used in Script
using Luxor
using Colors

# Auxiliary drawing functions
function make_drawing(width, height, img_path, bkg_color, origin_p)
    d = Drawing(width, height, img_path)
    background(bkg_color)
    origin(origin_p)
    return d
end

map2luxor(p) = Point(p.x, -p.y)

# Constants
## Canvas constants
width = 500
height = 500
path = "voronoi.png"
color = RGBA(.5, 0, .5)
op = Point(width / 2, height / 2)

## Point constants
p1 = (x = -125, y = 125)
p2 = (x = 125, y = 125)

# Setting up drawing
my_draw = make_drawing(width, height, path, color, op)

# Mathematical formulae for calculations
midpoint(x₁, y₁, x₂, y₂) = (x = (x₁ + x₂) / 2, y = (y₁ + y₂) / 2)
slope(x₁, y₁, x₂, y₂) = (y₂ - y₁) / (x₂ - x₁)
point_slope(x, x₁, y₁, m) = (x = x, y = m * (x - x₁) + y₁)
d(x1, y1, x2, y2, x, y) = (x - x1) * (y2 - y1) - (y - y1) * (x2 - x1)

#=

Steps to create Voronoi diagram

1. Place points in an xy axis
2. Starting at one point, determine the midpoint between these points
3. Calculate the slope between these two points
4. Calculate the slope's reverse reciprocal
5. Find two points that lie on the line formed by the reverse reciprocal and midpoint
6. Calculate all points that lie in either half plane formed by the perpendicular bisector

=#

sethue(RGBA(1, 1, 0))
circle(map2luxor(p1), 15; action = :fill)
circle(map2luxor(p2), 15; action = :fill)

sethue(RGBA(1, 1, 1))
line(map2luxor(p1), map2luxor(p2); action = :stroke)

mp = midpoint(p1.x, p1.y, p2.x, p2.y)

sethue(RGBA(0, 0, 0))
circle(map2luxor(mp), 7.5; action = :fill)

sp = slope(p1.x, p1.y, p2.x, p2.y)

rr = sp == 0 ? 0 : -1 / sp
sp1 = point_slope(width / 2, mp.x, mp.y, rr)
sp2 = point_slope(-width / 2, mp.x, mp.y, rr)

line(map2luxor(sp1), map2luxor(sp2); action = :stroke)

h1_plane = []
h2_plane = []
h3_plane = []
for x in (-width/2):(width/2), y in (-height/2):(height/2)
	if d(sp1.x, sp1.y, sp2.x, sp2.y, x, y) < 0
		p = (x = x, y = y)
		push!(h1_plane, map2luxor(p))
	elseif d(sp1.x, sp1.y, sp2.x, sp2.y, x, y) > 0
		p = (x = x, y = y)
		push!(h2_plane, map2luxor(p))
	else
		p = (x = x, y = y)
		push!(h3_plane, map2luxor(p))
	end
end

setcolor(1, 0, 0, 1)
for c in h1_plane
	circle(c, 1; action = :stroke)
end

setcolor(0, 0, 1, 1)
for c in h2_plane
	circle(c, 1; action = :stroke)
end

setcolor(1, 1, 1, 1)
for c in h3_plane
	circle(c, 1; action = :stroke)
end

setcolor(0, 0, 1, 1)
for plane in [h1_plane, h2_plane, h3_plane]
	if map2luxor(p1) in plane
		circle.(plane, 1; action = :stroke)
	end
end

finish()

max_stars_repo_name: TheCedarPrince/DataExplorers

max_issues_repo_name: TheCedarPrince/DataExplorers

max_forks_repo_name: TheCedarPrince/DataExplorers


text:
using TOML

using Test
using Dates

const jsnval = Dict{String,Function}(
    "string" =>identity,
    "float"    => (s -> Base.parse(Float64, s)),
    "integer"  => (s -> Base.parse(Int64, s)),
    "datetime" => (s -> Base.parse(DateTime, endswith(s, 'Z') ? chop(s) : s)),
    "datetime-local" => (s -> Base.parse(DateTime, endswith(s, 'Z') ? chop(s) : s)),
    "date-local" => (s -> Base.parse(DateTime, endswith(s, 'Z') ? chop(s) : s)),
    "time-local" => (s -> Base.parse(Time, s)),
    "array"    => (a -> map(jsn2data, a)),
    "bool"     => (b -> b == "true")
)

function jsn2data(jsn)
    if "type" in keys(jsn)
        jsnval[jsn["type"]](jsn["value"])
    elseif jsn isa Vector
        [jsn2data(v) for v in jsn]
    else
        Dict{String,Any}([k => jsn2data(v) for (k, v) in jsn])
    end
end

#########
# Valid #
#########

function check_valid(f)
    jsn = try jsn2data(@eval include($f * ".jl"))
    # Some files cannot be reprsented with julias DateTime (timezones)
    catch
        return false
    end
    tml = TOML.tryparsefile(f * ".toml")
    tml isa TOML.Internals.ParserError && return false
    return isequal(tml, jsn)
end

@testset "valid" begin

failures = [
    "testfiles/valid/spec-example-1-compact.toml",
    "testfiles/valid/spec-example-1.toml",
    "testfiles/valid/comment/tricky.toml",
    "testfiles/valid/datetime/datetime.toml",
    "testfiles/valid/datetime/milliseconds.toml",
    "testfiles/valid/datetime/timezone.toml",
    "testfiles/valid/float/zero.toml",
    "testfiles/valid/string/multiline-quotes.toml",
]

valid_test_folder = joinpath(@__DIR__, "testfiles", "valid")
for (root, dirs, files) in walkdir(valid_test_folder)
    for f in files
        if endswith(f, ".toml")
            file = joinpath(root, f)
            rel = relpath(file, @__DIR__)
            if Sys.iswindows()
                rel = replace(rel, '\\' => '/')
            end
            v = check_valid(splitext(file)[1])
            if rel in failures
                @test_broken v
            else
                @test v
            end
        end
    end
end

end # testset

###########
# Invalid #
###########

# TODO: Check error type
function check_invalid(f)
    tml = try TOML.tryparsefile(f)
    catch
        return false
    end
    return tml isa TOML.Internals.ParserError
end

@testset "invalid" begin

failures = [
    "testfiles/invalid/control/comment-del.toml",
    "testfiles/invalid/control/comment-lf.toml",
    "testfiles/invalid/control/comment-null.toml",
    "testfiles/invalid/control/comment-us.toml",
    "testfiles/invalid/control/multi-del.toml",
    "testfiles/invalid/control/multi-lf.toml",
    "testfiles/invalid/control/multi-null.toml",
    "testfiles/invalid/control/multi-us.toml",
    "testfiles/invalid/control/rawmulti-del.toml",
    "testfiles/invalid/control/rawmulti-lf.toml",
    "testfiles/invalid/control/rawmulti-null.toml",
    "testfiles/invalid/control/rawmulti-us.toml",
    "testfiles/invalid/control/rawstring-del.toml",
    "testfiles/invalid/control/rawstring-lf.toml",
    "testfiles/invalid/control/rawstring-null.toml",
    "testfiles/invalid/control/rawstring-us.toml",
    "testfiles/invalid/control/string-bs.toml",
    "testfiles/invalid/control/string-del.toml",
    "testfiles/invalid/control/string-lf.toml",
    "testfiles/invalid/control/string-null.toml",
    "testfiles/invalid/control/string-us.toml",
    "testfiles/invalid/encoding/bad-utf8-in-comment.toml",
    "testfiles/invalid/encoding/bad-utf8-in-string.toml",
    "testfiles/invalid/integer/negative-bin.toml",
    "testfiles/invalid/integer/negative-hex.toml",
    "testfiles/invalid/integer/negative-hex.toml",
    "testfiles/invalid/integer/negative-oct.toml",
    "testfiles/invalid/integer/positive-bin.toml",
    "testfiles/invalid/integer/positive-hex.toml",
    "testfiles/invalid/integer/positive-oct.toml",
    "testfiles/invalid/key/multiline.toml",
]

invalid_test_folder = joinpath(@__DIR__, "testfiles", "invalid")
for (root, dirs, files) in walkdir(invalid_test_folder)
    for f in files
        if endswith(f, ".toml")
            file = joinpath(root, f)
            rel = relpath(file, @__DIR__)
            if Sys.iswindows()
                rel = replace(rel, '\\' => '/')
            end
            v = check_invalid(file)
            if rel in failures
                @test_broken v
            else
                @test v
            end
        end
    end
end

end # testset


text:
IntegerVector{I<:Integer} = Vector{I}
RealVector{N<:Real} = Vector{N}
RealMatrix{N<:Real} = Matrix{N}

RealLowerTriangular{T, S} = LowerTriangular{T,S} where {T<:Real, S<:AbstractMatrix{T}}
RealLowerTriangular(m) = LowerTriangular(m)

RealNormal{N<:Real} = Normal{N}
MultivariateGMM{D<:MvNormal} = MultivariateMixture{Continuous, D}

multivecs{T}(::Type{T}, n::Int) = [T[] for _ =1:n]


text:
#=
Form factor + edge velocity calculation
=#

################################################
"""
$(SIGNATURES)

Calculate node volume  and voronoi surface contributions for cell.
"""

function cellfactors!(T::Type{Edge1D},::Type{Cartesian1D},coord,cellnodes,icell,nodefac,edgefac)
    en=local_celledgenodes(T)
    K=cellnodes[en[1,1],icell]
    L=cellnodes[en[2,1],icell]
    xK=coord[1,K]
    xL=coord[1,L]
    d=abs(xL[1]-xK[1])
    nodefac[1]=d/2
    nodefac[2]=d/2
    edgefac[1]=1/d
    nothing
end

function cellfactors!(::Type{Edge1D},::Type{<:Polar1D}, coord,cellnodes,icell,nodefac,edgefac)
    en=local_celledgenodes(T)
    K=cellnodes[en[1,1],icell]
    L=cellnodes[en[2,1],icell]
    xK=coord[1,K]
    xL=coord[1,L]
    r0=xK[1]
    r1=xL[1]
    if r1<r0
        r0=xL[1]
        r1=xK[1]
    end
    rhalf=0.5*(r1+r0);
    # cpar[1]= π*(r1*r1-r0*r0);         # circular volume
    nodefac[1]= π*(rhalf*rhalf-r0*r0);   # circular volume between midline and boundary
    nodefac[2]= π*(r1*r1-rhalf*rhalf);   # circular volume between midline and boundary
    edgefac[1]= 2.0*π*rhalf/(r1-r0);     # circular surface / width
    nothing
end

function cellfactors!(T::Type{Triangle2D},::Type{Cartesian2D},coord,cellnodes,icell,npar,epar)
    en=local_celledgenodes(T)
    @views n=SVector{3}(cellnodes[:,icell])
    V=@MMatrix zeros(2,3)
    dd=@MVector zeros(3)

    # Fill matrix of edge vectors
    V[1,3]= coord[1,n[en[1,3]]]- coord[1,n[en[2,3]]]
    V[2,3]= coord[2,n[en[1,3]]]- coord[2,n[en[2,3]]]

    V[1,2]= coord[1,n[en[1,2]]]- coord[1,n[en[2,2]]]
    V[2,2]= coord[2,n[en[1,2]]]- coord[2,n[en[2,2]]]

    V[1,1]= coord[1,n[en[1,1]]]- coord[1,n[en[2,1]]]
    V[2,1]= coord[2,n[en[1,1]]]- coord[2,n[en[2,1]]]

    # Compute determinant
    det=V[1,3]*V[2,2] - V[1,2]*V[2,3]
    vol=abs(0.5*det)
    ivol = 1.0/vol

    # squares of edge lengths
    dd[3]=V[1,3]*V[1,3]+V[2,3]*V[2,3]
    dd[2]=V[1,2]*V[1,2]+V[2,2]*V[2,2]
    dd[1]=V[1,1]*V[1,1]+V[2,1]*V[2,1]

    # contributions to \sigma_kl/h_kl
    epar[1]= (dd[2]+dd[3]-dd[1])*0.125*ivol
    epar[2]= (dd[3]+dd[1]-dd[2])*0.125*ivol
    epar[3]= (dd[1]+dd[2]-dd[3])*0.125*ivol

    # contributions to \omega_k

    npar.=0.0
    for i=1:3
        npar[en[1,i]]+=epar[i]*dd[i]*0.25
        npar[en[2,i]]+=epar[i]*dd[i]*0.25
    end

    nothing
end

function cellfactors!(T::Type{Triangle2D},::Type{<:Cylindrical2D},coord,cellnodes,icell,npar,epar)

    function area2d(coord1, coord2, coord3)
        V11= coord2[1]- coord1[1]
        V21= coord2[2]- coord1[2]

        V12= coord3[1]- coord1[1]
        V22= coord3[2]- coord1[2]

        V13= coord3[1]- coord2[1]
        V23= coord3[2]- coord2[2]

        # Compute determinant
        det=V11*V22 - V12*V21
        area=abs(0.5*det)
    end

    en=local_celledgenodes(T)

    @views n=SVector{3}(cellnodes[:,icell])
    V=@MMatrix zeros(2,3)
    dd=@MVector zeros(3)
    cc=@MVector zeros(2)
    emid=@MMatrix zeros(2,3)

    # Fill matrix of edge vectors
    V[1,3]= coord[1,n[en[1,3]]]- coord[1,n[en[2,3]]]
    V[2,3]= coord[2,n[en[1,3]]]- coord[2,n[en[2,3]]]

    V[1,2]= coord[1,n[en[1,2]]]- coord[1,n[en[2,2]]]
    V[2,2]= coord[2,n[en[1,2]]]- coord[2,n[en[2,2]]]

    V[1,1]= coord[1,n[en[1,1]]]- coord[1,n[en[2,1]]]
    V[2,1]= coord[2,n[en[1,1]]]- coord[2,n[en[2,1]]]

    # Compute determinant
    det=V[1,3]*V[2,2] - V[1,2]*V[2,3]
    area=abs(0.5*det)

    # Integrate R over triangle (via quadrature rule)
    vol=2.0*π*area*(coord[1,n[1]]+coord[1,n[2]]+coord[1,n[3]])/3.0

    # squares of edge lengths
    dd[3]=V[1,3]*V[1,3]+V[2,3]*V[2,3]
    dd[2]=V[1,2]*V[1,2]+V[2,2]*V[2,2]
    dd[1]=V[1,1]*V[1,1]+V[2,1]*V[2,1]

    for i=1:3
        emid[1,i]=0.5*(coord[1,n[en[1,i]]]+coord[1,n[en[2,i]]])
        emid[2,i]=0.5*(coord[2,n[en[1,i]]]+coord[2,n[en[2,i]]])
    end
    @views tricircumcenter!(cc,coord[:,n[1]],coord[:,n[2]],coord[:,n[3]])

    r(p)=p[1]
    z(p)=p[2]

    for i=1:3
        @views epar[i]= π*(r(cc)+r(emid[:,i]))*sqrt((r(cc)-r(emid[:,i]))^2+(z(cc)-z(emid[:,i]))^2)/sqrt(dd[i]);
    end

    rintegrate(coord1, coord2, coord3)=2.0*π*area2d(coord1,coord2,coord3)*(coord1[1]+coord2[1]+coord3[1])/3.0

    npar.=0.0
    for i=1:3
        @views npar[en[1,i]]+=rintegrate(coord[:,n[en[1,i]]],cc,emid[:,i])
        @views npar[en[2,i]]+=rintegrate(coord[:,n[en[2,i]]],cc,emid[:,i])
    end
    nothing
end

function cellfactors!(T::Type{Tetrahedron3D},::Type{Cartesian3D},coord,cellnodes,icell,npar,epar)
    # Transferred from WIAS/pdelib, (c) J. Fuhrmann, H. Langmach, I. Schmelzer
    es=@MVector zeros(6)
    dd=@MVector zeros(6)
    @views n=SVector{4}(cellnodes[:,icell])

    en=local_celledgenodes(T)

    @views pp1=SVector{6}(en[1,:])
    @views pp2=SVector{6}(en[2,:])

#    pp1=(1,1,1,2,2,3)
#    pp2=(2,3,4,3,4,4)

    pi1=(5,6,5,1)
    pi2=(6,3,1,4)
    pi3=(4,2,3,2)

    po1=(2,1,2,6)
    po2=(1,4,6,3)
    po3=(3,5,4,5)
#    pedge=(6,5,4,3,2,1)

    for i=1:6
        p1=cellnodes[pp1[i],icell]
        p2=cellnodes[pp2[i],icell]
        dx=coord[1,p1]-coord[1,p2]
        dy=coord[2,p1]-coord[2,p2]
        dz=coord[3,p1]-coord[3,p2]
        dd[i]=dx*dx+dy*dy+dz*dz
        es[i]=0.0
        epar[i]=0.0
    end

    x1=coord[1,n[2]]-coord[1,n[1]]
    y1=coord[2,n[2]]-coord[2,n[1]]
    z1=coord[3,n[2]]-coord[3,n[1]]

    x2=coord[1,n[3]]-coord[1,n[1]]
    y2=coord[2,n[3]]-coord[2,n[1]]
    z2=coord[3,n[3]]-coord[3,n[1]]

    x3=coord[1,n[4]]-coord[1,n[1]]
    y3=coord[2,n[4]]-coord[2,n[1]]
    z3=coord[3,n[4]]-coord[3,n[1]]

    det= (x1*(y2*z3-y3*z2) + x2*(y3*z1-y1*z3) + x3*(y1*z2-y2*z1))

    if det <0
        det=-det
    end
    vol=det/6
    vv= 96*6*vol

    for  i=1:4
        npar[i]=0.0
        i1=pi1[i]
        i2=pi2[i]
        i3=pi3[i]

        h1=dd[i1]*(dd[i2]+dd[i3]-dd[i1])
        h2=dd[i2]*(dd[i3]+dd[i1]-dd[i2])
        h3=dd[i3]*(dd[i1]+dd[i2]-dd[i3])

        df=h1+h2+h3
        vf=(h1*dd[po1[i]] + h2*dd[po2[i]] + h3*dd[po3[i]]-2*dd[i1]*dd[i2]*dd[i3]) / (vv*df)
        es[i1] += h1*vf
        es[i2] += h2*vf
        es[i3] += h3*vf
    end

    for i=1:6
        epar[i] = 6*es[i]/dd[i]
        npar[pp1[i]] += es[i]
        npar[pp2[i]] += es[i]
    end
    nothing
end

################################################
"""
$(SIGNATURES)

Calculate node volume  contributions for boundary face.
"""

function bfacefactors!(::Type{Vertex0D},::Type{Cartesian1D},coord,bfacenodes,ibface,nodefac,edgefac)
    nodefac[1]=1.0
    nothing
end

function bfacefactors!(::Type{Vertex0D},::Type{<:Polar1D},coord,bfacenodes,ibface,nodefac, edgefac)
    inode=bfacenodes[1,ibface]
    r=coord[1,inode]
    nodefac[1]=2*π*r
    nothing
end

function bfacefactors!(T::Type{Edge1D},::Type{<:Cartesian2D},coord,bfacenodes,ibface,nodefac, edgefac)
    en         = local_celledgenodes(T)
    i1         = bfacenodes[en[1,1],ibface]
    i2         = bfacenodes[en[2,1],ibface]
    dx         = coord[1,i1]-coord[1,i2]
    dy         = coord[2,i1]-coord[2,i2]
    d          = sqrt(dx*dx+dy*dy)
    nodefac[1] = d/2
    nodefac[2] = d/2
    edgefac[1] = 1/d
    nothing
end

function bfacefactors!(T::Type{Edge1D},::Type{<:Cylindrical2D},coord,bfacenodes,ibface,nodefac,edgefac)
    en=local_celledgenodes(T)
    i1=bfacenodes[en[1,1],ibface]
    i2=bfacenodes[en[2,1],ibface]
    r1=coord[1,i1]
    r2=coord[1,i2]
    z1=coord[2,i1]
    z2=coord[2,i2]
    dr=r1-r2
    rmid=(r1+r2)/2
    dz=z1-z2
    l=sqrt(dr*dr+dz*dz)
    nodefac[1]=π*(r1+rmid)*l/2
    nodefac[2]=π*(r2+rmid)*l/2
    nothing
end

function bfacefactors!(T::Type{Triangle2D}, ::Type{<:Cartesian3D}, coord, bfacenodes, ibface, npar, epar)

    en=local_celledgenodes(T)
    @views n=SVector{3}(bfacenodes[:,ibface])
    dd = @MVector zeros(3)

    # Transferred from WIAS/pdelib, (c) J. Fuhrmann, H. Langmach, I. Schmelzer
    for j=1:3
        d = coord[j,n[en[1,1]]] - coord[j,n[en[2,1]]];  dd[1]+= d*d;
        d = coord[j,n[en[1,2]]] - coord[j,n[en[2,2]]];  dd[2]+= d*d;
        d = coord[j,n[en[1,3]]] - coord[j,n[en[2,3]]];  dd[3]+= d*d;
    end

    # Kanten-Flaechenanteile (ohne Abschneiden); epar als Hilfsfeld benutzt
    epar[1] = (dd[2]+dd[3]-dd[1])*dd[1];
    epar[2] = (dd[3]+dd[1]-dd[2])*dd[2];
    epar[3] = (dd[1]+dd[2]-dd[3])*dd[3];
    vol     = sqrt(epar[1]+epar[2]+epar[3])*0.25;

    d = 1.0/(8*vol);

    # Knoten-Flaechenanteile (ohne Abschneiden)
    npar[1] = (epar[3]+epar[2])*d*0.25;
    npar[2] = (epar[1]+epar[3])*d*0.25;
    npar[3] = (epar[2]+epar[1])*d*0.25;

    # Kantengewichte
    epar[1] = epar[1]*d/dd[1];
    epar[2] = epar[2]*d/dd[2];
    epar[3] = epar[3]*d/dd[3];
    nothing
end

##################################################################

#
# TODO: this should be generalized for more quadrules
#
function integrate(coordl,coordr,hnormal,velofunc)
    wl=1.0/6.0
    wm=2.0/3.0
    wr=1.0/6.0
    coordm=0.5*(coordl+coordr)
    (vxl,vyl)=velofunc(coordl[1],coordl[2])
    (vxm,vym)=velofunc(coordm[1],coordm[2])
    (vxr,vyr)=velofunc(coordr[1],coordr[2])
    return (wl*vxl +wm*vxm+wr*vxr)*hnormal[1] + (wl*vyl +wm*vym+wr*vyr)*hnormal[2]
end

"""
$(SIGNATURES)

Project velocity onto grid edges,
"""
function edgevelocities(grid,velofunc)

    @assert dim_space(grid)<3

    cn=grid[CellNodes]
    ec=grid[EdgeCells]
    en=grid[EdgeNodes]
    coord=grid[Coordinates]

    velovec=zeros(Float64,num_edges(grid))
    if dim_space(grid)==1
        for iedge=1:num_edges(grid)
            K=en[1,iedge]
            L=en[2,iedge]
            elen=coord[1,L]-coord[1,K]
            vx,vy=velofunc((coord[1,K]+coord[1,L])/2,0)
            velovec[iedge]=-elen*vx
        end
    else
        for iedge=1:num_edges(grid)
            K=en[1,iedge]
            L=en[2,iedge]
            p1=@MVector zeros(2)
            p2=@MVector zeros(2)
            tricircumcenter!(p1,
                             coord[:,cn[1,ec[1,iedge]]],
                             coord[:,cn[2,ec[1,iedge]]],
                             coord[:,cn[3,ec[1,iedge]]])
            if ec[2,iedge]>0
                tricircumcenter!(p2,
                                 coord[:,cn[1,ec[2,iedge]]],
                                 coord[:,cn[2,ec[2,iedge]]],
                                 coord[:,cn[3,ec[2,iedge]]])
            else
                p2.=0.5*(coord[:,K]+coord[:,L])
            end
            hnormal=coord[:,K]-coord[:,L]
            velovec[iedge]=integrate(p1,p2,hnormal,velofunc)
        end
    end
    return velovec
end

function bfacevelocities(grid,velofunc)
    @assert dim_space(grid)<3
    bfacenodes=grid[BFaceNodes]
    coord=grid[Coordinates]
    bfacecells=grid[BFaceCells]
    bfacenormals=grid[BFaceNormals]
    bfr=grid[BFaceRegions]
    velovec=zeros(Float64,2,num_bfaces(grid))
    if dim_space(grid)==1
        for ibface=1:num_bfaces(grid)
            vx,vy=velofunc(coord[1,bfacenodes[1,ibface]])
            velovec[ibface]=vx*bfacenormals[1,ibface]
        end
    else
        for ibface=1:num_bfaces(grid)
            p1=coord[:,bfacenodes[1,ibface]]
            p2=coord[:,bfacenodes[2,ibface]]
            pm=0.5*(p1+p2)
            velovec[1,ibface]=integrate(p1,pm,bfacenormals[:,ibface],velofunc)
            velovec[2,ibface]=integrate(pm,p2,bfacenormals[:,ibface],velofunc)
        end
    end
    return velovec
end


text:
using ThreeBodyDecay
using Test
# using Plots

@testset "Scattering Length Approximation" begin
    mJψ = 3.0969;  mππ = 0.3;
    mD0 = 1.86483; mDstar0 = 2.00685

    cs = (c11 = -0.03, c12 = 0.07, c22=0.0)
    #
    ScA = ScattLenApproximation(
        Ms=(cs.c11,cs.c12,cs.c22),
        ms1=(mD0, mDstar0), ms2=(mJψ,mππ))
    #
    ev = range(-15,10, length=300) # MeV
    sv = ((mD0+mDstar0) .+ ev .* 1e-3).^2 # GeV
    #
    calv = map(x->amp(x, ScA), sv.+1e-7im) # amplitude
    v,i = findmax(imag.(calv))
    #
    # plot(sqrt.(sv), [real.(calv) imag.(calv)])
    # vline!([1.86483+2.00685])
    #
    # peaking structure
    @test abs(-1.705-ev[i]) < 0.001
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay