New metric: Kerr with Dark Matter

src/metrics/kerr-dark-matter.jl

@@ -0,0 +1,74 @@
+module __KerrDarkMatter
+using ..StaticArrays
+using ..MuladdMacro
+using ..Gradus: _smooth_interpolate
+
+@muladd @fastmath begin
+    Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2
+    Δ(r, R, a) = r^2 + a^2 - R * r
+
+    function G(r, Δr, rₛ)
+        dr = (r - rₛ) / Δr
+        (3 - 2 * dr) * dr^2
+    end
+
+    function dark_matter_mass(r, Δr, rₛ, M)
+        if r < rₛ
+            zero(r)
+        elseif r < rₛ + Δr
+            M * G(r, Δr, rₛ)
+        else
+            M
+        end
+    end
+
+    # the way this function must be defined is a little complex
+    # but helps with type-stability
+    function metric_components(M_bh, a, M_dm, Δr, rₛ, rθ)
+        (r, θ) = rθ
+
+        M = M_bh + dark_matter_mass(r, Δr, rₛ, M_dm)
+
+        R = 2M
+        sinθ2 = sin(θ)^2
+        cosθ2 = (1 - sinθ2)
+        # slightly faster, especially when considering AD evals
+        Σ₀ = r^2 + a^2 * cosθ2
+
+        tt = -(1 - (R * r) / Σ₀)
+        rr = Σ₀ / Δ(r, R, a)
+        θθ = Σ₀
+        ϕϕ = sinθ2 * (r^2 + a^2 + (sinθ2 * R * r * a^2) / Σ₀)
+
+        tϕ = (-R * r * a * sinθ2) / Σ₀
+        @SVector [tt, rr, θθ, ϕϕ, tϕ]
+    end
+end
+
+end # module
+
+"""
+    struct KerrDarkMatter
+
+$(FIELDS)
+
+https://arxiv.org/pdf/2003.06829.pdf
+"""
+@with_kw struct KerrDarkMatter{T} <: AbstractStaticAxisSymmetric{T}
+    @deftype T
+    "Black hole mass."
+    M = 1.0
+    "Black hole spin."
+    a = 0.0
+    "Dark matter mass."
+    M_dark_matter = 2.0
+    Δr = 20.0
+    rₛ = 10.0
+end
+
+# implementation
+metric_components(m::KerrDarkMatter{T}, rθ) where {T} =
+    __KerrDarkMatter.metric_components(m.M, m.a, m.M_dark_matter, m.Δr, m.rₛ, rθ)
+inner_radius(m::KerrDarkMatter{T}) where {T} = m.M + √(m.M^2 - m.a^2)
+
+export KerrDarkMatter

src/metrics/kerr-refractive-ad.jl

@@ -1,4 +1,5 @@
 module __KerrRefractiveAD
+using ..Gradus: _smooth_interpolate
 using ..StaticArrays
 
 @fastmath Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2
@@ -22,18 +23,8 @@ using ..StaticArrays
     # need to interpolate the refractve index boundary
     # else we don't have a gradient to compute, and we miss it
     # hemorraging energy in the process
-    δr = 2.0
-    if r ≤ corona_radius
-        # ...
-    elseif corona_radius ≤ r ≤ corona_radius + δr
-        t = (r - corona_radius) / δr
-        # use an arbitrarily steep smooth interpolation
-        # this one isn't perfect, but does a good job
-        k = atan(1e5t) * 2 / π
-        n = k + n * (1 - k)
-    else
-        n = 1.0
-    end
+    t = _smooth_interpolate(r, corona_radius)
+    n = t + (1 - t) * n
 
     @SVector [tt / n^2, rr, θθ, ϕϕ, tϕ / n]
 end

src/utils.jl

@@ -150,4 +150,20 @@ _rotate_about_spinaxis(n::SVector{3}, ϕ) = SVector(n[1] * cos(ϕ), n[1] * sin(
 
 _zero_if_nan(x::T) where {T} = isnan(x) ? zero(T) : x
 
+@inline function _smooth_interpolate(
+    x::T,
+    x₀;
+    δx = T(2.5),
+    smoothing_offset = T(1e4),
+) where {T}
+    if x ≤ x₀
+        one(T)
+    elseif x₀ ≤ x ≤ x₀ + δx
+        t = (x - x₀) / δx
+        atan(smoothing_offset * t) * 2 / π
+    else
+        zero(T)
+    end
+end
+
 export cartesian_squared_distance, cartesian_distance, spherical_to_cartesian