Merge pull request #187 from astro-group-bristol/fergus/noZ-metric

feat: implemented NoZ metric

src/Gradus.jl

@@ -435,6 +435,7 @@ 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/noz-metric.jl")
 include("metrics/dilaton-axion-ad.jl")
 include("metrics/bumblebee-ad.jl")
 include("metrics/kerr-newman-ad.jl")

src/geometry/discs/thin-disc.jl

@@ -25,4 +25,43 @@ function distance_to_disc(d::ThinDisc, x4; gtol)
     _spinaxis_project(x4, signed = false) - _gtol_error(gtol, x4)
 end
 
-export ThinDisc
+"""
+    struct WarpedThinDisc{T} <: AbstractAccretionDisc{T}
+    WarpedThinDisc(inner_radius::T, outer_radius::T, f::F)
+
+Similar to [`ThinDisc`](@ref), however the scale height as a function of the radial coordinate may be specified by an arbitrary function. The function should return the ``h = \\cos(\\theta)`` signed scaled height (c.f. [`ThickDisc`](@ref) which is unsigned).
+
+The function should be of the form
+
+```julia
+function scale_height(ρ)
+    # ...
+end
+```
+"""
+struct WarpedThinDisc{T,F} <: AbstractAccretionDisc{T}
+    f::F
+    inner_radius::T
+    outer_radius::T
+end
+
+WarpedThinDisc(f; inner_radius = 0.0, outer_radius = 500.0) = WarpedThinDisc(f, inner_radius, outer_radius)
+
+optical_property(::Type{<:WarpedThinDisc}) = OpticallyThin()
+
+inner_radius(disc::WarpedThinDisc) = disc.inner_radius
+
+function distance_to_disc(d::WarpedThinDisc, x4; gtol)
+    ρ = _equatorial_project(x4)
+    if ρ < d.inner_radius || ρ > d.outer_radius
+        return one(eltype(x4))
+    end
+
+    h = d.f(ρ)
+    γ = _spinaxis_project(x4, signed = true)
+
+    abs(h - γ) - _gtol_error(gtol, x4)
+end
+
+export ThinDisc, WarpedThinDisc
+

src/line-profiles.jl

@@ -88,6 +88,7 @@ function lineprofile(
     # this 5maxrₑ is a heuristic that is insulting
     # todo: make it friendly
     plane = PolarPlane(GeometricGrid(); Nr = 450, Nθ = 1300, r_max = 5maxrₑ),
+    callback = domain_upper_hemisphere(),
     solver_args...,
 ) where {T}
     progress_bar = init_progress_bar("Lineprofile: ", trajectory_count(plane), verbose)
@@ -101,7 +102,7 @@ function lineprofile(
         save_on = false,
         verbose = verbose,
         progress_bar = progress_bar,
-        callback = domain_upper_hemisphere(),
+        callback = callback,
         ensemble = EnsembleEndpointThreads(),
         solver_args...,
     )

src/metrics/noz-metric.jl

@@ -0,0 +1,100 @@
+
+module __NoZMetric
+using ..StaticArrays
+using ..MuladdMacro
+
+@muladd @fastmath begin
+
+    epsilon(M, a, ϵ, y) = ϵ * M * a * y
+
+    # the way this function must be defined is a little complex
+    # but helps with type-stability
+    function metric_components(M, a, ϵ, rθ)
+        (r, θ) = rθ
+        sinθ2 = sin(θ)^2
+
+        y = cos(θ)
+        _epsilon = epsilon(M, a, ϵ, y)
+
+        tt =
+            -1 + ((2M * r * (r^(2) + a^(2) * y^(2))) /
+             ((r^(2) + a^(2) * y^(2))^(2) + (r^(2) - 2M * r + a^(2) * y^(2)) * _epsilon))
+        ϕϕ =
+            (
+                (1 - y^2) * (r^2 + a^2 * y^2 + _epsilon) * (
+                    r^4 +
+                    a^4 * y^2 +
+                    r^2 * (a^2 + a^2 * y^2 + _epsilon) +
+                    a^2 * _epsilon +
+                    2M * r * (a^2 - a^2 * y^2 - _epsilon)
+                )
+            ) / ((r^2 + a^2 * y^2)^2 + (r^2 - 2M * r + a^2 * y^2) * _epsilon)
+        rr = (r^(2) + a^(2) * y^(2) + _epsilon) / (r^(2) - 2M * r + a^(2))
+        yy = (r^(2) + a^(2) * y^(2) + _epsilon) / (1 - y^(2))
+
+        tϕ =
+            -(2M * r * a*(1 - y^(2)) * (r^(2) + a^(2) * y^(2) + _epsilon)) /
+            ((r^(2) + a^(2) * y^(2))^(2) + (r^(2) - 2M * r + a^(2) * y^(2)) * _epsilon)
+
+        # include the coordinate transformation factor
+        #   y = cos(θ)  -->  dy^2 = dθ^2 * sin(θ)^2
+        @SVector [tt, rr, yy * sinθ2, ϕϕ, tϕ]
+    end
+end
+
+end # module
+
+# new structure for our spacetime
+"""
+    struct NoZMetric
+"""
+@with_kw struct NoZMetric{T} <: AbstractStaticAxisSymmetric{T}
+    @deftype T
+    "Black hole mass."
+    M = 1.0
+    "Black hole spin."
+    a = 0.0
+    "Deviation parameter"
+    ϵ = 0.0
+end
+
+# implementation
+metric_components(m::NoZMetric, rθ) = __NoZMetric.metric_components(m.M, m.a, m.ϵ, rθ)
+inner_radius(m::NoZMetric) = m.M + √(m.M^2 - m.a^2)
+
+function _solve_orbit_θ(m, r)
+    function _objective(θ)
+        rθ = SVector(r, θ)
+        _, J = Gradus.metric_jacobian(m, rθ)
+        ∂rg = J[:, 1]
+        ∂θg = J[:, 2]
+        Ωϕ = Gradus.CircularOrbits._Ω_analytic(∂rg, false)
+
+        ∂θg[1] + 2 * ∂θg[5] * Ωϕ + ∂θg[4] * Ωϕ^2
+    end
+
+    Gradus.Roots.find_zero(_objective, π / 2)
+end
+
+function isco(m::NoZMetric{T}) where {T}
+    kerr_isco = isco(KerrMetric(M = m.M, a = m.a))
+    rs = range(inner_radius(m), kerr_isco + 1.0, 100)
+    thetas = map(rs) do r
+       try
+          _solve_orbit_θ(m, r)
+       catch
+           zero(T)
+       end
+    end
+
+    # TODO: this is such a hack
+    replace!(thetas, zero(T) => thetas[findfirst(!=(zero(T)), thetas)])
+
+    interp = _make_interpolation(rs, thetas)
+
+    dE(r) = ForwardDiff.derivative(x -> CircularOrbits.energy(m, SVector(x, interp(x))), r)
+    # optimize from the Kerr equivalent metric
+    Roots.find_zero(dE, kerr_isco)
+end
+
+export NoZMetric