Merge pull request #177 from astro-group-bristol/fergus/gotcha-fix

fix: gotcha with thin disc lineprofiles

src/line-profiles.jl

@@ -1,3 +1,21 @@
+_warn_disc_integration_limits(::AbstractAccretionGeometry, _, _) = nothing
+function _warn_disc_integration_limits(d::ThinDisc, minrₑ, maxrₑ)
+    if !(d.inner_radius ≈ minrₑ) && (d.inner_radius != 0)
+        @warn """
+        `ThinDisc` inner radius does not equal `minrₑ`. Integration occurs only over 
+        `minrₑ`, and inner radius of the disc is ignored. Set the inner radius to `0` to
+        supress this message.
+        """
+    end
+    if !(d.inner_radius ≈ maxrₑ) && (d.outer_radius != Inf)
+        @warn """
+        `ThinDisc` outer radius does not equal `maxrₑ`. Integration occurs only over 
+        `maxrₑ`, and outer radius of the disc is ignored. Set the outer radius to `Inf` to
+        supress this message.
+        """
+    end
+end
+
 @inline function lineprofile(
     m::AbstractMetric,
     u,
@@ -48,6 +66,7 @@ function lineprofile(
     Nr = 1000,
     kwargs...,
 )
+    _warn_disc_integration_limits(d, minrₑ, maxrₑ)
     radii = Grids._inverse_grid(minrₑ, maxrₑ, numrₑ)
     itfs = interpolated_transfer_branches(m, u, d, radii; verbose = verbose, kwargs...)
     flux = integrate_lineprofile(ε, itfs, radii, bins; h = h, Nr = Nr)

src/transfer-functions/cunningham-transfer-functions.jl

@@ -1,3 +1,9 @@
+_promote_disc_for_transfer_functions(d::AbstractAccretionDisc) = (d, d)
+function _promote_disc_for_transfer_functions(::ThinDisc{T}) where {T}
+    plane = DatumPlane(zero(T))
+    plane, plane
+end
+
 struct _TransferFunctionSetup{T}
     h::T
     θ_offset::T
@@ -209,8 +215,16 @@ function _rear_workhorse(
     rₑ;
     kwargs...,
 )
-    workhorse, jacobian_function, _ =
-        _setup_workhorse_jacobian_with_kwargs(setup, m, x, d, rₑ; kwargs...)
+    disc, jacobian_disc = _promote_disc_for_transfer_functions(d)
+    workhorse, jacobian_function, _ = _setup_workhorse_jacobian_with_kwargs(
+        setup,
+        m,
+        x,
+        disc,
+        rₑ;
+        jacobian_disc = jacobian_disc,
+        kwargs...,
+    )
     function _disc_workhorse(θ::T)::NTuple{3,T} where {T}
         g, gp, r = workhorse(θ)
         α, β = _rθ_to_αβ(r, θ; α₀ = setup.α₀, β₀ = setup.β₀)

test/smoke-tests/cunningham-transfer-functions.jl

@@ -21,19 +21,19 @@ function measure_ctf(ctf)
 end
 
 # test for different angles
-@test measure_ctf(test_ctf1(0.998, 3, 4.0)) ≈ 0.12010092173620189 atol = 1e-4
-@test measure_ctf(test_ctf1(0.998, 35, 4.0)) ≈ 0.10321653667404884 atol = 1e-4
-@test measure_ctf(test_ctf1(0.998, 74, 4.0)) ≈ 0.053601943161397436 atol = 1e-4
-@test measure_ctf(test_ctf1(0.998, 85, 4.0)) ≈ 0.034591443363911956 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 3, 4.0)) ≈ 0.12161376109873144 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 35, 4.0)) ≈ 0.10362951554307089 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 74, 4.0)) ≈ 0.054070356518218586 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 85, 4.0)) ≈ 0.034811715215212875 atol = 1e-4
 
 # different radii
-@test measure_ctf(test_ctf1(0.998, 30, 4.0)) ≈ 0.1064780591310143 atol = 1e-4
-@test measure_ctf(test_ctf1(0.998, 30, 7.0)) ≈ 0.11934189586976277 atol = 1e-4
-@test measure_ctf(test_ctf1(0.998, 30, 10.0)) ≈ 0.12397743562664382 atol = 1e-4
-@test measure_ctf(test_ctf1(0.998, 30, 15.0)) ≈ 0.12741349055786652 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 30, 4.0)) ≈ 0.10779115390995794 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 30, 7.0)) ≈ 0.1202759989850966 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 30, 10.0)) ≈ 0.12461894214061674 atol = 1e-4
+@test measure_ctf(test_ctf1(0.998, 30, 15.0)) ≈ 0.1275864358885266 atol = 1e-4
 
 # large radii
-@test measure_ctf(test_ctf1(0.998, 30, 300.0)) ≈ 0.1336262903064427 rtol = 1e-2
+@test measure_ctf(test_ctf1(0.998, 30, 300.0)) ≈ 0.13191798015557799 rtol = 1e-2
 @test measure_ctf(test_ctf1(0.998, 30, 800.0)) ≈ 0.13481966054586172 rtol = 1e-2
 @test measure_ctf(test_ctf1(0.998, 30, 1000.0)) ≈ 0.13570698595034328 rtol = 1e-2
 

test/smoke-tests/reverberation.jl

@@ -39,8 +39,8 @@ model1 = LampPostModel()
 freq1, τ1 = calc_lag_freq(m, d, model1, radii, itb)
 
 # smoke test
-@test sum(freq1) ≈ 2449.8787687490535 atol = 1e-2
-@test τ1[132] ≈ 9.281677930459137 atol = 1e-2
+@test sum(freq1) ≈ 2449.8787687490535 rtol = 1e-2
+@test τ1[132] ≈ 9.281677930459137 rtol = 1e-2
 
 # check it works for thick discs too
 d_thick = ShakuraSunyaev(m)
@@ -49,4 +49,4 @@ itb_thick = Gradus.interpolated_transfer_branches(m, x, d_thick, radii; β₀ =
 freq2, τ2 = calc_lag_freq(m, d_thick, model1, radii, itb_thick)
 
 # should be the same at this inclination
-@test sum(freq2) ≈ 2449.8787687490535 atol = 1e-2
+@test sum(freq2) ≈ 2449.8787687490535 rtol = 1e-2