Gradus.jl + SpectralFitting.jl = <3

[fjebaker]

Using this issue to track what needs to get done for seamless integration of Gradus.jl with SpectralFitting.jl so that the models can easily be used and fit to data.

With #191 and #190 things are in good shape. Transfer functions for given parameters are now fast to compute (from ~20 seconds down to ~4 seconds), which means we don't need to be passing pre-computed tables around, but the user can compute their own dense table in about 10 minutes or so. The transfer function interpolation scheme is also ready, so we just need to provide some helper structures and examples to get everything working together.

Approximate things todo:

• Speed up transfer function integration #192
• SpectralFitting.jl package extension #193

The goal is to get spectral models working in a way that lets the user choose how to put the model together (whether they want different coronal models or a specific disc geometry).

Then useful to have but not pressing

• Transfer function IO #194
• Fast disc emissivities #196
• Thick disc transfer function calculations are slow and maybe erroneous #198
• Auto-diff through transfer function interpolation and integration #199

[fjebaker]

@phajy let me know if there's anything I'm missing!

[fjebaker]

Okay here's a fully working example of an additive model with a fixed analytic emissivity function:

using Gradus, Plots, SpectralFitting

struct LineProfile{D,T} <: AbstractTableModel{T,Additive}
    table::D
    K::T
    "Spin" 
    a::T
    θ::T 
    rmin::T 
    rout::T 
end

function LineProfile(
    profile,
    table::Gradus.CunninghamTransferTable;
    K = FitParam(1.0),
    a = FitParam(0.998),
    θ = FitParam(45.0),
    rmin = FitParam(1.0),
    rout = FitParam(100.0, upper_limit = 100.0)
)
    setup = integration_setup(profile, table((get_value(θ), get_value(a))))
    LineProfile((; setup = setup, table = table), K, a, θ, rmin, rout)
end

function SpectralFitting.invoke!(output, domain, model::LineProfile)
    grid = model.table.table((model.θ, model.a))
    rmin = if model.rmin < grid.r_grid[1]
        grid.r_grid[1]
    else
        model.rmin
    end
    Gradus.integrate_lineprofile!(output, model.table.setup, grid, domain; rmin = rmin, rmax = model.rout)
    output
end

function mapper(a, theta)
    @info a, theta
    m = KerrMetric(1.0, a)
    x = SVector(0.0, 10000.0, deg2rad(theta), 0.0)
    d = ThinDisc(0.0, Inf)

    radii = Gradus.Grids._inverse_grid(Gradus.isco(m) + 1e-2, 100.0, 100)
    @time Gradus.transfer_function_grid(m, x, d, radii; verbose = true)
end

# emissivity function
ε(r) = r^(-3)

function direct_computation(f, domain, a, θ; rmax = 50.0)
    x = SVector(0.0, 1000.0, deg2rad(θ), 0.0)
    d = ThinDisc(0.0, 400.0)
    m = KerrMetric(1.0, a)
    radii = Gradus.Grids._inverse_grid(Gradus.isco(m) + 1e-2, rmax, 100)
    tf_grid = @time Gradus.transfer_function_grid(m, x, d, radii; verbose = true)
    Gradus.integrate_lineprofile(f, tf_grid, domain; rmax = rmax)
end

as = collect(range(0.0, 0.998, 10))
thetas = collect(range(0.1, 89.0, 10))
# big_grid = [mapper(a, t) for a in as, t in thetas]

table = Gradus.CunninghamTransferTable((thetas, as), big_grid)

a = 0.998
θ = 59.0
rout = 50.0

model = LineProfile(ε, table)
set_value!(model.θ, θ)
set_value!(model.a, a)
set_value!(model.rout, rout)

domain = collect(range(0.1, 1.4, 100))
begin
    output = invokemodel(domain, model)
    plot(domain[1:end-1], output)
end

begin
    flux = direct_computation(ε, domain, a, θ; rmax = rout)
    plot!(domain, flux)
end

When the parameters match knots of the interpolation, the results are identical. In other regions there are errors, but that can be improved with more interpolation knots.

[fjebaker]

Low inclinations are poor performing. In the above that's because it's interpolating between 0.1 and 9.0, and at 0.1 the transfer functions are really noisy so it interpolates a negative value...

[fjebaker]

Yeah fix is simple, just don't go to 0.1 without also going to 1.0 to get a stable knot.

[fjebaker]

This has been more or less fully realized now so closing this issue.