Event horizons and Johannsen Psaltis (2011) metric (#41)

* added a function for plotting event horizons

* added Johannsen-Psaltis metric

* minimal docstring for JP metric

* updated docs with new metric and examples

* bump version

Project.toml

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

docs/Project.toml

@@ -1,3 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Gradus = "c5b7b928-ea15-451c-ad6f-a331a0f3e4b7"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

docs/make.jl

@@ -16,21 +16,22 @@ makedocs(
             "Point functions" => "overview/point-functions.md",
             # "Callbacks"
             "Available metrics" => "overview/metrics.md",
-            "Accretion geometry" => "overview/accretion-geometry.md"
+            "Accretion geometry" => "overview/accretion-geometry.md",
         ],
         # "Reverberation Lags"
         "Internals" => [
             "Geodesic integration" => "overview/geodesic-integration.md",
             "Implementing new metrics" => "internals/custom-metrics.md",
             "Special radii" => "internals/special-radii.md",
         ],
-        "Module API" => [
-            "Gradus" => "api-documentation/Gradus.md",
-            "GradusBase" => "api-documentation/GradusBase.md",
-            "GeodesicTracer" => "api-documentation/GeodesicTracer.md",
-            "AccretionGeometry" => "api-documentation/AccretionGeometry.md"
-        ] |> sort,
+        "Module API" =>
+            [
+                "Gradus" => "api-documentation/Gradus.md",
+                "GradusBase" => "api-documentation/GradusBase.md",
+                "GeodesicTracer" => "api-documentation/GeodesicTracer.md",
+                "AccretionGeometry" => "api-documentation/AccretionGeometry.md",
+            ] |> sort,
     ],
 )
 
-deploydocs(repo = "github.com/astro-group-bristol/Gradus.jl.git")
+# deploydocs(repo = "github.com/astro-group-bristol/Gradus.jl.git")

docs/src/internals/special-radii.md

@@ -9,4 +9,5 @@ isco
 r_ph
 r_mb
 r_s
+Gradus.event_horizon
 ```

docs/src/overview/examples.md

@@ -5,12 +5,20 @@ Pages = ["examples.md"]
 Depth = 3
 ```
 
+```@setup env
+using Gradus
+using StaticArrays
+using Plots
+ENV["GKSwstype"]="nul"
+Plots.default(show=false)
+```
+
 ## Redshift image
 
 !!! note
     The [`Gradus.ConstPointFunctions.redshift`](@ref) function is an analytic solution for redshift, which may not be implemented for every type of metric or disc geometry. See [Interpolating redshifts](@ref) for a more flexible numeric alternative.
 
-```julia
+```@example env
 using Gradus
 using StaticArrays
 using Plots
@@ -42,13 +50,11 @@ img = rendergeodesics(
 heatmap(img)
 ```
 
-![redshift](../assets/examples/redshift.png)
-
 ## Line profile
 
 Using the redshift example, we can bin a line-profile using [StatsBase.jl](https://juliastats.org/StatsBase.jl/stable/empirical/#StatsBase.Histogram). We'll calculate the iron line profile, with a delta-emission at 6.4 keV.
 
-```julia
+```@example env
 using StatsBase
 
 # remove nans and flatten the redshift image
@@ -63,13 +69,11 @@ lineprof = fit(Histogram, data, x_bins)
 plot(x_bins[1:end-1], lineprof.weights, seriestype = :steppre)
 ```
 
-![iron-line-profile](../assets/examples/line-profile.svg)
-
 ## Interpolating redshifts
 
 In cases where no analytic redshift solution is known, we can instead interpolate a numeric approximation. For example, interpolating the plunging region velocities and using the analytic solution for general static, axis symmetric metrics outside of the ISCO can be achieved with:
 
-```julia
+```@example env
 using Gradus
 using StaticArrays
 using Plots
@@ -103,6 +107,94 @@ img = rendergeodesics(
 heatmap(img)
 ```
 
-![redshift](../assets/examples/redshift-interpolated.png)
+For more complex disc geometry: TODO
+
+## ISCO
+
+The [Gradus.isco](@ref) may be calculated with a simple convenience function, as may the energy associated with a given stable circular orbit.
+
+```@example env
+using Gradus
+using Plots
+
+# prepare plot
+p = plot(legend=:bottomright, ylabel = "E", xlabel = "r", xscale = :log10)
+
+# choice of spin to plot energy curves for
+for a in [0.0, 0.4, 0.6]
+    m = BoyerLindquistAD(M = 1.0, a = a)
+
+    rs = range(Gradus.isco(m), 100.0, 500)
+    energy = map(rs) do r
+        Gradus.CircularOrbits.energy(m, r)
+    end
+
+    plot!(rs, energy, label = "a=$a")
+end
+
+# calculate the ISCO as a function of spin
+data = map(range(-1.0, 0.8, 100)) do a
+    m = BoyerLindquistAD(M = 1.0, a = a)
+    r = Gradus.isco(m)
+    Gradus.CircularOrbits.energy(m, r), r
+end
+
+# overlay onto plot
+plot!(last.(data), first.(data), color=:black, linestyle=:dash, label="ISCO")
+
+p
+```
+
+## Event horizons and naked singularities
+
+Here is an example of how to use [`Gradus.event_horizon`](@ref) to plot the shape of an event horizon in two dimensions. In the case of a naked singularity, as with the certain parameters combinations in the [`JohannsenPsaltisAD`](@ref) metric, we see a disconnected region in the plot.
+
+```@example env
+using Gradus
+using Plots
+
+function draw_horizon(p, m)
+    rs, θs = Gradus.event_horizon(m, resolution = 200)
+    radius = rs
+
+    x = @. radius * sin(θs)
+    y = @. radius * cos(θs)
+    plot!(p, x, y, label = "a = $(m.a)")
+end
+
+p = plot(aspect_ratio = 1)
+for a in [0.0, 0.5, 0.6, 0.7, 0.8]
+    m = JohannsenPsaltisAD(M = 1.0, a = a, ϵ3 = 2.0)
+    draw_horizon(p, m)
+end
+p
+```
+
+We can also calculate parameter combinations that lead to naked singularities, and plot the parameter space domains to show exclusion zones:
+
+```@example env
+function calc_exclusion(as, ϵs)
+    regions = [
+        Gradus.is_naked_singularity(JohannsenPsaltisAD(M = 1.0, a = a, ϵ3 = ϵ))
+        for a in as, ϵ in ϵs
+    ]
+
+    map(i -> i ? 1.0 : NaN, regions)
+end
+
+# define ranges (small in this example as a little computationally intense)
+as = range(0, 1.0, 40)
+ϵs = range(-10, 10, 40)
+
+img = calc_exclusion(as, ϵs)
+heatmap(
+    as, 
+    ϵs, 
+    img', 
+    color = :black, 
+    colorbar = false, 
+    xlabel = "a", 
+    ylabel = "ϵ"
+)
+```
 
-For more complex disc geometry: TODO

docs/src/overview/metrics.md

@@ -25,6 +25,7 @@ BoyerLindquistFO
 ```@docs
 BoyerLindquistAD
 JohannsenAD
+JohannsenPsaltisAD
 KerrRefractiveAD
 DilatonAxionAD
 MorrisThorneAD

src/Gradus.jl

@@ -78,7 +78,12 @@ export AbstractCoronaModel,
 include("metrics/metrics.jl")
 
 export BoyerLindquistAD,
-    BoyerLindquistFO, JohannsenAD, MorrisThorneAD, KerrRefractiveAD, DilatonAxionAD
+    BoyerLindquistFO,
+    JohannsenAD,
+    JohannsenPsaltisAD,
+    MorrisThorneAD,
+    KerrRefractiveAD,
+    DilatonAxionAD
 
 # downstream modules and work
 include("special-radii.jl")

src/metrics/johannsen-psaltis-ad.jl

@@ -0,0 +1,48 @@
+module __JohannsenPsaltisAD
+using ..StaticArrays
+
+@inline h(r, M, ϵ3, Σ) = ϵ3 * M^3 * r / Σ^2
+@inline Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2
+@inline Δ(r, M, a) = r^2 - 2M * r + a^2
+
+@fastmath function metric_components(M, a, ϵ3, rθ)
+    (r, θ) = rθ
+
+    Σ₀ = Σ(r, a, θ)
+    h₀ = h(r, M, ϵ3, Σ₀)
+    sinθ2 = sin(θ)^2
+
+    tt = -(1 + h₀) * (1 - 2M * r / Σ₀)
+    rr = Σ₀ * (1 + h₀) / (Δ(r, M, a) + a^2 * sinθ2 * h₀)
+    θθ = Σ₀
+
+    term1 = sinθ2 * (r^2 + a^2 + 2a^2 * M * r * sinθ2 / Σ₀)
+    term2 = h₀ * a^2 * (Σ₀ + 2M * r) * sinθ2^2 / Σ₀
+    ϕϕ = term1 + term2
+
+    tϕ = -2a * M * r * sinθ2 * (1 + h₀) / Σ₀
+
+    @SVector [tt, rr, θθ, ϕϕ, tϕ]
+end
+
+end # module
+
+"""
+    struct JohannsenPsaltisAD{T} <: AbstractAutoDiffStaticAxisSymmetricParams{T}
+
+Johannsen and Psaltis 2011
+$(FIELDS)
+"""
+@with_kw struct JohannsenPsaltisAD{T} <: AbstractAutoDiffStaticAxisSymmetricParams{T}
+    @deftype T
+    "Black hole mass."
+    M = 1.0
+    "Black hole spin."
+    a = 0.0
+    "``\\epsilon_{3}`` deviation parameter."
+    ϵ3 = 0.0
+end
+
+GeodesicTracer.metric_components(m::JohannsenPsaltisAD{T}, rθ) where {T} =
+    __JohannsenPsaltisAD.metric_components(m.M, m.a, m.ϵ3, rθ)
+GradusBase.inner_radius(m::JohannsenPsaltisAD{T}) where {T} = m.M + √(m.M^2 - m.a^2)

src/metrics/metrics.jl

@@ -1,6 +1,7 @@
 include("boyer-lindquist-ad.jl")
 include("boyer-lindquist-fo.jl")
 include("johannsen-ad.jl")
+include("johannsen-psaltis-ad.jl")
 include("morris-thorne-ad.jl")
 include("kerr-refractive-ad.jl")
 include("dilaton-axion-ad.jl")

src/special-radii.jl

@@ -75,3 +75,59 @@ Event horizon radius, often equivalent to [`GradusBase.inner_radius`](@ref), how
 distinct, such that the latter may still be an arbitrary chart cutoff.
 """
 r_s(m::AbstractMetricParams{T}) where {T} = error("Not implemented for $(typeof(m)).")
+
+"""
+    event_horizon(m::AbstractMetricParams; select = last, resolution = 100, θε = 1e-7, rmax = 5.0)
+
+Utility function for helping plot an event horizon shape. Returns a tuple containing the `r`
+and `θ` vectors that solve
+
+```math
+    g_{t\\phi}^2 - g_{tt} g_{\\phi \\phi} = 0.
+```
+
+A `NaN` value in the `r` vector indicates no solution for that particular ``\\theta``, i.e.
+that the metric describes a naked singularity.
+
+Often the equation will have multiple roots, in which case the keyword argument `select` may be
+assigned to select the desired root.
+"""
+event_horizon(m::AbstractMetricParams{T}; kwargs...) where {T} =
+    error("Not implemented for $(typeof(m)).")
+
+function _event_horizon_condition(m, r, θ)
+    g = metric_components(m, (r, θ))
+    g[5]^2 - g[1] * g[4]
+end
+
+function event_horizon(
+    m::AbstractAutoDiffStaticAxisSymmetricParams{T};
+    select = maximum,
+    resolution::Int = 100,
+    θε::T = T(1e-7),
+    rmax = 5.0,
+) where {T}
+    θs = range(θε, 2π - θε, resolution)
+    rs = map(θs) do θ
+        f(r) = _event_horizon_condition(m, r, θ)
+        r = Roots.find_zeros(f, 0.0, rmax)
+        if isempty(r)
+            push!(r, NaN)
+        end
+        select(r)
+    end
+    rs, θs
+end
+
+function is_naked_singularity(
+    m::AbstractAutoDiffStaticAxisSymmetricParams{T};
+    resolution::Int = 100,
+    θε::T = T(1e-7),
+    rmax = 5.0,
+) where {T}
+    any(range(θε, 2π - θε, resolution)) do θ
+        f(r) = _event_horizon_condition(m, r, θ)
+        r = Roots.find_zeros(f, 0.0, rmax)
+        isempty(r)
+    end
+end