Merge pull request #210 from astro-group-bristol/fergus/emda-fixes

Fixes to Dilaton Axion metric

.github/workflows/auto-register.yml

@@ -21,7 +21,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ['1.10']
+        julia-version: ['1.11']
         os: [ubuntu-latest]
     steps:
       - uses: actions/checkout@v2

.github/workflows/docs.yml

@@ -17,7 +17,7 @@ jobs:
       - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1.10'
+          version: '1.11'
       - run: |
           julia --project=docs -e '
             using Pkg

.github/workflows/test.yml

@@ -26,7 +26,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ['1.10']
+        julia-version: ['1.11']
         os: [ubuntu-latest]
 
     steps:

src/metrics/dilaton-axion-ad.jl

@@ -5,42 +5,41 @@ García et al. (1995).
 module __DilatonAxionAD
 using ..StaticArrays
 
-@fastmath Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2
-@fastmath Δ(r, a, rg) = r^2 - 2rg * r + a^2
-@fastmath A(r, a, Δ, θ) = (r^2 + a^2)^2 - a^2 * Δ * sin(θ)^2
+Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2
+Δ(r, R, a) = r^2 + a^2 - R * r
+# the a here is actually the dimensionless spin parameter a / R ?
+Σhat(Σ, β, b, r, βb, a, θ, R) = Σ - (β^2 + 2b * r) + R^2 * βb * (βb - 2 * a * cos(θ))
+Δhat(Δ, β, b, r, βb, R) = Δ - (β^2 + 2b * r) - R * (R + 2b) * βb^2
 
-@fastmath W(βab, θ, βa) = 1 + (βab * (2 * cos(θ) - βab) + βa^2) * csc(θ)^2
-@fastmath Σ̂(Σ, β, b, r, βb, a, θ, rg) =
-    Σ - (β^2 + 2b * r) + rg^2 * βb * (βb - 2 * (a / rg) * cos(θ))
-@fastmath Δ̂(Δ, β, b, r, βb, rg) = Δ - (β^2 + 2b * r) - rg * (rg + 2b) * βb^2
-@fastmath Â(δ, a, Δ̂, W, θ) = δ^2 - a^2 * Δ̂ * W^2 * sin(θ)^2
-@fastmath δ(r, b, a) = r^2 - 2b * r + a^2
+δ(r, b, a) = r^2 - 2b * r + a^2
+
+W(θ, βab, βa) = 1 + (βab * (2 * cos(θ) - βab) + βa^2) * csc(θ)^2
+A(a, θ, δ, W, Δhat) = δ^2 - Δhat * (W * a * sin(θ))^2
 
 @fastmath function metric_components(M, a, β, b, rθ)
     (r, θ) = rθ
 
-    rg = M
-
-    Δ₀ = Δ(r, a, rg)
-    Σ₀ = Σ(r, a, θ)
+    R = M
 
-    βa = β / a
+    # check these for the normalisation term R
     βb = β / b
-    βab = rg * β / (a * b)
+    βa = β / a
+    βab = β / (a * b)
 
-    Δ̂₀ = Δ̂(Δ₀, β, b, r, βb, rg)
-    Σ̂₀ = Σ̂(Σ₀, β, b, r, βb, a, θ, rg)
+    Σ₀ = Σ(r, a, θ)
+    Δ₀ = Δ(r, R, a)
+    Δhat₀ = Δhat(Δ₀, β, b, r, βb, R)
+    Σhat₀ = Σhat(Σ₀, β, b, r, βb, a, θ, R)
     δ₀ = δ(r, b, a)
-    W₀ = W(βab, θ, βa)
-
-    Â₀ = Â(δ₀, a, Δ̂₀, W₀, θ)
+    W₀ = W(θ, βab, βa)
+    A₀ = A(a, θ, δ₀, W₀, Δhat₀)
 
-    tt = -(Δ̂₀ - a^2 * sin(θ)^2) / Σ̂₀
-    rr = Σ̂₀ / Δ̂₀
-    θθ = Σ̂₀
-    ϕϕ = Â₀ * sin(θ)^2 / Σ̂₀
+    tt = -(Δhat₀ - a^2 * sin(θ)^2) / Σhat₀
+    rr = Σhat₀ / Δhat₀
+    θθ = Σhat₀
+    ϕϕ = A₀ * sin(θ)^2 / Σhat₀
 
-    tϕ = -a * (δ₀ - Δ̂₀ * W₀) * sin(θ)^2 / Σ̂₀
+    tϕ = -a * (δ₀ - Δhat₀ * W₀) * sin(θ)^2 / Σhat₀
     @SVector [tt, rr, θθ, ϕϕ, tϕ]
 end
 
@@ -60,7 +59,7 @@ Einstein-Maxwell-Dilaton-Axion metric.
     "Singularity mass."
     M = 1.0
     "Singularity spin."
-    a = 0.0
+    a = 0.5
     "Dilaton coupling strength."
     β = 0.0
     "Axion coupling strength."
@@ -69,6 +68,7 @@ end
 
 metric_components(m::DilatonAxion{T}, rθ) where {T} =
     __DilatonAxionAD.metric_components(m.M, m.a, m.β, m.b, rθ)
-inner_radius(m::DilatonAxion{T}) where {T} = m.M + √(m.M^2 - m.a^2)
+inner_radius(m::DilatonAxion{T}) where {T} =
+    m.M + m.b + √((m.M + m.b)^2 - m.a^2 + m.β^2 - (m.M - 2m.b) * m.M * (m.β / m.b)^2)
 
 export DilatonAxion

src/orbits/orbit-solving.jl

@@ -100,6 +100,7 @@ struct PlungingInterpolation{M,_interp_type}
     m::M
     t::_interp_type
     r::_interp_type
+    θ::_interp_type
     ϕ::_interp_type
 
     function PlungingInterpolation(m::M, sol) where {M<:AbstractMetric{T}} where {T}
@@ -110,13 +111,15 @@ struct PlungingInterpolation{M,_interp_type}
 
         vt = sol[5, :][I]
         vr = sol[6, :][I]
+        vθ = sol[7, :][I]
         vϕ = sol[8, :][I]
 
         r_interp = _make_interpolation(r, vt)
         new{M,typeof(r_interp)}(
             m,
             r_interp,
             _make_interpolation(r, vr),
+            _make_interpolation(r, vθ),
             _make_interpolation(r, vϕ),
         )
     end
@@ -130,21 +133,22 @@ function (pintrp::PlungingInterpolation)(r)
     r_bounded = _enforce_interpolation_bounds(r, pintrp)
     vt = pintrp.t(r_bounded)
     vr = pintrp.r(r_bounded)
+    vθ = pintrp.r(r_bounded)
     vϕ = pintrp.ϕ(r_bounded)
-    SVector(vt, vr, 0, vϕ)
+    SVector(vt, vr, vθ, vϕ)
 end
 
 function interpolate_plunging_velocities(
     m::AbstractMetric{T};
     max_time = 50_000,
     contra_rotating = false,
     reltol = 1e-9,
+    δr = (reltol * 10),
     kwargs...,
 ) where {T}
     isco = Gradus.isco(m)
 
     # rule of thumb to achieve desired error
-    δr = (reltol * 10)
     u = @SVector([0.0, isco - δr, deg2rad(90), 0.0])
     v = Gradus.CircularOrbits.plunging_fourvelocity(
         m,

test/integration/test-precision.jl

@@ -6,24 +6,24 @@ m = KerrMetric(M = 1.0, a = 1.0)
 u = SVector(0.0, 1000.0, π / 2, 0.0)
 
 # target position
-target = SVector(10.0, 0.001, 0.0)
+target = SVector(10.0, 0.005, 0.0)
 
 α, β, accuracy = Gradus.impact_parameters_for_target(target, m, u)
-@test α ≈ -0.0017523796800187875 rtol = 1e-3
+@test α ≈ -0.004013630261097743 rtol = 1e-3
 @test β ≈ 10.969606493445841 rtol = 1e-3
-@test accuracy ≈ 0.009091572709970781 rtol = 1e-3
+@test accuracy ≈ 0.004587323209289997 rtol = 1e-3
 
 # new target
 target = SVector(10.0, deg2rad(40), -π / 4)
 
 α, β, accuracy = Gradus.impact_parameters_for_target(target, m, u)
 @test α ≈ 4.848373364532467 rtol = 1e-3
 @test β ≈ 8.012499652493656 rtol = 1e-3
-@test accuracy ≈ 0.0049975644770909105 rtol = 1e-3
+@test accuracy ≈ 0.0030663177403400794 rtol = 1e-3
 
 # test geodesic point interface
 α, β, gp, accuracy = Gradus.optimize_for_target(target, m, u)
 @test gp.x[1] ≈ 1005.2700874611182 rtol = 1e-3
 @test α ≈ 4.848373364532467 rtol = 1e-3
 @test β ≈ 8.012499652493656 rtol = 1e-3
-@test accuracy ≈ 0.0049975644770909105 rtol = 1e-3
+@test accuracy ≈ 0.0030663177403400794 rtol = 1e-3

test/smoke-tests/special-radii.jl

@@ -12,7 +12,7 @@ using StaticArrays
         KerrMetric(1.0, -0.998),
         JohannsenMetric(M = 1.0, a = 0.998, α13 = 1.0),
         JohannsenMetric(M = 1.0, a = 0.998, α22 = 1.0),
-        DilatonAxion(M = 1.0, a = 0.998, β = 0.2, b = 1.0),
+        DilatonAxion(M = 1.0, a = 0.6, β = -0.5, b = 0.1),
     )
 
     @testset "iscos" begin
@@ -26,7 +26,7 @@ using StaticArrays
                 8.99437445480357,
                 2.8482863127671534,
                 1.1306596884484472,
-                6.880202293032178,
+                32.95283734688169,
             ],
         )
             isco_r = Gradus.isco(m)

test/transfer-functions/test-thick-disc.jl

@@ -8,15 +8,15 @@ d = ShakuraSunyaev(m)
 tf = cunningham_transfer_function(m, x, d, 3.0; β₀ = 1.0)
 
 total = sum(filter(!isnan, tf.f))
-@test total ≈ 12.253422180875667 atol = 1e-4
+@test total ≈ 12.255956053708246 atol = 1e-4
 
 m = KerrMetric(1.0, 0.2)
 x = SVector(0.0, 10_000, deg2rad(20), 0.0)
 d = ShakuraSunyaev(m; eddington_ratio = 0.2)
 
 tf = cunningham_transfer_function(m, x, d, 5.469668466100368; β₀ = 1.0)
 total = sum(filter(!isnan, tf.f))
-@test total ≈ 20.83469 atol = 1e-4
+@test total ≈ 20.83782174681574 atol = 1e-4
 
 # the transfer function here is pretty horrible as it's almost impossible to actually 
 # see; this is a test to make sure it doesn't error