From e6a614325adf3f6a91e4d1cd389615805e535208 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Sat, 9 Sep 2023 21:51:05 +0100 Subject: [PATCH 1/5] feat: add dot product check to visibility criteria Also compute the dot product between disc normal and velocity when trying to ascertain whether a point on the disc is visible or not. Caution against using just the projected radius, as this is biases towards visibility when the disc is steeply inclined. Prefer checking the Cartesian distance instead. --- src/corona/samplers.jl | 8 ++++ src/geometry/discs.jl | 2 +- src/geometry/discs/thick-disc.jl | 22 +++++++++ src/tracing/precision-solvers.jl | 46 ++++++++++++------ .../cunningham-transfer-functions.jl | 48 +++---------------- 5 files changed, 69 insertions(+), 57 deletions(-) diff --git a/src/corona/samplers.jl b/src/corona/samplers.jl index 37cd31dd..5663a89b 100644 --- a/src/corona/samplers.jl +++ b/src/corona/samplers.jl @@ -63,6 +63,14 @@ function _cart_to_spher_jacobian(θ, ϕ) ] end +function _spher_to_cart_jacobian(θ, ϕ, r) + @SMatrix [ + sin(θ)*cos(ϕ) r*cos(θ)*cos(ϕ) -r*sin(θ)*sin(ϕ) + sin(θ)*sin(ϕ) r*cos(θ)*sin(ϕ) r*sin(θ)*cos(ϕ) + cos(θ) -r*sin(θ) 0 + ] +end + function _cart_local_direction(θ, ϕ) @SVector [sin(θ) * cos(ϕ), sin(θ) * sin(ϕ), cos(θ)] end diff --git a/src/geometry/discs.jl b/src/geometry/discs.jl index 93c30755..4c852131 100644 --- a/src/geometry/discs.jl +++ b/src/geometry/discs.jl @@ -45,7 +45,7 @@ end cross_section(d::AbstractThickAccretionDisc, x4) = error("Not implemented for $(typeof(d)).") r_cross_section(d::AbstractThickAccretionDisc, r::Number) = - cross_section(d, SVector(0.0, r, π / 2, 0.0)) + cross_section(d, SVector(0, r, π / 2, 0)) struct EllipticalDisc{T} <: AbstractAccretionDisc{T} inner_radius::T diff --git a/src/geometry/discs/thick-disc.jl b/src/geometry/discs/thick-disc.jl index afbd213f..247338ec 100644 --- a/src/geometry/discs/thick-disc.jl +++ b/src/geometry/discs/thick-disc.jl @@ -61,5 +61,27 @@ function distance_to_disc(d::AbstractThickAccretionDisc, x4; gtol) abs(z) - height - (gtol * x4[2]) end +function _cartesian_tangent_vector(ρ, d::AbstractThickAccretionDisc) + function _parameterization(r) + xz_parameterize(d, r) + end + ∇f = ForwardDiff.derivative(_parameterization, ρ) + v = SVector(∇f[1], 0, ∇f[2]) + v ./ √(v[1]^2 + v[2]^2 + v[3]^2) +end + +function _cartesian_surface_normal(ρ, d::AbstractThickAccretionDisc) + ∇f = _cartesian_tangent_vector(ρ, d) + # rotate 90 degrees in about ϕ̂ + SVector(-∇f[3], ∇f[2], ∇f[1]) +end + +function _rotate_cartesian_about_z(n, ϕ) + SVector(n[1] * cos(ϕ), n[1] * sin(ϕ), n[3]) +end + +function _cartesian_surface_normal(ρ, ϕ, d::AbstractThickAccretionDisc) + _rotate_cartesian_about_z(_cartesian_surface_normal(ρ, d), ϕ) +end export ThickDisc diff --git a/src/tracing/precision-solvers.jl b/src/tracing/precision-solvers.jl index 5a0530a4..3e3797e3 100644 --- a/src/tracing/precision-solvers.jl +++ b/src/tracing/precision-solvers.jl @@ -182,7 +182,6 @@ function impact_parameters_for_radius_obscured( x::AbstractVector{T}, d::AbstractAccretionDisc, radius; - visible_tolerance = 1e-3, N = 500, β₀ = zero(T), α₀ = zero(T), @@ -192,26 +191,43 @@ function impact_parameters_for_radius_obscured( β = zeros(T, N) impact_parameters_for_radius!(α, β, m, x, d, radius; β₀ = β₀, α₀ = α₀, kwargs...) - I = map( - i -> - !_trace_is_visible( - m, - x, - α[i], - β[i], - d, - radius; - visible_tolerance = visible_tolerance, - kwargs..., - ), - eachindex(α), + gps = tracegeodesics( + m, + x, + i -> map_impact_parameters(m, x, α[i], β[i]), + datumplane(d, radius), + 2x[2]; + ensemble = EnsembleEndpointThreads(), + save_on = false, + trajectories = length(α), ) - + n = _cartesian_surface_normal(radius, d) + I = [!_is_visible(m, d, gp, n) for gp in gps] α[I] .= NaN β[I] .= NaN α, β end +""" +If dot product between surface normal and the final velocity vector is positive, +the point is visible. +""" +function _is_visible(m::AbstractMetric, d, gp::AbstractGeodesicPoint, n::SVector{3}) + gp_new = + tracegeodesics(m, gp.x_init, gp.v_init, d, gp.λ_max; save_on = false) |> + unpack_solution + # geodesics sufficiently close, test the angle + v = SVector(gp_new.v[2], gp_new.v[3], gp_new.v[4]) ./ gp_new.v[1] + v_geodesic = _spher_to_cart_jacobian(gp_new.x[3], gp_new.x[4], gp_new.x[2]) * v + n_rot = _rotate_cartesian_about_z(n, gp_new.x[4]) + + X1 = to_cartesian(gp.x) + X2 = to_cartesian(gp_new.x) + dist = sum(i -> i^2, X1 .- X2) + + return _fast_dot(v_geodesic, n_rot) < 0 && isapprox(dist, 0, atol = 1e-10) +end + function jacobian_∂αβ_∂gr( m, x, diff --git a/src/transfer-functions/cunningham-transfer-functions.jl b/src/transfer-functions/cunningham-transfer-functions.jl index 1520418a..f9e87eec 100644 --- a/src/transfer-functions/cunningham-transfer-functions.jl +++ b/src/transfer-functions/cunningham-transfer-functions.jl @@ -187,15 +187,15 @@ function _rear_workhorse_with_impact_parameters( redshift_pf = redshift_pf, tracer_kwargs..., ) - (g, J, gp.x[1], α, β) + (g, J, gp, α, β) end (_workhorse, tracer_kwargs) end function _rear_workhorse(m::AbstractMetric, x, d::AbstractAccretionDisc, rₑ; kwargs...) workhorse, _ = _rear_workhorse_with_impact_parameters(m, x, d, rₑ; kwargs...) function _disc_workhorse(θ::T)::NTuple{3,T} where {T} - g, J, t, _, _ = workhorse(θ) - g, J, t + g, J, gp, _, _ = workhorse(θ) + g, J, gp.x[1] end end @@ -205,7 +205,6 @@ function _rear_workhorse( d::AbstractThickAccretionDisc, rₑ; max_time = 2 * x[2], - visible_tolerance = 1e-3, kwargs..., ) plane = datumplane(d, rₑ) @@ -217,47 +216,14 @@ function _rear_workhorse( max_time = max_time, kwargs..., ) + n = _cartesian_surface_normal(rₑ, d) function _thick_workhorse(θ::T)::NTuple{4,T} where {T} - g, J, t, α, β = datum_workhorse(θ) - # check if the location is visible or not - is_visible = _trace_is_visible( - m, - x, - α, - β, - d, - rₑ; - visible_tolerance = visible_tolerance, - tracer_kwargs..., - ) - (g, J, t, is_visible) + g, J, gp, _, _ = datum_workhorse(θ) + is_visible = _is_visible(m, d, gp, n) + (g, J, gp.x[1], is_visible) end end -function _trace_is_visible( - m, - x, - α, - β, - d, - rₑ; - max_time = 2x[2], - visible_tolerance = 1e-3, - kwargs..., -) - sol = tracegeodesics( - m, - x, - map_impact_parameters(m, x, α, β), - d, - max_time; - save_on = false, - kwargs..., - ) - gp = unpack_solution(sol) - abs(rₑ - _equatorial_project(gp.x)) < visible_tolerance -end - function _cunningham_transfer_function!( data::_TransferDataAccumulator, workhorse, From 76808a7b88c1f54b8486470efcd6709e43dbd296 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Sat, 9 Sep 2023 21:53:52 +0100 Subject: [PATCH 2/5] fix: missing aux in unpack_solution_full --- src/solution-processing.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/solution-processing.jl b/src/solution-processing.jl index 2210f986..1127ac70 100644 --- a/src/solution-processing.jl +++ b/src/solution-processing.jl @@ -114,6 +114,7 @@ function unpack_solution_full( ui = SVector{4,T}(us[i][1:4]) vi = SVector{4,T}(us[i][5:8]) ti = ts[i] + aux = unpack_auxiliary(us[i]) GeodesicPoint( get_status_code(sol.prob.p), t_start, @@ -122,6 +123,7 @@ function unpack_solution_full( ui, v_start, vi, + aux, ) end end From f66fd673e186f4efdf12611c88cf5da6c3dde7a6 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Sat, 9 Sep 2023 21:54:17 +0100 Subject: [PATCH 3/5] fix: make plotting aware of status of geodesic when deciding the midpoint for interpolation --- src/plotting-recipes.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/plotting-recipes.jl b/src/plotting-recipes.jl index 34999eba..1472f3c5 100644 --- a/src/plotting-recipes.jl +++ b/src/plotting-recipes.jl @@ -1,7 +1,14 @@ using RecipesBase -function _extract_path(sol, n_points; projection = :none, t_span = 100.0) - mid_i = max(1, length(sol.u) ÷ 2) +function _extract_path(sol, n_points; projection = :none, t_span = 30.0) + status = Gradus.get_status_code(sol.prob.p) + mid_i = + if (status == StatusCodes.IntersectedWithGeometry) || + (status == StatusCodes.WithinInnerBoundary) + lastindex(sol.t) + else + max(1, length(sol.u) ÷ 2) + end start_t = max(sol.t[mid_i] - t_span, sol.t[1]) end_t = min(sol.t[mid_i] + t_span, sol.t[end]) @@ -31,13 +38,13 @@ end zlims --> _range aspect_ratio --> 1 - itr = if !(typeof(sol) <: SciMLBase.EnsembleSolution) - (; u = (sol,)) + itr = if typeof(sol) <: SciMLBase.EnsembleSolution + sol.u else sol end - for s in itr.u + for s in itr path = _extract_path(s, n_points, projection = :none, t_span = t_span) @series begin path From 0feb2a9128fe8e9a9f3ca207c0e265a175164404 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Sat, 9 Sep 2023 22:43:01 +0100 Subject: [PATCH 4/5] fix: the most annoying bug so far related to jacobians and thick discs --- src/transfer-functions/cunningham-transfer-functions.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/transfer-functions/cunningham-transfer-functions.jl b/src/transfer-functions/cunningham-transfer-functions.jl index f9e87eec..837423fc 100644 --- a/src/transfer-functions/cunningham-transfer-functions.jl +++ b/src/transfer-functions/cunningham-transfer-functions.jl @@ -147,11 +147,12 @@ function _rear_workhorse_with_impact_parameters( d::AbstractAccretionDisc, rₑ; max_time = 2 * x[2], - redshift_pf = ConstPointFunctions.redshift(m, x), offset_max = 0.4rₑ + 10, zero_atol = 1e-7, β₀ = 0, α₀ = 0, + redshift_pf = ConstPointFunctions.redshift(m, x), + jacobian_disc = d, tracer_kwargs..., ) function _workhorse(θ) @@ -180,7 +181,7 @@ function _rear_workhorse_with_impact_parameters( J = jacobian_∂αβ_∂gr( m, x, - d, + jacobian_disc, α, β, max_time; @@ -214,6 +215,7 @@ function _rear_workhorse( plane, rₑ; max_time = max_time, + jacobian_disc = d, kwargs..., ) n = _cartesian_surface_normal(rₑ, d) From 156458ec57a4309e43e0064fd2d4b14d43f63afe Mon Sep 17 00:00:00 2001 From: fjebaker Date: Sat, 9 Sep 2023 22:56:12 +0100 Subject: [PATCH 5/5] test: update thick disc test value --- test/transfer-functions/test-thick-disc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/transfer-functions/test-thick-disc.jl b/test/transfer-functions/test-thick-disc.jl index 077aa0cd..4f60c330 100644 --- a/test/transfer-functions/test-thick-disc.jl +++ b/test/transfer-functions/test-thick-disc.jl @@ -8,4 +8,4 @@ d = ShakuraSunyaev(m) tf = cunningham_transfer_function(m, x, d, 3.0; β₀ = 1.0) total = sum(filter(!isnan, tf.f)) -@test total ≈ 13.419539069804333 atol = 1e-4 +@test total ≈ 12.148257594740231 atol = 1e-4