Feat: disc emissivity

[fjebaker]

There's a mistake in the emissivity profile as it is currently being calculated; it should be

$$ \varepsilon(r) = \frac{N(r)}{g^2 A(r, \text{d}r)}, $$

but we currently are missing the divide by redshift, where here the redshift if for the disc element and X-ray source respectively

$$ g^{-1} = \frac{E\text{disc}}{E\text{source}}. $$

This is a little annoying to implement as currently the lag transfer functions don't track the model, but they probably should.

[fjebaker]

Ahh, we're also missing the term accounting for proper area when dividing the counts to get $\varepsilon$. The area should be

$$ \text{d}A = \sqrt{g{rr} g{\theta\theta}} \ \text{d} r \ \text{d} \theta, $$

or equivalently for static, axis-symmetric:

$$ A = 2\pi \sqrt{g{rr} g{\theta\theta}} \ \text{d} r. $$

[fjebaker]

From Wilkins and Fabian (2014):

Ours at fairly low resolution:

using Gradus
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))

model = LampPostModel(h = 10.0)
# plane doesn't matter, just need something to get the transfer function
# this will probably be separated in the future so it's easier to create
# emissivity profiles without having to go through the lagtransfer function
plane = PolarPlane(GeometricGrid(); Nr = 300, Nθ = 300)

tf = @time lagtransfer(
    model,
    m,
    u,
    plane,
    d,
    callback = domain_upper_hemisphere(),
    n_samples = 200_000,
    verbose = true,
)

profile = Gradus.RadialDiscProfile(tf)

# preview the emissivity
begin
    radii = range(Gradus.isco(m), 300.0, 2000)
    us = [(; u2 = SVector(0.0, r, π / 2, 0.0)) for r in radii]
    emiss = profile.f.(us)
    # normalize maximum to 1
    emiss = emiss ./ maximum(emiss)
    p = plot(
        radii,
        emiss,
        xscale = :log10,
        yscale = :log10,
        minorgrid = true,
        label = "ε",
        xlabel = "r",
        ylabel = "emissivity",
    )
    # renormalise for comparison
    cmp = radii .^ (-3)
    cmp = cmp ./ maximum(cmp)
    plot!(p, radii, cmp, label = "r⁻³")
    p
end

[codecov-commenter]

Codecov Report

Merging #91 (b74029d) into main (d763e67) will decrease coverage by 1.21%. The diff coverage is 26.66%.

:mega: This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@            Coverage Diff             @@
##             main      #91      +/-   ##
==========================================
- Coverage   72.10%   70.89%   -1.21%     
==========================================
  Files          45       46       +1     
  Lines        1606     1639      +33     
==========================================
+ Hits         1158     1162       +4     
- Misses        448      477      +29     

Impacted Files	Coverage Δ
src/corona-to-disc/corona-models.jl	50.00% <0.00%> (-27.78%)	:arrow_down:
...ransfer-functions/cunningham-transfer-functions.jl	97.43% <ø> (ø)
src/corona-to-disc/disc-profiles.jl	58.22% <24.00%> (-7.24%)	:arrow_down:
src/transfer-functions/types.jl	50.00% <50.00%> (ø)
src/transfer-functions/transfer-functions-2d.jl	68.75% <75.00%> (+2.75%)	:arrow_up:

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[fjebaker]

Using @phajy's patent-pending overlay method for estimating goodness of fit with unknown data:

[phajy]

Oh, that looks great!