export HH3DHyperSingularFDBIO

abstract type Helmholtz3DOp <: MaxwellOperator3D end
abstract type Helmholtz3DOpReg <: MaxwellOperator3DReg end

"""
```
∫_Γ dx ∫_Γ dy \left(α G g(x) n_x ⋅ n_y f(y) + β G \mbox{curl} g(x) ⋅ \mbox{curl} f(y) \right)
```

with ``G(x,y) = \frac{e^{-γ |x-y|}}{4 π |x-y|}``
"""
struct HH3DHyperSingularFDBIO{T,K} <: Helmholtz3DOp
    "coefficient of the weakly singular term"
    alpha::T
    "coefficient of the hyper singular term"
    beta::T
    "`im*κ` with `κ` the wave number"
    gamma::K
end

HH3DHyperSingularFDBIO(gamma) = HH3DHyperSingularFDBIO(gamma^2, one(gamma), gamma)
scalartype(op::HH3DHyperSingularFDBIO) = promote_type(typeof(op.alpha), typeof(op.beta), typeof(op.gamma))

"""
```math
a(u,v) = α ∬_{Γ×Γ} u(x) G_{γ}(|x-y|) v(y)
```

with ``G_{γ}(r) = \frac{e^{-γr}}{4πr}``.
"""
struct HH3DSingleLayerFDBIO{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

struct HH3DSingleLayerReg{T,K} <: Helmholtz3DOpReg
    alpha::T
    gamma::K
end

struct HH3DSingleLayerSng{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

struct HH3DDoubleLayer{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

struct HH3DDoubleLayerTransposed{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

function quaddata(op::Helmholtz3DOp, test_refspace::LagrangeRefSpace,
        trial_refspace::LagrangeRefSpace, test_elements, trial_elements)

    test_eval(x)  = test_refspace(x,  Val{:withcurl})
    trial_eval(x) = trial_refspace(x, Val{:withcurl})

    # The combinations of rules (6,7) and (5,7 are) BAAAADDDD
    # they result in many near singularity evaluations with any
    # resemblence of accuracy going down the drain! Simply don't!
    # (same for (5,7) btw...).
    # test_qp = quadpoints(test_eval,  test_elements,  (6,))
    # bssi_qp = quadpoints(trial_eval, trial_elements, (7,))

    test_qp = quadpoints(test_eval,  test_elements,  (4,))
    bsis_qp = quadpoints(trial_eval, trial_elements, (7,))

    return test_qp, bsis_qp
end

function quadrule(op::HH3DSingleLayerFDBIO, test_refspace, trial_refspace, i,
        test_element, j, trial_element, quadrature_data)

    tol, hits = sqrt(eps(eltype(eltype(test_element.vertices)))), 0
    for t in test_element.vertices
        for s in trial_element.vertices
            norm(t-s) < tol && (hits +=1; break)
    end end

    test_quadpoints  = quadrature_data[1]
    trial_quadpoints = quadrature_data[2]

    hits != 0 && return WiltonSEStrategy(
        test_quadpoints[1,i],
        DoubleQuadStrategy(
            test_quadpoints[1,i],
            trial_quadpoints[1,j]))

    return DoubleQuadStrategy(
        quadrature_data[1][1,i],
        quadrature_data[2][1,j])
end


function quadrule(op::Helmholtz3DOp, test_refspace, trial_refspace, i,
        test_element, j, trial_element, quadrature_data)

    test_quadpoints  = quadrature_data[1]
    trial_quadpoints = quadrature_data[2]

    return DoubleQuadStrategy(
        quadrature_data[1][1,i],
        quadrature_data[2][1,j])
end




function integrand(op::HH3DHyperSingularFDBIO,
        kernel, test_values, test_element, trial_values, trial_element)

    α = op.alpha
    β = op.beta

    G = kernel.green

    g, curlg = test_values
    f, curlf = trial_values

    nx = normal(test_element)
    ny = normal(trial_element)

    α*dot(nx,ny)*g*f*G + β*dot(curlg,curlf)*G
end





HH3DSingleLayerFDBIO(gamma) = HH3DSingleLayerFDBIO(one(gamma), gamma)

regularpart(op::HH3DSingleLayerFDBIO) = HH3DSingleLayerReg(op.alpha, op.gamma)
singularpart(op::HH3DSingleLayerFDBIO) = HH3DSingleLayerSng(op.alpha, op.gamma)

function integrand(op::Union{HH3DSingleLayerFDBIO,HH3DSingleLayerReg},
        kernel, test_values, test_element, trial_values, trial_element)

    α = op.alpha
    G = kernel.green
    g, curlg = test_values
    f, curlf = trial_values

    α*g*f*G
end


function innerintegrals!(op::HH3DSingleLayerSng,
        test_neighborhood, test_refspace, trial_refspace, test_elements,
        trial_element, zlocal, quadrature_rule::WiltonSEStrategy, dx)

    γ = op.gamma
    α = op.alpha

    s1, s2, s3 = trial_element.vertices

    x = cartesian(test_neighborhood)
    n = normalize((s1-s3)×(s2-s3))
    ρ = x - dot(x - s1, n) * n

    scal, vec = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
    ∫G = (scal[2] - γ*scal[3] + 0.5*γ^2*scal[4]) / (4π)

    zlocal[1,1] += α * ∫G * dx
    return nothing
end



function integrand(biop::HH3DDoubleLayer,
        kernel, fp, mp, fq, mq)

    nq = normal(mq)
    fp[1] * dot(nq, -kernel.gradgreen) * fq[1]
end



function integrand(biop::HH3DDoubleLayerTransposed,
        kernel, fp, mp, fq, mq)

    np = normal(mp)
    fp[1] * dot(np, kernel.gradgreen) * fq[1]
end


module Helmholtz3D
    using ..BEAST

    singlelayer(;
        gamma=error("propagation constant is a required argument"),
        alpha=one(gamma)) = BEAST.HH3DSingleLayerFDBIO(alpha,gamma)

    hypersingular(;
        gamma=error("propagation constant is a required argument"),
        alpha=gamma^2,
        beta=one(gamma)) = BEAST.HH3DHyperSingularFDBIO(alpha, beta, gamma)


    planewave(;
        direction=error("direction is a required argument"),
        wavenumber=error("wavenumber is a required arguement"),
        amplitude=one(eltype(direction))) = BEAST.HH3DPlaneWave(direction, wavenumber, amplitude)

    doublelayer(;gamma=error("gamma missing"), alpha=one(gamma)) =
        BEAST.HH3DDoubleLayer(alpha, gamma)

    doublelayer_transposed(;gamma=error("gamma missing"), alpha=one(gamma)) =
        BEAST.HH3DDoubleLayerTransposed(alpha, gamma)
end

export Helmholtz3D
