Lineprofiles: root finder is temperamental when observer close to equitorial plane

[fjebaker]

As described in the comments of

• 71

when the observer angle is close to 85, the root finder misses the disc underneath the equitorial plane.

Code to reproduce (v0.2.3):

using Gradus, StaticArrays, Plots
m = KerrMetric(M=1.0, a=0.998)
d = GeometricThinDisc(0.0, 3000.0, π/2)

angle = 85

u = @SVector [0.0, 1000.0, deg2rad(angle), 0.0]

# errors here
ctf = cunningham_transfer_function(m, u, d, 4.0)

mask = @. (ctf.g✶ > 0.001) & (ctf.g✶ < 0.999) 
plot(ctf.g✶[mask], ctf.f[mask], label="θᵢ=$angle", xlabel="g✶", ylabel="f") 

Setting e.g. offset_max = 15.0 seems to fix it, but produces visible errors in the resulting transfer function:

[fjebaker]

using Gradus
using StaticArrays

m = KerrMetric(M=1.0, a=0.998)
u = @SVector [0.0, 1000.0, deg2rad(80), 0.0]
d = GeometricThinDisc(Gradus.isco(m), 150.0, π/2)

r_target = 4.0
α, β = impact_parameters_for_radius(
    m,
    u,
    d,
    r_target;
    offset_max = 20.0,
    N = 100,
)

# visualise projection onto observer's plane
plot(α, β) |> display

# plot 3d
vs = [map_impact_parameters(m, u, a, b) for (a, b) in zip(α, β)]
us = fill(u, length(vs))
# retrace for full paths
sols = tracegeodesics(m, us, vs, d, (0.0, 2000.0))

function plot_event_horizon!(p, m; N = 32, kwargs...)
    R, ϕ = event_horizon(m, resolution = N)

    θ = range(0, stop=π, length=N)
    x = R .* cos.(ϕ) .* sin.(θ)'
    y = R .* sin.(ϕ) .* sin.(θ)'
    z = R .* repeat(cos.(θ)',outer=[N, 1])

    surface!(p, x, y, z; kwargs...)
end

begin
    LIM = 8
    p = plot(legend=false)
    for sol in sols.u
        tspan = range(950, Base.min(1050, sol.t[end]), 300)
        # interpolate at selected times
        points = sol.(tspan)
        # unpack 
        r = [i[2] for i in points]
        θ = [i[3] for i in points]
        ϕ = [i[4] for i in points]
        # transform to cartesian
        x = @. r * sin(θ) * cos(ϕ)
        y = @. r * sin(θ) * sin(ϕ)
        z = @. r * cos(θ)
        # plot
        plot!(p, x, y, z)
    end
    plot!(p, [0], [0], [0], ms=50)

    plot_event_horizon!(p, m)
    p = plot(p, xlims = (-LIM, LIM), ylims = (-LIM, LIM), zlims = (-LIM ÷ 2, LIM), legend = false, camera=(20, 10))

    p
end

For 80 degrees:

For 85 degrees it works fine too, but the Cunningham TF fails when it tries to calculate the Jacobian, presumably then it is slipping under the disc.

[fjebaker]

Bumping N=120 in cunningham_transfer_function fixes the jump:

[fjebaker]

This is most likely then a bug in the interpolated estimation of the g_star extrema!