feat!: rename GeometricThinDisc to ThinDisc + disc refactor

Put each disc in its own file and directory, so it's easy to see what belongs to what. Removed the inclination argument from ThinDisc, since with PrecessingDisc and the like, it would be better to add general modifiers.
astro-group-bristol · Sep 7, 2023 · 9ca4d9c · 9ca4d9c
1 parent f559756
commit 9ca4d9c
Show file tree

Hide file tree

Showing 25 changed files with 224 additions and 227 deletions.
diff --git a/src/Gradus.jl b/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")
 

diff --git a/src/corona/emissivity.jl b/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)

diff --git a/src/corona/profiles/voronoi.jl b/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

diff --git a/src/geometry/discs.jl b/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}
-    Ṁ::T
-    Ṁ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.Ṁ / d.Ṁ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
diff --git a/src/geometry/discs/datum-plane.jl b/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
diff --git a/src/geometry/polish-doughnut.jl → src/geometry/discs/polish-doughnut.jl b/src/geometry/polish-doughnut.jl → src/geometry/discs/polish-doughnut.jl
diff --git a/src/geometry/discs/shakura-sunyaev.jl b/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}
+    Ṁ::T
+    Ṁ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.Ṁ / d.Ṁ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