mwe.jl

module mwe

using LinearAlgebra: dot, norm
using Statistics: mean
using Cthulhu
using InteractiveUtils

cmplx(r, c) = r + im * c
size3(v) = norm(v)
ssize3(v) = v[1]^2 + v[2]^2 + v[3]^2
dot3(v, w) = dot(v, w)


function off_triangle_rule(ir)
    p, w = ir.param_coords, vec(ir.weights)
    xqoff = vec(p[:, 1])
    yqoff = vec(p[:, 2])
    wqoff = deepcopy(w) .*= 2
    return xqoff, yqoff, wqoff
end

function on_triangle_rule(ir)
    xqoff, yqoff, wqoff = off_triangle_rule(ir)
    nqoff = length(wqoff)
    nqon=3*nqoff
    z = 0.0
    xqon = fill(z, nqon)
    yqon = fill(z, nqon)
    wqon = fill(z, nqon)
    for iq in 1:nqoff
        xqon[iq] = xqoff[iq]*(1/3)+yqoff[iq]
        yqon[iq] = xqoff[iq]*(1/3)
        wqon[iq] = wqoff[iq]/3
        xqon[iq+nqoff] = xqoff[iq]*(1/3)
        yqon[iq+nqoff] = xqoff[iq]*(1/3)+yqoff[iq]
        wqon[iq+nqoff] = wqoff[iq]/3
        xqon[iq+2*nqoff] = (1/3)*(1+2*xqoff[iq]-yqoff[iq])
        yqon[iq+2*nqoff] = (1/3)*(1-xqoff[iq]+2*yqoff[iq])
        wqon[iq+2*nqoff] = wqoff[iq]/3
    end
    return xqon, yqon, wqon
end

function area(qa, qb, qc)
    a1, a2, a3 = qb[1] - qa[1], qb[2] - qa[2], qb[3] - qa[3]
    b1, b2, b3 = qc[1] - qa[1], qc[2] - qa[2], qc[3] - qa[3]
    n1 = a2*b3-a3*b2
    n2 = a3*b1-a1*b3
    n3 = a1*b2-a2*b1
    return sqrt(n1^2 + n2^2 + n3^2) / 2
end

function f(
    k,
    p,
    vecp,
    qa,
    qb,
    qc,
    lponel,
    xq,
    yq,
    wq,
    llk,
    lmk,
    lmkt,
    lnk,
    temps
)
    pa, normq, ar0, ara, aopp, q, rr, qbmqa, qcmqa, qcmqb, pmqa, pmqb, pmqc  = temps

    zero = 0.0
    one = 1.0
    two = 2.0
    three = 3.0
    four = 4.0
    third = 1.0 / 3.0
    fourpi = four * pi
    czero = cmplx(zero, zero)
    cone = cmplx(one, zero)
    cimag = cmplx(zero, one)
    ek = 1.0e-16
    egeom = 0.0
    lvalid = false

    @. qbmqa = qb - qa
    @. qcmqa = qc - qa
    @. qcmqb = qc - qb
    @. pmqa = p - qa
    @. pmqb = p - qb
    @. pmqc = p - qc
    qarea = area(qa, qb, qc)

    lfail = false
    lerror = false

    llorn = llk || lnk
    lmorn = lmk || lnk
    lmsorn = lmk || lmkt || lnk
    lmtorn = lmkt || lnk
    lmormt = lmk || lmkt
    if (lponel)
        lmksv = lmk
        lmktsv = lmkt
        lmk = false
        lmkt = false
    end

    if (lmorn)
        normq[1] = qbmqa[2] * qcmqa[3] - qbmqa[3] * qcmqa[2]
        normq[2] = qcmqa[1] * qbmqa[3] - qbmqa[1] * qcmqa[3]
        normq[3] = qbmqa[1] * qcmqa[2] - qbmqa[2] * qcmqa[1]
        acrss = size3(normq)
        normq ./= acrss
    end

    if (lnk)
        dnpnq = dot3(vecp, normq)
    end

    imk = imag(k)
    rek = real(k)
    sk = k * k
    sko2 = sk / two
    resko2 = real(sko2)
    kzero = false
    kreal = false
    kimag = false
    kcomp = false

    if (abs(rek) < ek && abs(imk) < ek)
        kzero = true
    elseif (abs(imk) < ek)
        kreal = true
    elseif (abs(rek) < ek)
        kimag = true
    else
        kcomp = true
    end

    csuml = czero
    csumm = czero
    csummt = czero
    csumn = czero
    if (lponel && llorn)
        rqbqc = size3(qcmqb)
        rqcqa = size3(qcmqa)
        rqaqb = size3(qbmqa)
        rqap = size3(pmqa)
        rqbp = size3(pmqb)
        rqcp = size3(pmqc)
        ar0 .= (rqap, rqbp, rqcp)
        ara .= (rqbp, rqcp, rqap)
        aopp .= (rqaqb, rqbqc, rqcqa)
        @inbounds for ii = 1:3
            r0 = ar0[ii]
            ra = ara[ii]
            opp = aopp[ii]
            if (r0 < ra)
                temp = ra
                ra = r0
                r0 = temp
            end
            sr0 = r0 * r0
            sra = ra * ra
            sopp = opp * opp
            a = acos((sra + sr0 - sopp) / two / ra / r0)
            b = atan(ra * sin(a) / (r0 - ra * cos(a)))
            csuml += cmplx(
                    (r0 * sin(b) * (log(tan((b + a) / two)) - log(tan(b / two)))) / qarea,
                    zero,
                )
            csumn += cmplx((cos(b + a) - cos(b)) / r0 / sin(b) / qarea, zero)
        end
    end
    csumn = csumn - sko2 * csuml

    if (lponel && kzero)
        return czero
    end

    nq = length(xq)
    for iq = 1:nq
        wqiq = wq[iq]
        @. q = qa + xq[iq] * qbmqa + yq[iq] * qcmqa
        @. rr = p - q
        sr = ssize3(rr)
        r = sqrt(sr)
        if (lnk)
            cr = r * sr
        end
        if (lmorn)
            rnq = -dot3(rr, normq) / r
        end
        if (lmtorn)
            rnp = dot3(rr, vecp) / r
        end
        if (lnk)
            rnprnq = rnp * rnq
        end
        if (lnk)
            rnpnq = -(dnpnq + rnprnq) / r
        end
        if (lponel && (!kzero))
            if (llorn)
                fpg0 = one / r
            end
            if (lmsorn)
                fpg0r = -one / sr
            end
            if (lnk)
                fpg0rr = two / cr
            end
        end

        fpg, fpgr, wfpgr = czero, czero, czero
        if (kzero)
            if (llk)
                fpg = cmplx(one / r, zero)
            end
            if (llk)
                csuml += cmplx(wqiq * real(fpg), zero)
            end
            if (lmsorn)
                fpgr = cmplx(-one / sr, zero)
            end
            if (lmormt)
                wfpgr = cmplx(wqiq * real(fpgr), zero)
            end
            if (lmk)
                csumm += cmplx(real(wfpgr) * rnq, zero)
            end
            if (lmkt)
                csummt += cmplx(real(wfpgr) * rnp, zero)
            end
            if (lnk)
                fpgrr = cmplx(two / cr, zero)
                csumn += cmplx(wqiq * (real(fpgr) * rnpnq + real(fpgrr) * rnprnq), zero)
            end
        elseif (kreal)
            kr = cmplx(real(k) * r, zero)
            ikr = cmplx(zero, real(kr))
            skr = cmplx(real(kr) * real(kr), zero)
            e = exp(ikr)
            if (llk)
                fpg = e / r
            end
            if (llk)
                if (!lponel)
                    csuml += wqiq * fpg
                else
                    csuml += wqiq * (fpg - fpg0)
                end
            end
            if (lmsorn)
                eosr = e / sr
                fpgr = im*(eosr) * real(kr) - eosr
            end
            if ((!lponel) && lmormt)
                wfpgr = wqiq * fpgr
                if (lmk)
                    csumm += wfpgr * rnq
                end
                if (lmkt)
                    csummt += + wfpgr * rnp
                end
            end
            if (lnk)
                fpgrr = e * (two - ikr - ikr - skr) / cr
                if (!lponel)
                    csumn += wqiq * (fpgr * rnpnq + fpgrr * rnprnq)
                else
                    csumn += wqiq *
                        ((fpgr - fpg0r) * rnpnq + (fpgrr - fpg0rr) * rnprnq + resko2 * fpg0)
                end
            end
        elseif (kimag)
            imkr = imk * r
            rekr = zero
            kr = cmplx(rekr, imkr)
            ikr = cmplx(-imkr, rekr)
            skr = cmplx(-imkr * imkr, zero)
            e = exp(ikr)
            if (llk)
                fpg = cmplx(real(e) / r, zero)
            end
            if (llk)
                if (!lponel)
                    csuml += cmplx(wqiq * real(fpg), zero)
                end
                if (lponel)
                    csuml += cmplx(wqiq * real(fpg - fpg0), zero)
                end
            end
            if (lmsorn)
                eosr = cmplx(real(e) / sr, zero)
                fpgr = cmplx(real(eosr) * real(ikr - one), zero)
            end
            if ((!lponel) && lmormt)
                wfpgr = cmplx(wqiq * real(fpgr), zero)
                if (lmk)
                    csumm += cmplx(real(wfpgr) * rnq, zero)
                end
                if (lmkt)
                    csummt += cmplx(real(wfpgr) * rnp, zero)
                end
            end
            if (lnk)
                fpgrr = cmplx(real(e) * real(two - ikr - ikr - skr) / cr, zero)
                if (!lponel)
                    csumn += cmplx(wqiq * (real(fpgr) * rnpnq + real(fpgrr) * rnprnq), zero)
                else
                    csumn += cmplx(
                            wqiq * (
                                real(fpgr - fpg0r) * rnpnq +
                                real(fpgrr - fpg0rr) * rnprnq +
                                resko2 * fpg0
                            ),
                            zero,
                        )
                end
            end
        elseif (kcomp)
            kr = k * r
            ikr = im * (kr)
            skr = sk * sr
            e = exp(ikr)
            if (llk)
                fpg = e / r
            end
            if (llk)
                if (!lponel)
                    csuml += wqiq * fpgv
                end
                if (lponel)
                    csuml += wqiq * (fpg - fpg0)
                end
            end
            if (lmsorn)
                eosr = e / sr
                fpgr = im * (eosr) * kr - eosr
            end
            if ((!lponel) && lmormt)
                wfpgr = wqiq * fpgr
                if (lmk)
                    csumm += wfpgr * rnq
                end
                if (lmkt)
                    csummt += wfpgr * rnp
                end
            end
            if (lnk)
                fpgrr = e * (two - ikr - ikr - skr) / cr
                if (!lponel)
                    csumn +=  wqiq * (fpgr * rnpnq + fpgrr * rnprnq)
                else
                    csumn += wqiq *
                        ((fpgr - fpg0r) * rnpnq + (fpgrr - fpg0rr) * rnprnq + sko2 * fpg0)
                end
            end
        end

    end

    qao4pi = qarea / fourpi

    if llk
        return dislk = qao4pi * csuml
    elseif lmk
        return dismk = qao4pi * csumm
    elseif lmkt
        return dismkt = qao4pi * csummt
    else
        return disnk = qao4pi * csumn
    end

end

k = 1.0 + 0.0 * im
p = fill(0.0, 3)
vecp = fill(0.0, 3); vecp[3] = 1
qa = fill(0.0, 3);
qb = fill(0.0, 3); qb[1] = 1
qc = fill(0.0, 3); qc[2] = 1
lponel = false
xq = fill(0.0, 3)
yq = fill(0.0, 3)
wq = fill(0.0, 3)
llk = true
lmk = false
lmkt = false
lnk = false
temps = fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3),
fill(zero(Float64), 3)

@descend f(
    k,
    p,
    vecp,
    qa,
    qb,
    qc,
    lponel,
    xq,
    yq,
    wq,
    llk,
    lmk,
    lmkt,
    lnk,
    temps
    )

end