2D transfer functions

[phajy]

I'm trying to run the 2D transfer function example you shared earlier, @fjebaker, but I can't get it to work because I think a few things have changed under the hood in Gradus. The two errors I get are that getgeodesicpoint not defined and radial_sampler not defined.

[fjebaker]

Yes, there is a new API (~;

using Gradus
using Plots
using StaticArrays

m = KerrMetric(M = 1.0, a = 1.0)
u = @SVector [0.0, 1e6, deg2rad(60), 0.0]
d = GeometricThinDisc(Gradus.isco(m), 500.0, deg2rad(90.0))

tf = @time lagtransfer(m, u, d, callback = domain_upper_hemisphere())

begin
    t, E, f = binflux(tf, N_E=900, N_t=900)
    # log data where not zero
    I = f .> 0
    f[I] .= log.(f[I])
    p = heatmap(
        t,
        E,
        f,
        xlabel = "Time (GM/c^3)",
        ylabel = "Energy (keV)",
        xrange = [0, 150],
        yrange = [0, 9],
        clims = (-20, -1),
    )
end

There's a good chance this won't work in main just yet, since I have some polish commits I am still finishing, but I'll link the PR for which this works when it's merged!

[phajy]

Perfect - thanks!!

One warning which I ignored:

┌ Warning: dt(2.3283064365386963e-10) <= dtmin(2.3283064365386963e-10) at t=1.7537188949127322e-5. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/hLrpl/src/integrator_interface.jl:529

I can now star experimenting with different transfer functions :)

[phajy]

P.S. The new API is a lot shorter than the old method (at least from the end user's point of view!)

[phajy]

If I wanted to "fill in" the gaps in the red wing at later times what would you recommend - increasing n_samples, Nr, or something else? These photons are from small $r$ in the disc so only take up a relatively small solid angle of the image plane at $\infty$.

[fjebaker]

As you say, the red wing at large lag is from an extremely small region deep in the potential well, which is poorly sampled. You could increase say Nr or Nθ to probe this region better in the image plane (the n_samples is for the source-to-disc which is pretty well defined in that region).

Another option is to change the image plane gridding, and use e.g. InverseGrid(), which samples r as 1/r and therefore samples this inner region more.

E.g.:

# have to now explicitly define the model and plane we want to use
model = LampPostModel(h = 10.0, θ = deg2rad(0.0001))
# bumped Nr from 800 to 1200
plane = PolarPlane(GeometricGrid(); Nr = 1200, Nθ = 800, r_max = 250.0)

# nb: the order of these parameters will probably change in upcoming versions
# to be more consistent with the (m, u, plane, d) of `tracegeodesics`, and have model
# as a kwarg
tf = @time lagtransfer(model, m, u, plane, d, callback = domain_upper_hemisphere())

Fills in the red wing a little better in t~[50-70].

plane = PolarPlane(InverseGrid(); Nr = 800, Nθ = 800, r_max = 250.0)

The tradeoff is obvious. If you can think of a better distribution of r, then please feel free to implement a new AbstractImpactParameterGrid (examples are in src/image-planes/grids.jl).

I've seen e.g. Dauser's Relline table is more sigmoid-skewed, which might be a good option too? Or some blending between inverse at low r and geometric at larger r?

[fjebaker]

@phajy can we continue the discussion in #79, so that it is easier to search through for information later on?

[phajy]

Yes - see you in #79