Feat: Bumblebee metric (#93)

* feat: bumblebee metric in slow rotation regime

* chore: bump version

* style: formatting

* test: bumblebee metric

* fix: reach under the hood for missing function

* fix: remove show

Project.toml

@@ -1,7 +1,7 @@
 name = "Gradus"
 uuid = "c5b7b928-ea15-451c-ad6f-a331a0f3e4b7"
 authors = ["fjebaker <[email protected]>"]
-version = "0.2.10"
+version = "0.2.11"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

src/Gradus.jl

@@ -190,6 +190,7 @@ 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("special-radii.jl")
 include("redshift.jl")

src/metrics/bumblebee-ad.jl

@@ -0,0 +1,52 @@
+module __BumblebeeAD
+using ..StaticArrays
+using ..MuladdMacro
+
+@muladd @fastmath begin
+    Δ(r, M, l) = (r^2 - (2M * r)) / (l + 1)
+
+    function metric_components(M, a, l, rθ)
+        (r, θ) = rθ
+        sinθ2 = sin(θ)^2
+        Δ₀ = Δ(r, M, l)
+
+        tt = -(1 - (2M / r))
+        rr = r^2 / Δ₀
+        θθ = r^2
+        ϕϕ = r^2 * sinθ2
+
+        tϕ = -2M * a * sinθ2 / r
+        @SVector [tt, rr, θθ, ϕϕ, tϕ]
+    end
+end
+
+end # module
+
+struct BumblebeeMetric{T} <: AbstractAutoDiffStaticAxisSymmetricParams{T}
+    "Black hole mass."
+    M::T
+    "Black hole spin."
+    a::T
+    "LSB parameter."
+    l::T
+    function BumblebeeMetric(M::T, a, l) where {T}
+        if l <= -1.0
+            error("l must be >-1")
+        end
+        if abs(a) > 0.3
+            error(
+                "This metric is for the slow rotation approximation only, and requires |a| < 0.3.",
+            )
+        end
+        new{T}(T(M), T(a), T(l))
+    end
+end
+BumblebeeMetric(; M = 1.0, a = 0.0, l = 0.0) = BumblebeeMetric(M, a, l)
+
+# implementation
+metric_components(m::BumblebeeMetric{T}, rθ) where {T} =
+    __BumblebeeAD.metric_components(m.M, m.a, m.l, rθ)
+inner_radius(m::BumblebeeMetric{T}) where {T} = m.M + √(m.M^2 - m.a^2)
+
+
+export BumblebeeMetric

src/tracing/method-implementations/auto-diff.jl

@@ -261,8 +261,8 @@ but non-allocating.
 function metric_jacobian(m::AbstractAutoDiffStaticAxisSymmetricParams, rθ)
     f = x -> metric_components(m, x)
     T = typeof(ForwardDiff.Tag(f, eltype(rθ)))
-    ydual = ForwardDiff.static_dual_eval(T, f, rθ)
-    ForwardDiff.value.(T, ydual), ForwardDiff.extract_jacobian(T, ydual, rθ)
+    ydual = ForwardDiff.ForwardDiffStaticArraysExt.static_dual_eval(T, f, rθ)
+    ForwardDiff.value.(T, ydual), ForwardDiff.ForwardDiffStaticArraysExt.extract_jacobian(T, ydual, rθ)
 end
 
 @inbounds function geodesic_eq(

test/smoke-tests/rendergeodesics.jl

@@ -18,12 +18,13 @@ end
     u = @SVector [0.0, 100.0, deg2rad(85), 0.0]
 
     # last computed 21/01/2023: shrink resolution
-    expected = (8969.1564582409967, 8969.15634220181, 8977.502920124776, 413.4963434133726)
+    expected = (8969.1564582409967, 8969.15634220181, 8977.502920124776, 413.4963434133726, 8969.155411207657)
     result = map((
         KerrMetric(),
         JohannsenMetric(),
         KerrSpacetimeFirstOrder(),
         MorrisThorneWormhole(),
+        BumblebeeMetric(),
     )) do m
         α, β, img = rendergeodesics(
             m,

test/smoke-tests/tracegeodesics.jl

@@ -30,6 +30,7 @@ using StaticArrays
             JohannsenMetric(),
             KerrSpacetimeFirstOrder(),
             MorrisThorneWormhole(),
+            BumblebeeMetric()
         ]
             test_single(m, u, v)
             test_many(m, u, v)