fix!: thick discs only receive projected radius

src/geometry/discs.jl

@@ -44,8 +44,6 @@ end
 """
 cross_section(d::AbstractThickAccretionDisc, x4) =
     error("Not implemented for $(typeof(d)).")
-r_cross_section(d::AbstractThickAccretionDisc, r::Number) =
-    cross_section(d, SVector(0, r, π / 2, 0))
 
 struct EllipticalDisc{T} <: AbstractAccretionDisc{T}
     inner_radius::T

src/geometry/discs/datum-plane.jl

@@ -8,13 +8,11 @@ function distance_to_disc(d::DatumPlane{T}, x4; gtol) where {T}
     abs(h) - d.height - (gtol * x4[2])
 end
 
-function datumplane(disc::AbstractThickAccretionDisc, r::T) where {T}
-    h = r_cross_section(disc, r)
-    DatumPlane(h)
-end
-function datumplane(disc::AbstractThickAccretionDisc, x::SVector{4})
-    h = cross_section(disc, x)
+function datumplane(disc::AbstractThickAccretionDisc, ρ::T) where {T}
+    h = cross_section(disc, ρ)
     DatumPlane(h)
 end
+datumplane(disc::AbstractThickAccretionDisc, x::SVector{4}) =
+    datumplane(disc, _equatorial_project(x))
 
 export DatumPlane

src/geometry/discs/polish-doughnut.jl

@@ -113,13 +113,11 @@ function PolishDoughnut(m::AbstractMetric; rₖ = 12.0, n = 0.21, init_r = 5.0,
     PolishDoughnut(inner_radius, outer_radius, rₖ, n, f)
 end
 
-function cross_section(d::PolishDoughnut, u4)
-    let r = u4[2]
-        if d.inner_radius ≤ r ≤ d.outer_radius
-            d.f(r)
-        else
-            0.0
-        end
+function cross_section(d::PolishDoughnut, ρ)
+    if d.inner_radius ≤ ρ ≤ d.outer_radius
+        d.f(ρ)
+    else
+        zero(typeof(ρ))
     end
 end
 

src/geometry/discs/shakura-sunyaev.jl

@@ -26,12 +26,11 @@ struct ShakuraSunyaev{T} <: AbstractThickAccretionDisc{T}
     risco::T
 end
 
-@fastmath function cross_section(d::ShakuraSunyaev, x)
-    r = _equatorial_project(x)
-    if r < d.risco
-        return -one(typeof(r))
+@fastmath function cross_section(d::ShakuraSunyaev, ρ)
+    if ρ < d.risco
+        return -one(typeof(ρ))
     end
-    H = (3 / 2) * inv(d.η) * (d.Ṁ / d.Ṁedd) * (1 - sqrt(d.risco / r))
+    H = (3 / 2) * inv(d.η) * (d.Ṁ / d.Ṁedd) * (1 - sqrt(d.risco / ρ))
     2H
 end
 

src/geometry/discs/thick-disc.jl

@@ -46,42 +46,36 @@ function ThickDisc(cross_section::F, parameters::P) where {F,P}
     ThickDisc{Float64,F,P}(cross_section, parameters)
 end
 
-cross_section(d::ThickDisc, x4) = d.f(x4, d.params...)
-cross_section(d::ThickDisc{T,F,Nothing}, x4) where {T,F} = d.f(x4)
+cross_section(d::ThickDisc, ρ) = d.f(ρ, d.params...)
+cross_section(d::ThickDisc{T,F,Nothing}, ρ) where {T,F} = d.f(ρ)
 
-xz_parameterize(d::AbstractAccretionDisc, ρ) =
-    SVector(ρ, cross_section(d, SVector(0, ρ, π / 2, 0)))
+xz_parameterize(d::AbstractAccretionDisc, ρ) = SVector(ρ, cross_section(d, ρ))
 
 function distance_to_disc(d::AbstractThickAccretionDisc, x4; gtol)
-    height = cross_section(d, x4)
+    height = cross_section(d, _equatorial_project(x4))
     if height <= 0
         return one(eltype(x4))
     end
-    z = @inbounds x4[2] * cos(x4[3])
+    z = _spinaxis_project(x4)
     abs(z) - height - (gtol * x4[2])
 end
 
-function _cartesian_tangent_vector(ρ, d::AbstractThickAccretionDisc)
-    function _parameterization(r)
-        xz_parameterize(d, r)
+function _cartesian_tangent_vector(d::AbstractThickAccretionDisc, ρ)
+    function _parameterization(ρ̄)
+        xz_parameterize(d, ρ̄)
     end
     ∇f = ForwardDiff.derivative(_parameterization, ρ)
     v = SVector(∇f[1], 0, ∇f[2])
     v ./ √(v[1]^2 + v[2]^2 + v[3]^2)
 end
 
-function _cartesian_surface_normal(ρ, d::AbstractThickAccretionDisc)
-    ∇f = _cartesian_tangent_vector(ρ, d)
+function _cartesian_surface_normal(d::AbstractThickAccretionDisc, ρ)
+    ∇f = _cartesian_tangent_vector(d, ρ)
     # rotate 90 degrees in about ϕ̂
     SVector(-∇f[3], ∇f[2], ∇f[1])
 end
 
-function _rotate_cartesian_about_z(n, ϕ)
-    SVector(n[1] * cos(ϕ), n[1] * sin(ϕ), n[3])
-end
-
-function _cartesian_surface_normal(ρ, ϕ, d::AbstractThickAccretionDisc)
-    _rotate_cartesian_about_z(_cartesian_surface_normal(ρ, d), ϕ)
-end
+_cartesian_surface_normal(d::AbstractThickAccretionDisc, ρ, ϕ) =
+    _rotate_about_spinaxis(_cartesian_surface_normal(d, ρ), ϕ)
 
 export ThickDisc

src/tracing/precision-solvers.jl

@@ -201,7 +201,7 @@ function impact_parameters_for_radius_obscured(
         save_on = false,
         trajectories = length(α),
     )
-    n = _cartesian_surface_normal(radius, d)
+    n = _cartesian_surface_normal(d, radius)
     I = [!_is_visible(m, d, gp, n) for gp in gps]
     α[I] .= NaN
     β[I] .= NaN
@@ -219,7 +219,7 @@ function _is_visible(m::AbstractMetric, d, gp::AbstractGeodesicPoint, n::SVector
     # geodesics sufficiently close, test the angle
     v = SVector(gp_new.v[2], gp_new.v[3], gp_new.v[4]) ./ gp_new.v[1]
     v_geodesic = _spher_to_cart_jacobian(gp_new.x[3], gp_new.x[4], gp_new.x[2]) * v
-    n_rot = _rotate_cartesian_about_z(n, gp_new.x[4])
+    n_rot = _rotate_about_spinaxis(n, gp_new.x[4])
 
     X1 = to_cartesian(gp.x)
     X2 = to_cartesian(gp_new.x)

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

@@ -218,7 +218,7 @@ function _rear_workhorse(
         jacobian_disc = d,
         kwargs...,
     )
-    n = _cartesian_surface_normal(rₑ, d)
+    n = _cartesian_surface_normal(d, rₑ)
     function _thick_workhorse(θ::T)::NTuple{4,T} where {T}
         g, J, gp, _, _ = datum_workhorse(θ)
         is_visible = _is_visible(m, d, gp, n)

src/utils.jl

@@ -127,6 +127,12 @@ _equatorial_project(r, θ) = r * sin(abs(θ))
 _equatorial_project(x::SVector{4}) = _equatorial_project(x[2], x[3])
 _equatorial_project(x::SVector{8}) = _equatorial_project(x[2], x[3])
 
+_spinaxis_project(r, θ) = r * cos(abs(θ))
+_spinaxis_project(x::SVector{4}) = _spinaxis_project(x[2], x[3])
+_spinaxis_project(x::SVector{8}) = _spinaxis_project(x[2], x[3])
+
+_rotate_about_spinaxis(n::SVector{3}, ϕ) = SVector(n[1] * cos(ϕ), n[1] * sin(ϕ), n[3])
+
 _zero_if_nan(x::T) where {T} = isnan(x) ? zero(T) : x
 
 export cartesian_squared_distance, cartesian_distance, spherical_to_cartesian

test/discs/test-geometry.jl

@@ -0,0 +1,12 @@
+using Gradus, Test
+
+m = KerrMetric(1.0, 0.9)
+d = ShakuraSunyaev(m)
+n = Gradus._cartesian_tangent_vector(d, 2.6)
+@test n ≈ [0.689693957000099, 0.0, 0.724100991352412] atol = 1e-5
+
+n = Gradus._cartesian_tangent_vector(d, 6.6)
+@test n ≈ [0.9679192396299138, 0.0, 0.2512615083021063] atol = 1e-5
+
+n = Gradus._cartesian_tangent_vector(d, 1.0)
+@test n ≈ [1, 0, 0] atol = 1e-5

test/discs/test-polish-doughnut.jl

@@ -6,7 +6,7 @@ d = PolishDoughnut(m, rₖ = 12.0, n = 0.21);
 d.inner_radius;
 
 r = collect(range(10.0, 15.0, 200))
-h = map(x -> Gradus.r_cross_section(d, x), r)
+h = map(x -> Gradus.cross_section(d, x), r)
 
 # fingerprint the cross section map
 @test sum(h) ≈ 219.97440610254944 atol = 1e-5

test/runtests.jl

@@ -57,6 +57,7 @@ end
 
 @time @testset "geometry" verbose = true begin
     include("discs/test-polish-doughnut.jl")
+    include("discs/test-geometry.jl")
 end
 
 # little bit of aqua

test/smoke-tests/rendergeodesics.jl

@@ -4,12 +4,11 @@ using StaticArrays
 
 # Tests to make sure the basic `rendergeodesics` function works for (ideally) all metrics.
 
-function _thick_disc(u)
-    r = u[2]
-    if r < 9.0 || r > 11.0
+function _thick_disc(ρ)
+    if ρ < 9.0 || ρ > 11.0
         return -1.0
     else
-        x = r - 10.0
+        x = ρ - 10.0
         sqrt(1 - x^2)
     end
 end