Debugging coronal beaming emissivities

[fjebaker]

This PR is primarily in relation to trying to work out what's going on with the emissivity profiles when we have a moving source. @wiktoriatarnopolska has added the model from the literature (e.g. Wilkins+12, Dauser+13, Gonzalez+17; PR pending):

# structing a moving source
struct MovingSource{T} <: Gradus.AbstractCoronaModel{T}
    r::T
    β::T
end

# helper functions
drdt(g, β) = β * √(-g[1] / g[2])
drdt(m::AbstractMetric, x, β) = drdt(metric_components(m, SVector(x[2], x[3])), β)
function Gradus.sample_position_velocity(m::AbstractMetric, model::MovingSource)
    x = SVector{4}(0, model.r, 1e-4, 0)
    v = SVector(1, drdt(m, x, model.β), 0, 0)
    x, v
end

# blazing fast, exploiting symmetry
Gradus.emissivity_profile(
    ::Nothing,
    m::AbstractMetric,
    d::AbstractAccretionDisc,
    model::MovingSource,
    spec::Gradus.AbstractCoronalSpectrum;
    kwargs...,
) = Gradus._point_source_symmetric_emissivity_profile(m, d, model, spec; kwargs...)

This PR adds utility functions and small little changes that attempts to help get some kind of consensus between the results published in the literature. I will use the comments in this PR to document progress.

Relevant figures from the literature:

Wilkins+12 Figure 11:

Dauser+13 Figure 6:

Gonzalez+17 Fig 6:

[codecov-commenter]

Codecov Report

Merging #148 (01c0159) into main (50015ff) will decrease coverage by 0.74%.
Report is 21 commits behind head on main.
The diff coverage is 60.43%.

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

@@            Coverage Diff             @@
##             main     #148      +/-   ##
==========================================
- Coverage   68.81%   68.07%   -0.74%     
==========================================
  Files          58       59       +1     
  Lines        2559     2622      +63     
==========================================
+ Hits         1761     1785      +24     
- Misses        798      837      +39     

Files Changed	Coverage Δ
src/Gradus.jl	26.08% <0.00%> (-6.53%)	⬇️
src/metrics/minkowski.jl	0.00% <0.00%> (ø)
src/plotting-recipes.jl	0.00% <0.00%> (ø)
src/tracing/method-implementations/first-order.jl	70.21% <0.00%> (ø)
src/tracing/constraints.jl	72.72% <40.00%> (-27.28%)	⬇️
src/corona/flux-calculations.jl	32.55% <75.00%> (-4.66%)	⬇️
src/orbits/circular-orbits.jl	73.61% <84.61%> (+1.94%)	⬆️
src/corona/profiles/radial.jl	90.90% <85.00%> (+5.45%)	⬆️
src/corona/emissivity.jl	93.75% <100.00%> (+0.13%)	⬆️
src/corona/models/lamp-post.jl	100.00% <100.00%> (ø)
... and 6 more

... and 1 file with indirect coverage changes

[fjebaker]

Ostensibly, all of these authors are using the same beaming model, where the source velocity is described as

$$ v_\text{source}^\mu = \left( 1, \beta \frac{\text{d} r }{\text{d} t}, 0, 0 \right), $$

where $\beta$ is the coronal velocity.

Gonzalez+17 derive $\text{d} r / \text{d} t$ from $\text{d}s^2 = 0$ for $\text{d} \theta = \text{d} \phi = 0$, and therefore

$$ \frac{\text{d} r}{\text{d} t} = \sqrt{\frac{g_{tt}}{g_{rr}}}. $$

This $\text{d} r / \text{d} t$ is the coordinate velocity in a ZAMO frame co-located with the corona. A tetrad basis is then constructed from this velocity, with the seed vector

$$ \mathbf{e}^\mu_{(t)} = \sqrt{\frac{1}{g_{tt} + v^2 g_{rr}}} \left( 1, v, 0, 0 \right), $$

using the shorthand $v = \beta \text{d} r / \text{d} t$.

... I just noticed we are missing that normalization in the velocity prescription of the MovingSource model.

[fjebaker]

With the normalization on the velocity, i.e.

function Gradus.sample_position_velocity(m::AbstractMetric, model::MovingSource)
    x = SVector{4}(0, model.r, 1e-4, 0)
    g = metric_components(m, SVector(x[2], x[3]))
    u = drdt(g, model.β)
    v = √(inv(-g[1] - u^2 * g[2])) .* SVector(1, u, 0, 0)
    x, v
end

Without:

Clearly the normalization corrects for the difference between the $\beta = 0$ and LP case, which is nice 👍

[fjebaker]

Note the fact that ours do not quite asymptotically go to $\alpha = 3$ is because I am only using a few thousand samples to keep computation time really short whilst I debug. Increasing samples fixes the large $r$ behaviour.

[fjebaker]

Recreating Fig 8 from Fukumura & Kazanas 2007:

Ours:

I had to guess the number of geodesics but it turns out it is exactly 37.

[fjebaker]

• The difference the velocity normalization makes:

• Vs not normalizing: