diff --git a/src/Gradus.jl b/src/Gradus.jl index 4e45bcd6..a574f74c 100644 --- a/src/Gradus.jl +++ b/src/Gradus.jl @@ -162,12 +162,14 @@ include("accretion-geometry/discs.jl") include("accretion-geometry/meshes.jl") include("accretion-geometry/bootstrap.jl") -include("cunningham-transfer.jl") +include("transfer-functions/cunningham-transfer-functions.jl") +include("transfer-functions/integration.jl") +include("transfer-functions/precision-solvers.jl") include("corona-to-disc/sky-geometry.jl") include("corona-to-disc/corona-models.jl") include("corona-to-disc/disc-profiles.jl") -include("corona-to-disc/transfer-functions.jl") +include("transfer-functions/transfer-functions-2d.jl") include("corona-to-disc/flux-calculations.jl") include("metrics/boyer-lindquist-ad.jl") diff --git a/src/const-point-functions.jl b/src/const-point-functions.jl index c923b338..54534607 100644 --- a/src/const-point-functions.jl +++ b/src/const-point-functions.jl @@ -24,8 +24,10 @@ const filter_early_term = A [`FilterPointFunction`](@ref) that filters geodesics which intersected with the accretion disc. Default: `NaN`. """ -const filter_intersected = - FilterPointFunction((m, gp, max_time; kwargs...) -> gp.status == StatusCodes.IntersectedWithGeometry, NaN) +const filter_intersected = FilterPointFunction( + (m, gp, max_time; kwargs...) -> gp.status == StatusCodes.IntersectedWithGeometry, + NaN, +) """ affine_time(m::AbstractMetricParams, gp::AbstractGeodesicPoint, max_time) diff --git a/src/corona-to-disc/disc-profiles.jl b/src/corona-to-disc/disc-profiles.jl index 792ef156..7360c369 100644 --- a/src/corona-to-disc/disc-profiles.jl +++ b/src/corona-to-disc/disc-profiles.jl @@ -106,7 +106,11 @@ function VoronoiDiscProfile( d::AbstractAccretionDisc{T}, simsols::SciMLBase.EnsembleSolution{T}, ) where {T} - VoronoiDiscProfile(m, d, filter(i -> i.prob.p.status == StatusCodes.IntersectedWithGeometry, simsols.u)) + VoronoiDiscProfile( + m, + d, + filter(i -> i.prob.p.status == StatusCodes.IntersectedWithGeometry, simsols.u), + ) end @noinline function findindex(vdp::VoronoiDiscProfile{D,V}, p) where {D,V} diff --git a/src/cunningham-transfer.jl b/src/cunningham-transfer.jl deleted file mode 100644 index 506a3411..00000000 --- a/src/cunningham-transfer.jl +++ /dev/null @@ -1,564 +0,0 @@ -@with_kw struct CunninghamTransferFunction{T} - rₑ::T - gs::Vector{T} - gstar::Vector{T} - f::Vector{T} -end - -""" - integrate_single_geodesic(m::AbstractMetricParams, u, d::AbstractAccretionDisc, rₒ, θₒ; kwargs...) - -Integrate a single geodesic with impact parameters calculated via - -```math -\\begin{align} - \\alpha &= r_\\text{o} \\cos (\\theta_\\text{o}), \\ - \\beta &= r_\\text{o} \\sin (\\theta_\\text{o}). -\\end{align} -``` - -Returns an [AbstractGeodesicPoint](@ref), depending on `m`. - -""" -function integrate_single_geodesic( - m::AbstractMetricParams, - u, - d::AbstractAccretionDisc, - rₒ, - θₒ; - max_time = 2e3, - kwargs..., -) - α = rₒ * cos(θₒ) - β = rₒ * sin(θₒ) - v = map_impact_parameters(m, u, α, β) - sol = tracegeodesics(m, u, v, d, (0.0, max_time); save_on = false, kwargs...) - getgeodesicpoint(m, sol) -end - -function find_offset_for_radius( - m::AbstractMetricParams, - u, - d::AbstractAccretionDisc, - radius, - θₒ; - zero_atol = 1e-7, - offset_max = 20.0, - kwargs..., -) - f = r -> begin - gp = integrate_single_geodesic(m, u, d, r, θₒ; kwargs...) - r = if gp.status == StatusCodes.IntersectedWithGeometry - gp.u2[2] * sin(gp.u2[3]) - else - 0.0 - end - radius - r - end - r = Roots.find_zero(f, (0.0, offset_max); atol = zero_atol) - if !isapprox(f(r), 0.0, atol = 10 * zero_atol) - return NaN - end - r -end - -function _find_extremal_redshift_with_guess( - m, - u, - d, - rₑ, - g_guess, - θ_guess; - redshift_pf = Gradus.ConstPointFunctions.redshift, - δθ = 1e-1, - minimal = true, - zero_atol = 1e-7, - offset_max = 20.0, - kwargs..., -) - f = - θ -> begin - r_offset = find_offset_for_radius( - m, - u, - d, - rₑ, - θ; - zero_atol = zero_atol, - offset_max = offset_max, - kwargs..., - ) - if !isnan(r_offset) - gp = integrate_single_geodesic(m, u, d, r_offset, θ; kwargs...) - g = redshift_pf(m, gp, gp.t2) - minimal ? g : -g - else - NaN - # minimal ? g_guess : -g_guess - end - end - - res = optimize(f, θ_guess - δθ, θ_guess + δθ) - θopt = Optim.minimizer(res) - ( - res, - find_offset_for_radius( - m, - u, - d, - rₑ, - θopt; - zero_atol = zero_atol, - offset_max = offset_max, - kwargs..., - ), - ) -end - -function impact_parameters_for_radius!( - α, - β, - m::AbstractMetricParams, - u::AbstractVector{T}, - d::AbstractAccretionDisc, - radius; - kwargs..., -) where {T} - if size(α) != size(β) - throw(DimensionMismatch("α, β must have the same dimensions and size.")) - end - θs = range(0, 2π, length(α)) - @inbounds @threads for i in eachindex(θs) - θ = θs[i] - r = find_offset_for_radius(m, u, d, radius, θ; kwargs...) - α[i] = r * cos(θ) - β[i] = r * sin(θ) - end - (α, β) -end - -function impact_parameters_for_radius( - m::AbstractMetricParams, - u::AbstractVector{T}, - d::AbstractAccretionDisc, - radius; - N = 500, - kwargs..., -) where {T} - α = zeros(T, N) - β = zeros(T, N) - impact_parameters_for_radius!(α, β, m, u, d, radius; kwargs...) - (filter(!isnan, α), filter(!isnan, β)) -end - -function redshift_ratio!( - gs, - m::AbstractMetricParams, - u, - d::AbstractAccretionGeometry, - max_time, - αs, - βs; - redshift_pf = Gradus.ConstPointFunctions.redshift, - kwargs..., -) - if size(gs) != size(αs) - throw(DimensionMismatch("αs, gs must have the same dimensions and size.")) - end - - velfunc(i) = map_impact_parameters(m, u, αs[i], βs[i]) - simsols = tracegeodesics( - m, - u, - velfunc, - d, - (0.0, max_time); - trajectories = length(αs), - save_on = false, - kwargs..., - ) - - @inbounds for (i, sol) in enumerate(simsols.u) - gp = getgeodesicpoint(m, sol) - g = redshift_pf(m, gp, max_time) - gs[i] = g - end - gs -end - -function redshift_ratio( - m::AbstractMetricParams, - u, - d::AbstractAccretionGeometry, - max_time, - αs::AbstractVector{T}, - βs; - kwargs..., -) where {T} - gs = zeros(T, length(αs)) - redshift_ratio!(gs, m, u, d, max_time, αs, βs; kwargs...) - gs -end - -function jacobian_∂αβ_∂gr!( - Js, - m, - u, - d, - max_time, - gs, - gmin, - gmax, - αs, - βs; - order = 5, - redshift_pf = Gradus.ConstPointFunctions.redshift, - kwargs..., -) - if size(gs) != size(Js) || size(Js) != size(αs) || size(αs) != size(βs) - throw( - DimensionMismatch("αs, βs, gs, and Js must have the same dimensions and size."), - ) - end - - f = - ((α, β),) -> begin - v = map_impact_parameters(m, u, α, β) - sol = - tracegeodesics(m, u, v, d, (0.0, max_time); save_on = false, kwargs...) - gp = getgeodesicpoint(m, sol) - g = redshift_pf(m, gp, 2000.0) - # return r and g* - @SVector [gp.u2[2], gstar(g, gmin, gmax)] - end - - # choice between FiniteDifferences and FiniteDiff is tricky - # since FiniteDiff is so much faster, but seems to yield really bad jacobians - # for this specific problem, so instead stenciling with a given order - cfdm = FiniteDifferences.central_fdm(order, 1) - @inbounds @threads for i in eachindex(αs) - x = @SVector [αs[i], βs[i]] - j = FiniteDifferences.jacobian(cfdm, f, x) |> first - Js[i] = abs(inv(det(j))) - end - Js -end - -function jacobian_∂αβ_∂gr( - m, - u, - d, - max_time, - gs::AbstractArray{T}, - gmin, - gmax, - αs::AbstractArray{T}, - βs; - kwargs..., -) where {T} - Js = zeros(T, size(αs)) - jacobian_∂αβ_∂gr!(Js, m, u, d, max_time, gs, gmin, gmax, αs, βs; kwargs...) - Js -end - -gstar(g::AbstractArray) = gstar(g, extrema(g)...) -function gstar(g, gmin, gmax) - Δg = gmax - gmin - @. (g - gmin) / Δg -end - -g_to_gstar(args...) = gstar(args...) -function gstar_to_g(gs, gmin, gmax) - Δg = gmax - gmin - @. gs * Δg + gmin -end - -@muladd cunningham_transfer_function( - rₑ::Number, - gs::AbstractArray, - gstars::AbstractArray, - js::AbstractArray, -) = @. (1 / (π * rₑ)) * gs * √(gstars * (1 - gstars)) * js - -function cunningham_transfer_function!( - αs, - βs, - Js::AbstractVector{T}, - m::AbstractMetricParams, - u, - d::AbstractAccretionGeometry, - rₑ, - max_time; - finite_diff_order = 6, - redshift_pf::PointFunction = ConstPointFunctions.redshift, - offset_max = 20.0, - zero_atol = 1e-7, - tracer_kwargs..., -) where {T} - impact_parameters_for_radius!( - αs, - βs, - m, - u, - d, - rₑ; - offset_max = offset_max, - zero_atol = zero_atol, - tracer_kwargs..., - ) - # filter and trim arrays - mask = @. !(isnan(αs) | isnan(βs)) - _αs = @views αs[mask] - _βs = @views βs[mask] - - N = length(_αs) - @assert N > 0 "No valid impact parameters found (all are NaN)." - _Js = @views Js[1:N] - - gs = zeros(T, N) - redshift_ratio!( - gs, - m, - u, - d, - max_time, - _αs, - _βs; - redshift_pf = redshift_pf, - tracer_kwargs..., - ) - - gmin, gmax = extrema(gs) - jacobian_∂αβ_∂gr!( - _Js, - m, - u, - d, - max_time, - gs, - gmin, - gmax, - _αs, - _βs; - order = finite_diff_order, - redshift_pf = redshift_pf, - tracer_kwargs..., - ) - - gstars = gstar(gs, gmin, gmax) - f = cunningham_transfer_function(rₑ, gs, gstars, _Js) - - # package and return - CunninghamTransferFunction(rₑ, gs, gstars, f) -end - -function cunningham_transfer_function( - m::AbstractMetricParams{T}, - u, - d::AbstractAccretionGeometry, - rₑ, - max_time; - num_points = 1000, - kwargs..., -) where {T} - α = zeros(T, num_points) - β = zeros(T, num_points) - Js = zeros(T, num_points) - cunningham_transfer_function!(α, β, Js, m, u, d, rₑ, max_time; kwargs...) -end - -function _split_branches(gstars::AbstractArray{T}, f::AbstractArray{T}) where {T} - decreasing = true - gprev = 1.0 - upper = Tuple{T,T}[] - lower = Tuple{T,T}[] - for (i, g) in enumerate(gstars) - if gprev > g - decreasing = true - else - decreasing = false - end - if decreasing - push!(upper, (g, f[i])) - else - push!(lower, (g, f[i])) - end - gprev = g - end - lower, upper -end - -struct InterpolatedCunninghamTransferBranches{T,I} - lower::I - upper::I - radius::T - g_extrema::Tuple{T,T} - g_limits::Tuple{T,T} -end - -_gs_in_limits(ictb, gs) = ictb.g_limits[2] > gs > ictb.g_limits[1] - -function _make_sorted_interpolation(g, f) - I = sortperm(g) - _g = @inbounds g[I] - _f = @inbounds f[I] - # Interpolations.deduplicate_knots!(_g) - # linear_interpolation(_g, _f, extrapolation_bc = NaN) - DataInterpolations.LinearInterpolation(_f, _g) -end - -function _interpolate_branches(ctf::CunninghamTransferFunction; offset = 1e-4) - # avoid extrema - mask = @. (ctf.gstar > offset) & (ctf.gstar < 1 - offset) - - gstars = @inbounds @views ctf.gstar[mask] - f = @inbounds @views ctf.f[mask] - - lower, upper = _split_branches(gstars, f) - - gs_lower = first.(lower) - gs_upper = first.(upper) - # valid interval - gs_min_limit = max(minimum(gs_lower), minimum(gs_upper)) - gs_max_limit = min(maximum(gs_lower), maximum(gs_upper)) - - # interpolate - f1 = _make_sorted_interpolation(gs_lower, last.(lower)) - f2 = _make_sorted_interpolation(gs_upper, last.(upper)) - - InterpolatedCunninghamTransferBranches( - f1, - f2, - ctf.rₑ, - extrema(ctf.gs), - (gs_min_limit, gs_max_limit), - ) -end - -@muladd _cunning_integrand(f, g, rₑ, gs, gmin, gmax) = - f * g^3 * π * rₑ / (√(gs * (1 - gs)) * (gmax - gmin)) - -function integrate_transfer_functions( - ε, - ictbs::Vector{<:InterpolatedCunninghamTransferBranches{T1}}, - bins, -) where {T1} - # global min and max - ggmin = maximum(i -> first(i.g_limits), ictbs) - ggmax = minimum(i -> last(i.g_limits), ictbs) - - radii = map(i -> i.radius, ictbs) - - r_low, r_high = extrema(radii) - - segbuf = QuadGK.alloc_segbuf(T1, eltype(bins)) - y = map(eachindex(@view(bins[1:end-1]))) do index - # bin edges - bin_low = bins[index] - bin_high = bins[index+1] - # calculate area under bin - integrated_bin_for_radii = - _areas_under_transfer_functions(ε, ictbs, bin_low, bin_high, ggmin, ggmax) - - # interpolate across different radii - intp = DataInterpolations.LinearInterpolation(integrated_bin_for_radii, radii) - # use (smooth) interpolation to integrate fully - res::T1, _ = quadgk(intp, r_low, r_high; segbuf = segbuf, order = 4) - res - end - - (bins, y) -end - -@muladd function _cunningham_line_profile_integrand(ictb) - gs -> begin - if !Gradus._gs_in_limits(ictb, gs) - 0.0 - else - gmin, gmax = ictb.g_extrema - g = gstar_to_g(gs, gmin, gmax) - w = g^3 / √(gs * (1 - gs)) - f = ictb.lower(gs) + ictb.upper(gs) - w * f / (gmax - gmin) - end - end -end - -function _calculate_interpolated_transfer_branches( - m::AbstractMetricParams{T}, - u, - d, - radii; - num_points = 100, - verbose = false, - offset = 1e-7, - kwargs..., -) where {T} - # pre-allocate arrays - αs = zeros(T, num_points) - βs = zeros(T, num_points) - Js = zeros(T, num_points) - - progress_bar = init_progress_bar("Transfer functions:", length(radii), verbose) - - # IILF for calculating the interpolated branches - 𝔉 = - rₑ -> begin - # this feels like such a bad practice - # calling GC on young objects to clear the temporarily allocated memory - # but we get a significant speedup - GC.gc(false) - ctf = cunningham_transfer_function!( - αs, - βs, - Js, - m, - u, - d, - rₑ, - 2000.0; - offset_max = 2rₑ + 20.0, - kwargs..., - ) - ProgressMeter.next!(progress_bar) - _interpolate_branches(ctf; offset = offset) - end - - # calculate interpolated transfer functions for each emission radius - map(𝔉, radii) -end - -function _areas_under_transfer_functions( - ε, - ictbs::Vector{<:InterpolatedCunninghamTransferBranches{T1}}, - bin_low::T2, - bin_high::T2, - ggmin, - ggmax, -) where {T1,T2} - segbuf = QuadGK.alloc_segbuf(T1, T2) - map(ictbs) do ictb - 𝒮 = _cunningham_line_profile_integrand(ictb) - Sg = g -> begin - gs = gstar(g, ictb.g_extrema...) - if ggmin < gs < ggmax - 𝒮(gs) - else - 0.0 - end - end - res, _ = quadgk(Sg, bin_low, bin_high; order = 2, segbuf = segbuf) - r = ictb.radius - res * r * ε(r) - end -end - -export impact_parameters_for_radius, - redshift_ratio, - jacobian_∂αβ_∂gr, - gstar, - g_to_gstar, - gstar_to_g, - cunningham_transfer_function, - CunninghamTransferFunction, - InterpolatedCunninghamTransferBranches diff --git a/src/line-profiles.jl b/src/line-profiles.jl index ab0532cc..8d09146c 100644 --- a/src/line-profiles.jl +++ b/src/line-profiles.jl @@ -3,6 +3,7 @@ struct CunninghamLineProfile <: AbstractLineProfileAlgorithm end struct BinnedLineProfile <: AbstractLineProfileAlgorithm end @inline function lineprofile( + bins, ε::Function, m::AbstractMetricParams, u, @@ -10,51 +11,28 @@ struct BinnedLineProfile <: AbstractLineProfileAlgorithm end algorithm::AbstractLineProfileAlgorithm = CunninghamLineProfile(), kwargs..., ) - lineprofile(algorithm, m, u, d, ε; kwargs...) -end - -function _change_interval(f, a, b) - α = (b - a) / 2 - β = (a + b) / 2 - (α, x -> f(α * x + β)) + lineprofile(bins, algorithm, m, u, d, ε; kwargs...) end function lineprofile( + bins, ::CunninghamLineProfile, m::AbstractMetricParams{T}, u, d::AbstractAccretionGeometry, ε; - num_points = 100, - min_re = isco(m) + 1e-2, # delta to avoid numerical instabilities - max_re = 20, - num_re = 100, - bins = range(0.0, 1.5, 100), + minrₑ = isco(m) + 1e-2, # delta to avoid numerical instabilities + maxrₑ = 50, + numrₑ = 100, verbose = false, - offset = 1e-7, + Ng✶ = 500, + h = 2e-8, kwargs..., ) where {T} - # this is just a placeholder: desire a distribution that favours - # small radii over large radii, and this one does that quite well - # radii here are emission radii rₑ - radii = exp.(range(log(1), log(1000), num_re)) - _max_radii = maximum(radii) - # rescale inplace - @. radii = (radii / _max_radii) * (max_re - min_re) + min_re - - ictbs = _calculate_interpolated_transfer_branches( - m, - u, - d, - radii; - num_points = num_points, - verbose = verbose, - offset = offset, - kwargs..., - ) - - integrate_transfer_functions(ε, ictbs, bins) + radii = weighted_rₑ_grid(minrₑ, maxrₑ, numrₑ) + itfs = interpolated_transfer_branches(m, u, d, radii; verbose = verbose, kwargs...) + flux = integrate_drdg✶(ε, itfs, radii, bins; h = h, Ng✶ = Ng✶) + bins, flux end - export AbstractLineProfileAlgorithm, BinnedLineProfile, CunninghamLineProfile, lineprofile diff --git a/src/tracing/callbacks.jl b/src/tracing/callbacks.jl index ea3c5d94..2eaa921a 100644 --- a/src/tracing/callbacks.jl +++ b/src/tracing/callbacks.jl @@ -41,7 +41,10 @@ end # predefined callbacks function domain_upper_hemisphere(δ = 1e-3) - DiscreteCallback((u, t, integrator) -> u[2] * cos(u[3]) < δ, terminate_with_status!(StatusCodes.OutOfDomain)) + DiscreteCallback( + (u, t, integrator) -> u[2] * cos(u[3]) < δ, + terminate_with_status!(StatusCodes.OutOfDomain), + ) end export domain_upper_hemisphere diff --git a/src/tracing/method-implementations/first-order.jl b/src/tracing/method-implementations/first-order.jl index 57c8e4da..1b1b56ed 100644 --- a/src/tracing/method-implementations/first-order.jl +++ b/src/tracing/method-implementations/first-order.jl @@ -81,9 +81,8 @@ mutable struct FirstOrderIntegrationParameters{T} <: AbstractIntegrationParamete status::StatusCodes.T end -make_parameters(L, Q, sign_θ, ::Type{T}) where {T} = FirstOrderIntegrationParameters{T}( - L, Q, -1, sign_θ, [0.0, 0.0], StatusCodes.NoStatus -) +make_parameters(L, Q, sign_θ, ::Type{T}) where {T} = + FirstOrderIntegrationParameters{T}(L, Q, -1, sign_θ, [0.0, 0.0], StatusCodes.NoStatus) function integrator_problem( m::AbstractFirstOrderMetricParams{T}, diff --git a/src/transfer-functions/cunningham-transfer-functions.jl b/src/transfer-functions/cunningham-transfer-functions.jl new file mode 100644 index 00000000..aa78f39d --- /dev/null +++ b/src/transfer-functions/cunningham-transfer-functions.jl @@ -0,0 +1,185 @@ +@with_kw struct CunninghamTransferFunction{T} + "``g^\\ast`` values." + g✶::Vector{T} + "Transfer function data." + f::Vector{T} + gmin::T + gmax::T + "Emission radius." + rₑ::T +end + +@with_kw struct InterpolatedCunninghamTransferFunction{T,U,L} + upper_f::U + lower_f::L + gmin::T + gmax::T + rₑ::T +end + +function _make_sorted_interpolation(g, f) + I = sortperm(g) + _g = @inbounds g[I] + _f = @inbounds f[I] + DataInterpolations.LinearInterpolation(_f, _g) +end + +function splitbranches(ctf::CunninghamTransferFunction{T}) where {T} + upper_f = T[] + lower_f = T[] + upper_g✶ = T[] + lower_g✶ = T[] + N = (length(ctf.f) ÷ 2) + 1 + sizehint!(upper_f, N) + sizehint!(lower_f, N) + sizehint!(upper_g✶, N) + sizehint!(lower_g✶, N) + + decreasing = true + g✶previous = 0.0 + for (i, g✶) in enumerate(ctf.g✶) + if g✶previous > g✶ + decreasing = true + else + decreasing = false + end + if decreasing + push!(lower_f, ctf.f[i]) + push!(lower_g✶, g✶) + else + push!(upper_f, ctf.f[i]) + push!(upper_g✶, g✶) + end + g✶previous = g✶ + end + + (lower_g✶, lower_f, upper_g✶, upper_f) +end + +function interpolate_transfer_function(ctf::CunninghamTransferFunction{T}) where {T} + (lower_g✶, lower_f, upper_g✶, upper_f) = splitbranches(ctf) + lower_branch = _make_sorted_interpolation(lower_g✶, lower_f) + upper_branch = _make_sorted_interpolation(upper_g✶, upper_f) + InterpolatedCunninghamTransferFunction( + upper_branch, + lower_branch, + ctf.gmin, + ctf.gmax, + ctf.rₑ, + ) +end + +@muladd function _calculate_transfer_function(rₑ, g, g✶, J) + @. (1 / (π * rₑ)) * g * √(g✶ * (1 - g✶)) * J +end + +function cunningham_transfer_function( + m::AbstractMetricParams{T}, + u, + d, + rₑ; + max_time = 2e3, + diff_order = 5, + redshift_pf = ConstPointFunctions.redshift, + offset_max = 20.0, + zero_atol = 1e-7, + N = 100, + tracer_kwargs..., +) where {T} + Js = zeros(T, N) + gs = zeros(T, N) + + θs = range(π / 2, 2π + π / 2, N) + @inbounds for i in eachindex(θs) + θ = θs[i] + r, gp = find_offset_for_radius( + m, + u, + d, + rₑ, + θ; + zero_atol = zero_atol, + offset_max = offset_max, + max_time = max_time, + tracer_kwargs..., + ) + α = r * cos(θ) + β = r * sin(θ) + g = redshift_pf(m, gp, max_time) + gs[i] = g + Js[i] = jacobian_∂αβ_∂gr( + m, + u, + d, + α, + β, + max_time; + diff_order = diff_order, + redshift_pf = redshift_pf, + tracer_kwargs..., + ) + end + + gmin, gmax = infer_extremal(gs, θs, π, 2π) + # convert from ∂g to ∂g✶ + @. Js = (gmax - gmin) * Js + @inbounds for i in eachindex(Js) + g✶ = g_to_g✶(gs[i], gmin, gmax) + # Js is now storing f + Js[i] = _calculate_transfer_function(rₑ, gs[i], g✶, Js[i]) + # gs is now storing g✶ + gs[i] = g✶ + end + + CunninghamTransferFunction(gs, Js, gmin, gmax, rₑ) +end + +function infer_extremal(y, x, x0, x1) + _, y0 = interpolate_extremal(y, x, x0) + _, y1 = interpolate_extremal(y, x, x1) + if y0 > y1 + return y1, y0 + else + return y0, y1 + end +end + +function interpolate_extremal(y, x, x0) + interp = DataInterpolations.CubicSpline(y, x) + ∂(f) = x -> ForwardDiff.derivative(f, x) + x̄ = find_zero(∂(interp), x0) + x̄, interp(x̄) +end + +function interpolated_transfer_branches( + m::AbstractMetricParams{T}, + u, + d, + radii; + verbose = false, + kwargs..., +) where {T} + # IILF for calculating the interpolated branches + 𝔉 = + rₑ -> begin + ctf = cunningham_transfer_function( + m, + u, + d, + rₑ, + ; + offset_max = 2rₑ + 20.0, + kwargs..., + ) + interpolate_transfer_function(ctf) + end + + # calculate interpolated transfer functions for each emission radius + ThreadsX.map(𝔉, radii) +end + +export CunninghamTransferFunction, +InterpolatedCunninghamTransferFunction, +splitbranches, +interpolate_transfer_function, +cunningham_transfer_function \ No newline at end of file diff --git a/src/transfer-functions/integration.jl b/src/transfer-functions/integration.jl new file mode 100644 index 00000000..8909557d --- /dev/null +++ b/src/transfer-functions/integration.jl @@ -0,0 +1,152 @@ +g_to_g✶(g, gmin, gmax) = @. (g - gmin) / (gmax - gmin) +g✶_to_g(g✶, gmin, gmax) = @. (gmax - gmin) * g✶ + gmin + +function interpolate_over_radii(itfs, g✶_grid) + radii = map(itf -> itf.rₑ, itfs) + + if !issorted(radii) + I = sortperm(radii) + radii = radii[I] + itfs = itfs[I] + end + + fr_interp = map(g✶_grid) do g✶ + f = map(itfs) do itf + itf.lower_f(g✶) + itf.upper_f(g✶) + end + DataInterpolations.LinearInterpolation(f, radii) + end + + gmin_interp = DataInterpolations.LinearInterpolation(map(itf -> itf.gmin, itfs), radii) + gmax_interp = DataInterpolations.LinearInterpolation(map(itf -> itf.gmax, itfs), radii) + fr_interp, gmin_interp, gmax_interp +end + +function integrate_edge(S, h, lim, gmin, gmax, low::Bool) + if low + gh = g✶_to_g(h, gmin, gmax) + a = √gh - √lim + else + gh = g✶_to_g(1 - h, gmin, gmax) + a = √lim - √gh + end + 2 * S(gh) * √h * a * (gmax - gmin) +end + +function _cunningham_integrand(f, g, gs, gmin, gmax) + f * π * (g)^3 / (√(gs * (1 - gs)) * (gmax - gmin)) +end + +function integrate_bin(S, gmin, gmax, lo, hi; h = 2e-8) + glo = clamp(lo, gmin, gmax) + ghi = clamp(hi, gmin, gmax) + + lum = 0.0 + if glo == ghi + return lum + end + + g✶lo = g_to_g✶(lo, gmin, gmax) + g✶hi = g_to_g✶(hi, gmin, gmax) + + if g✶lo < h + if g✶hi > h + lum += integrate_edge(S, h, glo, gmin, gmax, true) + glo = g✶_to_g(h, gmin, gmax) + else + return lum # integrate_edge(S, h, glo, ghi, gmin, gmax, true) + end + end + + if g✶hi > 1 - h + if g✶lo < 1 - h + lum += integrate_edge(S, h, ghi, gmin, gmax, false) + ghi = g✶_to_g(1 - h, gmin, gmax) + else + return lum # integrate_edge(S, h, glo, ghi, gmin, gmax) + end + end + + res, _ = quadgk(S, glo, ghi) + lum += res + return lum +end + +function _wrap_cunningham_interpolations(fr_interp, gmin, gmax, r, g✶_grid) + f = map(fr_interp) do interp + interp(r) + end + S = DataInterpolations.LinearInterpolation(f, g✶_grid) + g -> begin + g✶ = g_to_g✶(g, gmin, gmax) + _cunningham_integrand(S(g✶), g, g✶, gmin, gmax) + end +end + +function weighted_rₑ_grid(min, max, N) + Iterators.reverse(inv(r) for r in range(1 / max, 1 / min, N)) +end + +function _build_g✶_grid(Ng✶, h) + g✶low = h + g✶high = 1 - h + map(1:Ng✶) do i + g✶low + (g✶high - g✶low) / (Ng✶ - 1) * (i - 1) + end +end + +function integrate_drdg✶( + ε, + itfs::Vector{<:InterpolatedCunninghamTransferFunction{T}}, + radii, + g_grid; + Ng✶ = 10, + h = 1e-8, +) where {T} + # pre-allocate output + flux = zeros(T, length(g_grid)) + + # init params + minrₑ, maxrₑ = extrema(radii) + g✶_grid = _build_g✶_grid(Ng✶, h) + fr_interp, gmin_interp, gmax_interp = interpolate_over_radii(itfs, g✶_grid) + g_grid_view = @views g_grid[1:end-1] + + # build fine radial grid for trapezoidal integration + fine_rₑ_grid = weighted_rₑ_grid(minrₑ, maxrₑ, 1000) |> collect + @inbounds for (i, rₑ) in enumerate(fine_rₑ_grid) + gmin = gmin_interp(rₑ) + gmax = gmax_interp(rₑ) + S = _wrap_cunningham_interpolations(fr_interp, gmin, gmax, rₑ, g✶_grid) + + # trapezoidal integration weight + if i == 1 + Δrₑ = (fine_rₑ_grid[i+1]) - (rₑ) + elseif i == lastindex(fine_rₑ_grid) + Δrₑ = (rₑ) - (fine_rₑ_grid[i-1]) + else + Δrₑ = (fine_rₑ_grid[i+1]) - (fine_rₑ_grid[i-1]) + end + + weight = Δrₑ * rₑ * ε(rₑ) + + for j in eachindex(g_grid_view) + glo = g_grid[j] + ghi = g_grid[j+1] + flux[j] += integrate_bin(S, gmin, gmax, glo, ghi; h = h) * weight + end + end + + Σflux = zero(T) + @inbounds for i in eachindex(g_grid_view) + glo = g_grid[i] + ghi = g_grid[i+1] + ḡ = (ghi + glo) + flux[i] = flux[i] / ḡ + Σflux += flux[i] + end + @. flux = flux / Σflux + flux +end + +export integrate_drdg✶ diff --git a/src/transfer-functions/precision-solvers.jl b/src/transfer-functions/precision-solvers.jl new file mode 100644 index 00000000..4fd16e13 --- /dev/null +++ b/src/transfer-functions/precision-solvers.jl @@ -0,0 +1,185 @@ +""" + integrate_single_geodesic(m::AbstractMetricParams, u, d::AbstractAccretionDisc, rₒ, θₒ; kwargs...) + +Integrate a single geodesic with impact parameters calculated via + +```math +\\begin{align} + \\alpha &= r_\\text{o} \\cos (\\theta_\\text{o}), \\ + \\beta &= r_\\text{o} \\sin (\\theta_\\text{o}). +\\end{align} +``` + +Returns an [AbstractGeodesicPoint](@ref), depending on `m`. +""" +function integrate_single_geodesic( + m::AbstractMetricParams, + u, + d::AbstractAccretionDisc, + rₒ, + θₒ; + max_time = 2e4, + kwargs..., +) + α = rₒ * cos(θₒ) + β = rₒ * sin(θₒ) + v = map_impact_parameters(m, u, α, β) + sol = tracegeodesics(m, u, v, d, (0.0, max_time); save_on = false, kwargs...) + getgeodesicpoint(m, sol) +end + +""" + find_offset_for_radius( + m::AbstractMetricParams, + u, + d::AbstractAccretionDisc, + radius, + θₒ; + zero_atol = 1e-7, + offset_max = 20.0, + kwargs..., + ) + +Find the offset ``r_\\text{o}`` on the observer's image plane that gives impact parameters +```math +\\alpha = r_\\text{o} \\cos \\theta_\\text{o}, +\\quad \\text{and} \\quad +\\beta = r_\\text{o} \\sin \\theta_\\text{o}, +``` +that trace a geodesic to intersect the disc geometry at an emission radius ``r_\\text{e}``. + +Returns ``NaN`` if no offset exists. This may occur of the geometry is not present at this location +(though this more commonly gives a bracketing interval error), or if the emission radius is within +the event horizon. +""" +function find_offset_for_radius( + m::AbstractMetricParams, + u, + d::AbstractAccretionDisc, + rₑ, + θₒ; + zero_atol = 1e-7, + offset_max = 20.0, + kwargs..., +) + measure = gp -> begin + r = if gp.status == StatusCodes.IntersectedWithGeometry + gp.u2[2] * sin(gp.u2[3]) + else + 0.0 + end + rₑ - r + end + + f = r -> begin + gp = integrate_single_geodesic(m, u, d, r, θₒ; kwargs...) + measure(gp) + end + + r0 = Roots.find_zero(f, (0.0, offset_max); atol = zero_atol) + + gp0 = integrate_single_geodesic(m, u, d, r0, θₒ; kwargs...) + if !isapprox(measure(gp0), 0.0, atol = 10 * zero_atol) + return NaN, gp0 + end + r0, gp0 +end + +""" + impact_parameters_for_radius!( + αs::AbstractVector, + βs::AbstractVector, + m::AbstractMetricParams, + u::AbstractVector{T}, + d::AbstractAccretionDisc, + rₑ; + kwargs..., + ) + +Finds ``\\alpha`` and ``\\beta`` impact parameters which trace geodesics to a given emission radius +`rₑ` on the accretion disc. + +This function assigns the values in-place to `αs` and `βs`. +The keyword arguments are forwarded to [`find_offset_for_radius`](@ref). + +This function is threaded. +""" +function impact_parameters_for_radius!( + αs::AbstractVector, + βs::AbstractVector, + m::AbstractMetricParams, + u::AbstractVector{T}, + d::AbstractAccretionDisc, + rₑ; + kwargs..., +) where {T} + if size(α) != size(β) + throw(DimensionMismatch("α, β must have the same dimensions and size.")) + end + θs = range(0, 2π, length(α)) + @inbounds @threads for i in eachindex(θs) + θ = θs[i] + r, _ = find_offset_for_radius(m, u, d, rₑ, θ; kwargs...) + α[i] = r * cos(θ) + β[i] = r * sin(θ) + end + (α, β) +end + +""" + function impact_parameters_for_radius( + m::AbstractMetricParams, + u::AbstractVector{T}, + d::AbstractAccretionDisc, + radius; + N = 500, + kwargs..., + ) + +Pre-allocating version of [`impact_parameters_for_radius!`](@ref), which allocates +for `N` impact paramter pairs. Filters for `NaN` and returns `(αs, βs)`. +""" +function impact_parameters_for_radius( + m::AbstractMetricParams, + u::AbstractVector{T}, + d::AbstractAccretionDisc, + radius; + N = 500, + kwargs..., +) where {T} + α = zeros(T, N) + β = zeros(T, N) + impact_parameters_for_radius!(α, β, m, u, d, radius; kwargs...) + (filter(!isnan, α), filter(!isnan, β)) +end + +function jacobian_∂αβ_∂gr( + m, + u, + d, + α, + β, + max_time; + diff_order = 5, + redshift_pf = ConstPointFunctions.redshift, + kwargs..., +) + # map impact parameters to r, g + f = + ((α, β),) -> begin + v = map_impact_parameters(m, u, α, β) + sol = + tracegeodesics(m, u, v, d, (0.0, max_time); save_on = false, kwargs...) + gp = getgeodesicpoint(m, sol) + g = redshift_pf(m, gp, 2000.0) + # return r and g* + @SVector [gp.u2[2], g] + end + + # choice between FiniteDifferences and FiniteDiff is tricky + # since FiniteDiff is so much faster, but seems to yield really bad jacobians + # for this specific problem, so instead stenciling with a given order + cfdm = FiniteDifferences.central_fdm(diff_order, 1) + j = FiniteDifferences.jacobian(cfdm, f, @SVector([α, β])) |> first + abs(inv(det(j))) +end diff --git a/src/corona-to-disc/transfer-functions.jl b/src/transfer-functions/transfer-functions-2d.jl similarity index 100% rename from src/corona-to-disc/transfer-functions.jl rename to src/transfer-functions/transfer-functions-2d.jl diff --git a/test/smoke-tests/cunningham-transfer-functions.jl b/test/smoke-tests/cunningham-transfer-functions.jl index ca5c9e74..e3e79431 100644 --- a/test/smoke-tests/cunningham-transfer-functions.jl +++ b/test/smoke-tests/cunningham-transfer-functions.jl @@ -6,8 +6,8 @@ # test for different angles for (angle, expected) in zip( [3, 35, 74, 85], - # values last updated: 16/09/2022 - reduced number of points in test set - [49.09620827339265, 45.16844822761859, 27.107307569529434, 19.664521288378978], + # values last updated: 10/11/2022 - reimplemented algorithm + [24.953515365137005, 22.91076063485734, 13.795031677382394, 11.086660268049126], ) u = @SVector [0.0, 1000.0, deg2rad(angle), 0.0] ctf = cunningham_transfer_function( @@ -15,19 +15,16 @@ u, d, 4.0, - 2000.0; - finite_diff_order = 5, - num_points = 200, ) s = sum(ctf.f) - @test isapprox(s, expected, atol = 1e-2, rtol = 0.0) + @test isapprox(s, expected, atol = 1e-1, rtol = 0.0) end # different radii for (r, expected) in zip( [4.0, 7.0, 10.0, 15.0], - # values last updated: 16/09/2022 - reduced number of points in test set - [46.446844434569144, 51.08146707424876, 52.788771093295686, 53.78540573544522], + # values last updated: 10/11/2022 - reimplemented algorithm + [23.491668065237214, 25.997534896876644, 26.73859117429336, 27.186165046726007], ) u = @SVector [0.0, 1000.0, deg2rad(30), 0.0] ctf = cunningham_transfer_function( @@ -35,9 +32,6 @@ u, d, r, - 2000.0; - finite_diff_order = 5, - num_points = 200, ) s = sum(ctf.f) @test isapprox(s, expected, atol = 1e-2, rtol = 0.0) diff --git a/test/smoke-tests/disc-profiles.jl b/test/smoke-tests/disc-profiles.jl index 1963b3f5..70487aab 100644 --- a/test/smoke-tests/disc-profiles.jl +++ b/test/smoke-tests/disc-profiles.jl @@ -15,7 +15,8 @@ sampler = EvenSampler(domain = LowerHemisphere()), ) - intersected_simsols = filter(i -> i.prob.p.status == StatusCodes.IntersectedWithGeometry, simsols.u) + intersected_simsols = + filter(i -> i.prob.p.status == StatusCodes.IntersectedWithGeometry, simsols.u) sd_endpoints = map(sol -> getgeodesicpoint(m, sol), intersected_simsols) # test ensemble solution constructor diff --git a/test/smoke-tests/line-profiles.jl b/test/smoke-tests/line-profiles.jl index b575ec24..14fbe97e 100644 --- a/test/smoke-tests/line-profiles.jl +++ b/test/smoke-tests/line-profiles.jl @@ -6,17 +6,16 @@ @testset "cunningham" begin x = range(0.1, 1.2, 20) - bins, lp = lineprofile( - (re) -> re^(-3), + bins = range(0.0, 1.5, 200) + _, lp = lineprofile( + bins, + (r) -> r^(-3), m, u, d; - num_points = 40, - bins = x, - num_re = 40, - max_re = 50, - finite_diff_order = 5, + N = 40, + numrₑ = 30 ) - @test isapprox(0.20282590591902194, sum(lp), atol = 1e-4) + @test isapprox(1.0, sum(lp), atol = 1e-4) end end