Automated circular orbit tracing

Project.toml

@@ -11,10 +11,11 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeometricalPredicates = "fd0ad045-b25c-564e-8f9c-8ef5c5f21267"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
+Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"

src/AccretionFormulae/AccretionFormulae.jl

@@ -6,10 +6,15 @@ import ..GradusBase: AbstractMetricParams, metric
 using ..FirstOrderMethods: FirstOrderGeodesicPoint
 using ..Rendering: PointFunction
 
+using Optim: optimize, minimizer, GoldenSection, Brent
+using DocStringExtensions
+using StaticArrays
+
 include("redshift.jl")
+include("orbit-discovery.jl")
 
 const redshift = PointFunction(_redshift_guard)
 
-export redshift
+export solve_equitorial_circular_orbit, trace_equitorial_circular_orbit, redshift
 
 end # module

src/AccretionFormulae/orbit-discovery.jl

@@ -0,0 +1,97 @@
+
+"""
+    $(TYPEDSIGNATURES)
+
+Quality of stability measure, which has a minima for circular orbits. Effectively
+a sum of the normalised residuals.
+"""
+Qs(rs) = sqrt(sum((rs ./ rs[1] .- 1.0) .^ 2) / length(rs))
+
+function trace_single_orbit(m, r, vϕ; max_time = 300.0, μ = 1.0, tracer_args...)
+    # fixed in equitorial plane
+    u = @SVector [0.0, r, deg2rad(90.0), 0.0]
+    v = @SVector [0.0, 0.0, 0.0, vϕ]
+    Gradus.tracegeodesics(m, u, v, (0.0, max_time); μ = μ, tracer_args...)
+end
+
+function measure_stability(m::AbstractMetricParams{T}, r, vϕ; tracer_args...) where {T}
+    sol = trace_single_orbit(m, r, vϕ; tracer_args...)
+    rs = selectdim(sol, 1, 6)
+    Qs(rs)
+end
+
+function __solve_equitorial_circular_orbit(
+    m::AbstractMetricParams{T},
+    r,
+    optimizer,
+    lower_bound,
+    upper_bound;
+    tracer_args...,
+) where {T}
+    res = optimize(
+        vϕ -> measure_stability(m, r, vϕ; tracer_args...),
+        lower_bound,
+        upper_bound,
+        optimizer,
+    )
+    minimizer(res)
+end
+
+function solve_equitorial_circular_orbit(
+    m::AbstractMetricParams{T},
+    r::Number;
+    lower_bound = 0.0,
+    upper_bound = 1.0,
+    optimizer = GoldenSection(),
+    tracer_args...,
+) where {T}
+    __solve_equitorial_circular_orbit(m, r, optimizer, lower_bound, upper_bound)
+end
+
+function sliding_window(func, N, lower_bound, upper_bound, lower_rate, upper_rate)
+    map(1:N) do i
+        res = func((i, lower_bound, upper_bound))
+        lower_bound = res * lower_rate
+        upper_bound = res * upper_rate
+        res
+    end
+end
+
+
+function solve_equitorial_circular_orbit(
+    m::AbstractMetricParams{T},
+    r_range::Union{AbstractRange,AbstractArray};
+    lower_bound = 0.0,
+    upper_bound = 1.0,
+    lower_rate = 0.98,
+    upper_rate = 1.5,
+    kwargs...,
+) where {T}
+    r_range_reverse = sort(r_range) |> reverse
+    candidate_vϕ = sliding_window(
+        length(r_range_reverse),
+        lower_bound,
+        upper_bound,
+        lower_rate,
+        upper_rate,
+    ) do (i, lower_bound, upper_bound)
+        r = r_range_reverse[i]
+        solve_equitorial_circular_orbit(
+            m,
+            r,
+            lower_bound = lower_bound,
+            upper_bound = upper_bound,
+        )
+    end
+    reverse!(candidate_vϕ)
+end
+
+function trace_equitorial_circular_orbit(
+    m::AbstractMetricParams{T},
+    rs;
+    kwargs...,
+) where {T}
+    map(zip(rs, solve_equitorial_circular_orbit(m, rs; kwargs...))) do (r, vϕ)
+        trace_single_orbit(m, r, vϕ; kwargs...)
+    end
+end

src/GeodesicTracer/GeodesicTracer.jl

@@ -61,8 +61,8 @@ function tracegeodesics(
         time_domain,
         solver;
         callback = nothing,
-        abstol = 1e-8,
-        reltol = 1e-8,
+        abstol = 1e-9,
+        reltol = 1e-9,
         solver_opts...,
     )
 end

src/Gradus.jl

@@ -73,6 +73,6 @@ include("special-radii.jl")
 include("AccretionFormulae/AccretionFormulae.jl")
 
 using .AccretionFormulae
-export redshift
+export solve_equitorial_circular_orbit, trace_equitorial_circular_orbit
 
 end # module

src/special-radii.jl

@@ -1,6 +1,27 @@
 """
-$(TYPEDSIGNATURES)
+    $(TYPEDSIGNATURES)
 
 Innermost stable circular orbit.
 """
 isco(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).")
+
+"""
+    $(TYPEDSIGNATURES)
+
+Photon orbit.
+"""
+r_ph(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).")
+
+"""
+    $(TYPEDSIGNATURES)
+
+Marginally bound orbit.
+"""
+r_mb(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).")
+
+"""
+    $(TYPEDSIGNATURES)
+
+Event horizon.
+"""
+r_s(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).")

test/smoke-tests/circular-orbits.jl

@@ -0,0 +1,30 @@
+# smoke test to make sure circular orbits work
+using Test, Gradus
+
+@testset "circular-orbits" begin
+
+    @testset "solve_equitorial_circular_orbit" begin
+        # only implemented for the BoyerLindquist metrics at the moment
+        # expected is sum of the circular orbit vϕ
+        # last updated: 22 Apr 2022
+        r_range = 6.0:0.5:10.0
+        for (m, expected) in [
+            (BoyerLindquistAD(M = 1.0, a = 0.0), 0.5432533297869712),
+            (BoyerLindquistAD(M = 1.0, a = 1.0), 0.5016710246454921),
+            (BoyerLindquistAD(M = 1.0, a = -1.0), 0.5993458160081419),
+            (JohannsenAD(M = 1.0, a = 1.0, α22 = 1.0), 0.4980454719932759),
+        ]
+            vϕs = solve_equitorial_circular_orbit(m, r_range)
+            @test isapprox(sum(vϕs), expected, atol = 1e-6, rtol = 0.0)
+        end
+    end
+
+    @testset "trace_equitorial_circular_orbit" begin
+        m = BoyerLindquistAD(M = 1.0, a = 0.0)
+        Gradus.isco(m)
+        r_range = 6.0:0.1:10.0
+        simsols = @time trace_equitorial_circular_orbit(m, r_range)
+        # smoke test passed
+        @test true
+    end
+end