Q.jl

module Q

using Cubature
using QuadGK
using GSL

function fi(ξ::Real, α::Real, E::Number, k4::Real, kp::Real)
    p1 = (k4^2 + E/2)^2 + kp^2 + α^2
    p2 = (k4^2 - E/2)^2 + kp^2 + (1-α)^2
    2 * 2kp *
        (α * (1+p1)^1.5ξ * (1+p2)^0.5ξ) /
        ((1+p1*(1+p1)^2ξ) * (1+p2*(1+p2)^2ξ))
end

function f(ξ::Real, α::Real, E::Real)
    function ii(x)
        k4, kp = x[1]/(1-x[1]), x[2]/(1-x[2])
        # @show x, k4, kp
        fi(ξ, α, E, k4, kp) / ((1-x[1]) * (1-x[2]))^2
    end
    val, err = hcubature(ii, (0.0, 0.0), (1.0, 1.0), reltol=0, abstol=1e-8, maxevals=10^6)
    # @show val, err
    val
end

function f(ξ::Real, α::Real, E::Complex)
    function ii(x, out)
        k4, kp = x[1]/(1-x[1]), x[2]/(1-x[2])
        # @show x, k4, kp
        w = fi(ξ, α, E, k4, kp) / ((1-x[1]) * (1-x[2]))^2
        out[1] = real(w)
        out[2] = imag(w)
    end
    val, err = hcubature(2, ii, (0.0, 0.0), (1.0, 1.0), reltol=0, abstol=1e-8, maxevals=10^6)
    # @show val, err
    val[1] + val[2] * im
end

function g1(ξ::Real, α::Real, t::Real, e0::Real)
    const θ = 0
    γ, δ = cis(θ), cis(π - θ)
    quadgk(u ->
           γ * cis(t * γ * u) / ((γ*u)^2 + e0^2) * f(ξ, α, γ*u) -
           δ * cis(t * δ * u) / ((δ*u)^2 + e0^2) * f(ξ, α, δ*u),
           0, Inf)
end

function gslfun(x::Cdouble, params::Ptr{Void})
    (unsafe_pointer_to_objref(params)::Function)(x)::Cdouble
end
gslfun_c = cfunction(gslfun, Cdouble, (Cdouble, Ptr{Void}))
Base.convert(::Type{GSL.gsl_function}, f::Function) = GSL.gsl_function(gslfun_c, pointer_from_objref(f))

function g2(ξ::Real, α::Real, t::Real, e0::Real)
    limit = 1024
    ws1 = GSL.integration_workspace_alloc(limit)
    ws2 = GSL.integration_workspace_alloc(1024)

    qawo_cos = ccall((:gsl_integration_qawo_table_alloc, GSL.libgsl),
                     Ptr{GSL.gsl_integration_qawo_table},
                     (Cdouble, Cdouble, Cint, Csize_t),
                     t, 1.0, 0, 10)
    qawo_sin = ccall((:gsl_integration_qawo_table_alloc, GSL.libgsl),
                     Ptr{GSL.gsl_integration_qawo_table},
                     (Cdouble, Cdouble, Cint, Csize_t),
                     t, 1.0, 1, 10)

    f1 = E -> (f(ξ, α, E) + f(ξ, α, -E))/(E^2+e0^2)
    f2 = E -> (f(ξ, α, E) - f(ξ, α, -E))/(E^2+e0^2)
    result1, result2, abserr = Ref(NaN), Ref(NaN), Ref(NaN)
    ret = ccall((:gsl_integration_qawf, GSL.libgsl), Cint,
                (Ref{gsl_function}, Cdouble, Cdouble, Csize_t, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_qawo_table}, Ref{Cdouble}, Ref{Cdouble}),
                f1, 0.0, 1e-8, limit, ws1, ws2, qawo_cos, result1, abserr)
    @show ret
    @show result1, abserr
    ret == 0 || warn("FAILED: gsl_integration_qawf")
    ret = ccall((:gsl_integration_qawf, GSL.libgsl), Cint,
                (Ref{gsl_function}, Cdouble, Cdouble, Csize_t, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_workspace}, Ptr{GSL.gsl_integration_qawo_table}, Ref{Cdouble}, Ref{Cdouble}),
                f2, 0.0, 1e-8, limit, ws1, ws2, qawo_sin, result2, abserr)
    @show ret
    @show result2, abserr
    ret == 0 || warn("FAILED: gsl_integration_qawf")

    result1[] + im * result2[]
end

end