import HCubature: hquadrature as integrate
using SpecialFunctions, Plots

ωmin, ωmax = 1e-10, 1000; imsmax = 100;
tgrid = range(-1, 1, length = 666)

fω(ω, a, b) = 1/(im*b*abs(ω)^a - ω)

function fs(s, a, b, sgnt)
    a1 = sgnt > 0 ? 2 - a : a
    s1 = s/abs(a - 1)
    s2 = (1/2) * abs(a1) * s1
    gammas = gamma(s) * (gamma(s1) * gamma(1 - s1)) / (gamma(s2) * gamma(1 - s2))
    return - 2*π*im/abs(a - 1) * sign(a1) * b^(s/(a - 1)) * gammas
end

function f(t, a, b)
    fωt(ω) = exp(im*ω*t)*fω(ω, a, b)
    valp, errp = integrate(fωt,  ωmin,  ωmax)
    valm, errm = integrate(fωt, -ωmax, -ωmin)
    return valp + valm
end

function 🦊H(t, a, b)
    c = abs(a - 1)/2
    sgnt = sign(t)
    fst(ims) = (s = c + ims*im; ds = im; 1/(2*π*im) * abs(t)^(-s) * fs(s, a, b, sgnt) * ds)
    val, err = integrate(fst,  -imsmax,  imsmax)
    return val - π*im*(a > 1)*sign(t)
end

function makeplot(a, b)
    fourier = [imag(f(t, a + 1e-6, b)) for t ∈ tgrid]
    mellin = [imag(🦊H(t, a + 1e-6, b)) for t ∈ tgrid]
    plot(tgrid, [fourier mellin],
        label = ["original integral, a = $a, b = $b" "uglier integral, a = $a, b = $b"],
        framestyle = :zerolines)
end

animate_being = @animate for a ∈ -1:0.1:3 makeplot(a, 1) end
gif(animate_being, fps = 10)