Feat!: disc refactoring

benchmark/integrator/benchmark-tracing.jl

@@ -6,7 +6,7 @@ tracing_suite = BenchmarkGroup(["tracing"])
 m = KerrMetric(M = 1.0, a = 0.0)
 u = SVector(0.0, 1000.0, deg2rad(75.0), 0.0)
 v = map_impact_parameters(m, u, 4.0, 0.0)
-d = GeometricThinDisc(0.0, 1000.0, π / 2)
+d = ThinDisc(0.0, 1000.0)
 time_span = (0.0, 2000.0)
 
 # single position, single velocity

src/Gradus.jl

@@ -424,12 +424,22 @@ include("rendering/rendering.jl")
 include("rendering/utility.jl")
 
 include("tracing/method-implementations/first-order.jl")
-
 include("point-functions.jl")
 
 include("orbits/circular-orbits.jl")
 include("orbits/orbit-solving.jl")
 
+include("metrics/kerr-metric.jl")
+include("metrics/kerr-metric-first-order.jl")
+include("metrics/johannsen-ad.jl")
+include("metrics/johannsen-psaltis-ad.jl")
+include("metrics/morris-thorne-ad.jl")
+include("metrics/kerr-refractive-ad.jl")
+include("metrics/dilaton-axion-ad.jl")
+include("metrics/bumblebee-ad.jl")
+include("metrics/kerr-newman-ad.jl")
+include("metrics/minkowski.jl")
+
 include("geometry/geometry.jl")
 include("geometry/intersections.jl")
 include("geometry/discs.jl")
@@ -451,24 +461,12 @@ include("corona/flux-calculations.jl")
 include("corona/emissivity.jl")
 include("corona/spectra.jl")
 
-include("metrics/kerr-metric.jl")
-include("metrics/kerr-metric-first-order.jl")
-include("metrics/johannsen-ad.jl")
-include("metrics/johannsen-psaltis-ad.jl")
-include("metrics/morris-thorne-ad.jl")
-include("metrics/kerr-refractive-ad.jl")
-include("metrics/dilaton-axion-ad.jl")
-include("metrics/bumblebee-ad.jl")
-include("metrics/kerr-newman-ad.jl")
-include("metrics/minkowski.jl")
-
 include("special-radii.jl")
 include("redshift.jl")
 include("const-point-functions.jl")
 
 include("line-profiles.jl")
 include("reverberation.jl")
-include("geometry/polish-doughnut.jl")
 
 include("plotting-recipes.jl")
 

src/corona/emissivity.jl

@@ -128,7 +128,7 @@ All other keyword arguments are forwarded to [`tracegeodesics`](@ref).
 
 ```julia
 m = KerrMetric()
-d = GeometricThinDisc(Gradus.isco(m), 1000.0, π/2)
+d = ThinDisc(Gradus.isco(m), 1000.0)
 model = LampPostModel(h = 10.0)
 
 profile = emissivity_profile(m, d, model; n_samples = 128)

src/corona/profiles/voronoi.jl

@@ -10,11 +10,6 @@ struct VoronoiDiscProfile{D,V,G} <: AbstractDiscProfile
         gen::Vector{V},
         gps::Vector{G},
     ) where {D<:AbstractAccretionDisc,V<:AbstractArray,G}
-        if !isapprox(d.inclination, π / 2)
-            return error(
-                "Currently only supported for discs in the equatorial plane (θ=π/2).",
-            )
-        end
         new{D,V,G}(d, polys, gen, gps)
     end
 end

src/geometry/discs.jl

@@ -47,173 +47,6 @@ cross_section(d::AbstractThickAccretionDisc, x4) =
 r_cross_section(d::AbstractThickAccretionDisc, r::Number) =
     cross_section(d, SVector(0.0, r, π / 2, 0.0))
 
-"""
-    struct GeometricThinDisc{T} <: AbstractAccretionDisc{T}
-    GeometricThinDisc(inner_radius::T, outer_radius::T, inclination::T)
-
-Simple geometrically thin accretion disc spanning from `inner_radius` to `outer_radius` in
-gravitational units. Inclination of the disc is relative to spin axis, with ``90^\\circ``
-being perpendicular to the spin axis.
-"""
-@with_kw struct GeometricThinDisc{T} <: AbstractAccretionDisc{T}
-    inner_radius::T
-    outer_radius::T
-    inclination::T
-end
-optical_property(::Type{<:GeometricThinDisc}) = OpticallyThin()
-
-@fastmath function distance_to_disc(d::GeometricThinDisc{T}, x4; gtol) where {T}
-    p = @inbounds let r = x4[2], θ = x4[3], ϕ = x4[4]
-        if r < d.inner_radius || r > d.outer_radius
-            return 1.0
-        end
-        sinθ = sin(θ)
-        @SVector [r * sinθ * cos(ϕ), r * sinθ * sin(ϕ), r * cos(θ)]
-    end
-    n = @SVector [T(0), cos(d.inclination), sin(d.inclination)]
-    # project u into normal vector n
-    k = p ⋅ n
-    abs(k) - (gtol * x4[2])
-end
-
-struct DatumPlane{T} <: AbstractAccretionDisc{T}
-    height::T
-end
-optical_property(::Type{<:DatumPlane}) = OpticallyThin()
-
-function distance_to_disc(d::DatumPlane{T}, x4; gtol) where {T}
-    h = @inbounds let r = x4[2], θ = x4[3], ϕ = x4[4]
-        r * cos(θ)
-    end
-    abs(h) - d.height - (gtol * x4[2])
-end
-
-function datumplane(disc::AbstractThickAccretionDisc, r::T) where {T}
-    h = cross_section(disc, SVector{4,T}(0, r, T(π / 2), 0))
-    DatumPlane(h)
-end
-function datumplane(disc::AbstractThickAccretionDisc, x::SVector{4})
-    h = cross_section(disc, x)
-    DatumPlane(h)
-end
-
-"""
-    ThickDisc{T,F,P} <: AbstractThickAccretionDisc{T}
-    ThickDisc(f, params=nothing; T = Float64)
-
-A standard wrapper for creating custom disc profiles from height cross-section function `f`. This function
-is given the disc parameters as unpacked arguments:
-
-```julia
-d.f(u, d.params...)
-```
-
-If no parameters are specified, none will be passed to `f`.
-
-## Example
-
-Specifying a toroidal disc centered on `r=10` with radius 1:
-
-```julia
-d = ThickDisc() do u
-    r = u[2]
-    if r < 9.0 || r > 11.0
-        return -1.0
-    else
-        x = r - 10.0
-        sqrt(1-x^2)
-    end
-end
-```
-
-"""
-struct ThickDisc{T,F,P} <: AbstractThickAccretionDisc{T}
-    f::F
-    params::P
-end
-
-function ThickDisc(cross_section::F) where {F}
-    # todo: float bit generic???
-    ThickDisc{Float64,F,Nothing}(cross_section, nothing)
-end
-
-function ThickDisc(cross_section::F, parameters::P) where {F,P}
-    # todo: float bit generic???
-    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)
-
-function distance_to_disc(d::AbstractThickAccretionDisc, x4; gtol)
-    height = cross_section(d, x4)
-    if height <= 0.0
-        return 1.0
-    end
-    z = @inbounds x4[2] * cos(x4[3])
-    abs(z) - height - (gtol * x4[2])
-end
-
-# common thick disc models
-
-"""
-    ShakuraSunyaev{T} <: AbstractThickAccretionDisc{T}
-    ShakuraSunyaev(
-        m::AbstractMetric;
-        eddington_ratio = 0.3,
-        η = nothing,
-        contra_rotating = false,
-    )
-
-The classic Shakura & Sunyaev (1973) accretion disc model, with height given by ``2H``, where
-
-```math
-H = \\frac{3}{2} \\frac{1}{\\eta} \\left( \\frac{\\dot{M}}{\\dot{M}_\\text{Edd}} \\right) \\left( 1 - \\sqrt{\\frac{r_\\text{isco}}{\\rho}} \\right)
-```
-
-Here ``\\eta`` is the radiative efficiency, which, if unspecified, is determined by the circular orbit energy at the ISCO:
-
-```math
-\\eta = 1 - E_\\text{isco}
-```
-"""
-struct ShakuraSunyaev{T} <: AbstractThickAccretionDisc{T}
-    Ṁ::T
-    Ṁedd::T
-    η::T
-    risco::T
-end
-
-@fastmath function cross_section(d::ShakuraSunyaev, u)
-    @inbounds let r = u[2], θ = u[3]
-        if r < d.risco
-            return -1.0
-        end
-        ρ = r * sin(θ)
-        H = (3 / 2) * inv(d.η) * (d.Ṁ / d.Ṁedd) * (1 - sqrt(d.risco / ρ))
-        2H
-    end
-end
-
-function ShakuraSunyaev(
-    m::AbstractMetric{T};
-    eddington_ratio = 0.3,
-    η = nothing,
-    contra_rotating = false,
-) where {T}
-    r_isco = isco(m)
-    radiative_efficiency = if isnothing(η)
-        1 - CircularOrbits.energy(
-            m,
-            SVector{2}(r_isco, π / 2);
-            contra_rotating = contra_rotating,
-        )
-    else
-        η
-    end
-    ShakuraSunyaev(T(eddington_ratio), 1.0, radiative_efficiency, r_isco)
-end
-
 struct EllipticalDisc{T} <: AbstractAccretionDisc{T}
     inner_radius::T
     semi_major::T
@@ -253,10 +86,10 @@ function distance_to_disc(d::PrecessingDisc, x4; gtol)
     distance_to_disc(d.disc, x4prime; gtol)
 end
 
-export AbstractThickAccretionDisc,
-    ThickDisc,
-    ShakuraSunyaev,
-    GeometricThinDisc,
-    cross_section,
-    PrecessingDisc,
-    EllipticalDisc
+include("discs/datum-plane.jl")
+include("discs/polish-doughnut.jl")
+include("discs/shakura-sunyaev.jl")
+include("discs/thick-disc.jl")
+include("discs/thin-disc.jl")
+
+export AbstractThickAccretionDisc, cross_section, PrecessingDisc, EllipticalDisc

src/geometry/discs/datum-plane.jl

@@ -0,0 +1,22 @@
+struct DatumPlane{T} <: AbstractAccretionDisc{T}
+    height::T
+end
+optical_property(::Type{<:DatumPlane}) = OpticallyThin()
+
+function distance_to_disc(d::DatumPlane{T}, x4; gtol) where {T}
+    h = @inbounds let r = x4[2], θ = x4[3], ϕ = x4[4]
+        r * cos(θ)
+    end
+    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)
+    DatumPlane(h)
+end
+
+export DatumPlane

src/geometry/discs/shakura-sunyaev.jl

@@ -0,0 +1,59 @@
+"""
+    ShakuraSunyaev{T} <: AbstractThickAccretionDisc{T}
+    ShakuraSunyaev(
+        m::AbstractMetric;
+        eddington_ratio = 0.3,
+        η = nothing,
+        contra_rotating = false,
+    )
+
+The classic Shakura & Sunyaev (1973) accretion disc model, with height given by ``2H``, where
+
+```math
+H = \\frac{3}{2} \\frac{1}{\\eta} \\left( \\frac{\\dot{M}}{\\dot{M}_\\text{Edd}} \\right) \\left( 1 - \\sqrt{\\frac{r_\\text{isco}}{\\rho}} \\right)
+```
+
+Here ``\\eta`` is the radiative efficiency, which, if unspecified, is determined by the circular orbit energy at the ISCO:
+
+```math
+\\eta = 1 - E_\\text{isco}
+```
+"""
+struct ShakuraSunyaev{T} <: AbstractThickAccretionDisc{T}
+    Ṁ::T
+    Ṁedd::T
+    η::T
+    risco::T
+end
+
+@fastmath function cross_section(d::ShakuraSunyaev, u)
+    @inbounds let r = u[2], θ = u[3]
+        if r < d.risco
+            return -1.0
+        end
+        ρ = r * sin(θ)
+        H = (3 / 2) * inv(d.η) * (d.Ṁ / d.Ṁedd) * (1 - sqrt(d.risco / ρ))
+        2H
+    end
+end
+
+function ShakuraSunyaev(
+    m::AbstractMetric{T};
+    eddington_ratio = 0.3,
+    η = nothing,
+    contra_rotating = false,
+) where {T}
+    r_isco = isco(m)
+    radiative_efficiency = if isnothing(η)
+        1 - CircularOrbits.energy(
+            m,
+            SVector{2}(r_isco, π / 2);
+            contra_rotating = contra_rotating,
+        )
+    else
+        η
+    end
+    ShakuraSunyaev(T(eddington_ratio), 1.0, radiative_efficiency, r_isco)
+end
+
+export ShakuraSunyaev