From e07b61e5cb42566aae98d43b035181d56e3c518d Mon Sep 17 00:00:00 2001 From: fjebaker Date: Sun, 14 Apr 2024 19:19:11 +0100 Subject: [PATCH] feat: implemented NoZ metric --- src/Gradus.jl | 1 + src/geometry/discs/thin-disc.jl | 41 ++++++++++++- src/line-profiles.jl | 3 +- src/metrics/noz-metric.jl | 100 ++++++++++++++++++++++++++++++++ 4 files changed, 143 insertions(+), 2 deletions(-) create mode 100644 src/metrics/noz-metric.jl diff --git a/src/Gradus.jl b/src/Gradus.jl index 8af1b769..b3123db1 100644 --- a/src/Gradus.jl +++ b/src/Gradus.jl @@ -435,6 +435,7 @@ include("metrics/johannsen-ad.jl") include("metrics/johannsen-psaltis-ad.jl") include("metrics/morris-thorne-ad.jl") include("metrics/kerr-refractive-ad.jl") +include("metrics/noz-metric.jl") include("metrics/dilaton-axion-ad.jl") include("metrics/bumblebee-ad.jl") include("metrics/kerr-newman-ad.jl") diff --git a/src/geometry/discs/thin-disc.jl b/src/geometry/discs/thin-disc.jl index 2951d2f8..f4e7af56 100644 --- a/src/geometry/discs/thin-disc.jl +++ b/src/geometry/discs/thin-disc.jl @@ -25,4 +25,43 @@ function distance_to_disc(d::ThinDisc, x4; gtol) _spinaxis_project(x4, signed = false) - _gtol_error(gtol, x4) end -export ThinDisc +""" + struct WarpedThinDisc{T} <: AbstractAccretionDisc{T} + WarpedThinDisc(inner_radius::T, outer_radius::T, f::F) + +Similar to [`ThinDisc`](@ref), however the scale height as a function of the radial coordinate may be specified by an arbitrary function. The function should return the ``h = \\cos(\\theta)`` signed scaled height (c.f. [`ThickDisc`](@ref) which is unsigned). + +The function should be of the form + +```julia +function scale_height(ρ) + # ... +end +``` +""" +struct WarpedThinDisc{T,F} <: AbstractAccretionDisc{T} + f::F + inner_radius::T + outer_radius::T +end + +WarpedThinDisc(f; inner_radius = 0.0, outer_radius = 500.0) = WarpedThinDisc(f, inner_radius, outer_radius) + +optical_property(::Type{<:WarpedThinDisc}) = OpticallyThin() + +inner_radius(disc::WarpedThinDisc) = disc.inner_radius + +function distance_to_disc(d::WarpedThinDisc, x4; gtol) + ρ = _equatorial_project(x4) + if ρ < d.inner_radius || ρ > d.outer_radius + return one(eltype(x4)) + end + + h = d.f(ρ) + γ = _spinaxis_project(x4, signed = true) + + abs(h - γ) - _gtol_error(gtol, x4) +end + +export ThinDisc, WarpedThinDisc + diff --git a/src/line-profiles.jl b/src/line-profiles.jl index b74268b8..631385c4 100644 --- a/src/line-profiles.jl +++ b/src/line-profiles.jl @@ -88,6 +88,7 @@ function lineprofile( # this 5maxrₑ is a heuristic that is insulting # todo: make it friendly plane = PolarPlane(GeometricGrid(); Nr = 450, Nθ = 1300, r_max = 5maxrₑ), + callback = domain_upper_hemisphere(), solver_args..., ) where {T} progress_bar = init_progress_bar("Lineprofile: ", trajectory_count(plane), verbose) @@ -101,7 +102,7 @@ function lineprofile( save_on = false, verbose = verbose, progress_bar = progress_bar, - callback = domain_upper_hemisphere(), + callback = callback, ensemble = EnsembleEndpointThreads(), solver_args..., ) diff --git a/src/metrics/noz-metric.jl b/src/metrics/noz-metric.jl new file mode 100644 index 00000000..75ff50b9 --- /dev/null +++ b/src/metrics/noz-metric.jl @@ -0,0 +1,100 @@ + +module __NoZMetric +using ..StaticArrays +using ..MuladdMacro + +@muladd @fastmath begin + + epsilon(M, a, ϵ, y) = ϵ * M * a * y + + # the way this function must be defined is a little complex + # but helps with type-stability + function metric_components(M, a, ϵ, rθ) + (r, θ) = rθ + sinθ2 = sin(θ)^2 + + y = cos(θ) + _epsilon = epsilon(M, a, ϵ, y) + + tt = + -1 + ((2M * r * (r^(2) + a^(2) * y^(2))) / + ((r^(2) + a^(2) * y^(2))^(2) + (r^(2) - 2M * r + a^(2) * y^(2)) * _epsilon)) + ϕϕ = + ( + (1 - y^2) * (r^2 + a^2 * y^2 + _epsilon) * ( + r^4 + + a^4 * y^2 + + r^2 * (a^2 + a^2 * y^2 + _epsilon) + + a^2 * _epsilon + + 2M * r * (a^2 - a^2 * y^2 - _epsilon) + ) + ) / ((r^2 + a^2 * y^2)^2 + (r^2 - 2M * r + a^2 * y^2) * _epsilon) + rr = (r^(2) + a^(2) * y^(2) + _epsilon) / (r^(2) - 2M * r + a^(2)) + yy = (r^(2) + a^(2) * y^(2) + _epsilon) / (1 - y^(2)) + + tϕ = + -(2M * r * a*(1 - y^(2)) * (r^(2) + a^(2) * y^(2) + _epsilon)) / + ((r^(2) + a^(2) * y^(2))^(2) + (r^(2) - 2M * r + a^(2) * y^(2)) * _epsilon) + + # include the coordinate transformation factor + # y = cos(θ) --> dy^2 = dθ^2 * sin(θ)^2 + @SVector [tt, rr, yy * sinθ2, ϕϕ, tϕ] + end +end + +end # module + +# new structure for our spacetime +""" + struct NoZMetric +""" +@with_kw struct NoZMetric{T} <: AbstractStaticAxisSymmetric{T} + @deftype T + "Black hole mass." + M = 1.0 + "Black hole spin." + a = 0.0 + "Deviation parameter" + ϵ = 0.0 +end + +# implementation +metric_components(m::NoZMetric, rθ) = __NoZMetric.metric_components(m.M, m.a, m.ϵ, rθ) +inner_radius(m::NoZMetric) = m.M + √(m.M^2 - m.a^2) + +function _solve_orbit_θ(m, r) + function _objective(θ) + rθ = SVector(r, θ) + _, J = Gradus.metric_jacobian(m, rθ) + ∂rg = J[:, 1] + ∂θg = J[:, 2] + Ωϕ = Gradus.CircularOrbits._Ω_analytic(∂rg, false) + + ∂θg[1] + 2 * ∂θg[5] * Ωϕ + ∂θg[4] * Ωϕ^2 + end + + Gradus.Roots.find_zero(_objective, π / 2) +end + +function isco(m::NoZMetric{T}) where {T} + kerr_isco = isco(KerrMetric(M = m.M, a = m.a)) + rs = range(inner_radius(m), kerr_isco + 1.0, 100) + thetas = map(rs) do r + try + _solve_orbit_θ(m, r) + catch + zero(T) + end + end + + # TODO: this is such a hack + replace!(thetas, zero(T) => thetas[findfirst(!=(zero(T)), thetas)]) + + interp = _make_interpolation(rs, thetas) + + dE(r) = ForwardDiff.derivative(x -> CircularOrbits.energy(m, SVector(x, interp(x))), r) + # optimize from the Kerr equivalent metric + Roots.find_zero(dE, kerr_isco) +end + +export NoZMetric