From 803cfa7ec2d9623c1a7c14164e96944f36ca7310 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Mon, 24 Jul 2023 14:59:52 +0100 Subject: [PATCH 1/9] feat!: coronal emissivity API rework The `emissivity_profile` function has now been superseded by a more descriptive `tracecorona`, which returns a `CoronaGeodesics` struct that may subsequently be passed to a profile struct. The `emissivity_profile` function now instead always assumes radial symmetry and returns a `RadialDiscProfile` struct, with the intention that this function may be optimized better for point source coronas. Tests have also been updated to accommodate these API changes. --- src/corona/corona-models.jl | 19 +----- src/corona/disc-profiles.jl | 15 ++++- src/corona/emissivity.jl | 65 ++++++++++++------- src/corona/flux-calculations.jl | 1 - src/corona/models/lamp-post.jl | 13 ++++ .../transfer-functions-2d.jl | 2 +- test/transfer-functions/test-2d.jl | 3 +- 7 files changed, 73 insertions(+), 45 deletions(-) create mode 100644 src/corona/models/lamp-post.jl diff --git a/src/corona/corona-models.jl b/src/corona/corona-models.jl index afc1e34a..cb0d64d4 100644 --- a/src/corona/corona-models.jl +++ b/src/corona/corona-models.jl @@ -80,22 +80,6 @@ function sample_local_velocity( sky_angles_to_velocity(m, x, v, θ, ϕ) end -# model implementations - -@with_kw struct LampPostModel{T} <: AbstractCoronaModel{T} - @deftype T - h = 5.0 - θ = 0.01 - ϕ = 0.0 -end - -function sample_position_velocity(m::AbstractMetric, model::LampPostModel{T}) where {T} - x = SVector{4,T}(0, model.h, model.θ, model.ϕ) - gcomp = metric_components(m, SVector(x[2], x[3])) - v = inv(√(-gcomp[1])) * SVector{4,T}(1, 0, 0, 0) - x, v -end - # bootstrap tracing function for convenience function tracegeodesics( m::AbstractMetric, @@ -109,4 +93,7 @@ function tracegeodesics( tracegeodesics(m, xs, vs, args...; kwargs...) end + +include("models/lamp-post.jl") + export LampPostModel diff --git a/src/corona/disc-profiles.jl b/src/corona/disc-profiles.jl index dfc4c52f..c2c32230 100644 --- a/src/corona/disc-profiles.jl +++ b/src/corona/disc-profiles.jl @@ -46,6 +46,13 @@ function RadialDiscProfile( RadialDiscProfile(m, radii, times, gs, intensities; kwargs...) end +""" +The relativistic correction is calculated via + +```math +\\tilde{A} = A \\sqrt{g_{\\mu,\\nu}(x)} +``` +""" function RadialDiscProfile( m::AbstractMetric, radii::AbstractArray{T}, @@ -80,11 +87,15 @@ function RadialDiscProfile( R = bins[i] r = i == 1 ? 0 : bins[i-1] + x = SVector(R, π / 2) + gcomp = metric_components(m, x) + # 2π comes from integrating over dϕ dr = R - r - A = 2π * dr + A = 2π * dr * √(gcomp[2] * gcomp[4]) + # I now stores emissivity - I[i] = source_to_disc_emissivity(m, I[i], A, SVector(R, π / 2), eratios(R)) + I[i] = source_to_disc_emissivity(m, I[i], A, x, eratios(R)) end # create interpolations diff --git a/src/corona/emissivity.jl b/src/corona/emissivity.jl index b8d7f360..88627b13 100644 --- a/src/corona/emissivity.jl +++ b/src/corona/emissivity.jl @@ -1,31 +1,23 @@ """ source_to_disc_emissivity(m, 𝒩, A, x, g) -Compute the emissivity (in arbitrary units) in the disc element with area `A`, photon +Compute the emissivity (in arbitrary units) in the disc element with proper area `A`, photon count `𝒩`, central position `x`, and redshift `g`. Evaluates ```math -\\varepsilon = \\frac{\\mathscr{N}}{\\tilde{A} g^2}, -``` - -where ``\\tilde{A}`` is the relativistically corrected area of `A`. The relativistic correction -is calculated via - -```math -\\tilde{A} = A \\sqrt{g_{\\mu,\\nu}(x)} +\\varepsilon = \\frac{\\mathscr{N}}{A g^2}. ``` """ -function source_to_disc_emissivity(m::AbstractStaticAxisSymmetric, 𝒩, A, x, g) - gcomp = metric_components(m, x) +function source_to_disc_emissivity(m::AbstractStaticAxisSymmetric, 𝒩, A, x, g; Γ = 2) v = CircularOrbits.fourvelocity(m, x) - # account for relativistic effects in area + # account for relativistic effects in area due to lorentz shift γ = lorentz_factor(m, SVector(0, x[1], x[2], 0), v) - A_corrected = A * √(gcomp[2] * gcomp[4]) # divide by area to get number density - 𝒩 / (g^2 * A_corrected * γ) + 𝒩 / (g^Γ * A * γ) end -struct CoronalEmissivity{T,M,G,C,P,V} + +struct CoronaGeodesics{T,M,G,C,P,V} trace::T metric::M geometry::G @@ -34,7 +26,7 @@ struct CoronalEmissivity{T,M,G,C,P,V} source_velocity::V end -function emissivity_profile( +function tracecorona( m::AbstractMetric, g::AbstractAccretionGeometry, model::AbstractCoronaModel; @@ -42,9 +34,9 @@ function emissivity_profile( n_samples = 1024, sampler = EvenSampler(domain = BothHemispheres(), generator = RandomGenerator()), trace = TraceGeodesic(), + callback = domain_upper_hemisphere(), kwargs..., ) - xs, vs, source_vels = sample_position_direction_velocity(m, model, sampler, n_samples) gps = tracegeodesics( m, @@ -55,14 +47,14 @@ function emissivity_profile( trace = trace, save_on = false, ensemble = EnsembleEndpointThreads(), - callback = domain_upper_hemisphere(), + callback = callback, kwargs..., ) mask = [i.status == StatusCodes.IntersectedWithGeometry for i in gps] - CoronalEmissivity(trace, m, g, model, gps[mask], source_vels[mask]) + CoronaGeodesics(trace, m, g, model, gps[mask], source_vels[mask]) end -function RadialDiscProfile(ce::CoronalEmissivity; kwargs...) +function RadialDiscProfile(ce::CoronaGeodesics; kwargs...) J = sortperm(ce.geodesic_points; by = i -> i.x[2]) @views RadialDiscProfile( ce.metric, @@ -72,7 +64,7 @@ function RadialDiscProfile(ce::CoronalEmissivity; kwargs...) ) end -function RadialDiscProfile(ce::CoronalEmissivity{<:TraceRadiativeTransfer}; kwargs...) +function RadialDiscProfile(ce::CoronaGeodesics{<:TraceRadiativeTransfer}; kwargs...) J = sortperm(ce.geodesic_points; by = i -> i.x[2]) @views RadialDiscProfile( ce.metric, @@ -83,7 +75,7 @@ function RadialDiscProfile(ce::CoronalEmissivity{<:TraceRadiativeTransfer}; kwar ) end -function RadialDiscProfile(ε, ce::CoronalEmissivity; kwargs...) +function RadialDiscProfile(ε, ce::CoronaGeodesics; kwargs...) J = sortperm(ce.geodesic_points; by = i -> i.x[2]) radii = @views [i.x[2] for i in ce.geodesic_points[J]] times = @views [i.x[1] for i in ce.geodesic_points[J]] @@ -101,4 +93,31 @@ function RadialDiscProfile(ε, ce::CoronalEmissivity; kwargs...) RadialDiscProfile(_emissivity_wrapper, _delay_wrapper) end -export emissivity_profile +""" +function emissivity_profile( + m::AbstractMetric, + d::AbstractAccretionGeometry, + model::AbstractCoronaModel; + N = 100, + kwargs..., +end + +Calculates a [`RadialDiscProfile`](@ref). + +This function assumes axis symmetry, and therefore always interpolates the emissivity +as a function of the radial coordinate on the disc. If non-symmetric profiles are +desired, consider using [`VoronoiDiscProfile`](@ref). +""" +function emissivity_profile( + m::AbstractMetric, + d::AbstractAccretionGeometry, + model::AbstractCoronaModel; + grid = GeometricGrid(), + N = 100, + kwargs..., +) + RadialDiscProfile(tracecorona(m, d, model; kwargs...); grid = grid, N = N) +end + + +export emissivity_profile, tracecorona diff --git a/src/corona/flux-calculations.jl b/src/corona/flux-calculations.jl index 948b8375..77548864 100644 --- a/src/corona/flux-calculations.jl +++ b/src/corona/flux-calculations.jl @@ -105,7 +105,6 @@ function energy_ratio(m, gp, v_src) e_src / e_disc end - function flux_source_to_disc( m::AbstractMetric, model::AbstractCoronaModel, diff --git a/src/corona/models/lamp-post.jl b/src/corona/models/lamp-post.jl new file mode 100644 index 00000000..82b40195 --- /dev/null +++ b/src/corona/models/lamp-post.jl @@ -0,0 +1,13 @@ +@with_kw struct LampPostModel{T} <: AbstractCoronaModel{T} + @deftype T + h = 5.0 + θ = 0.01 + ϕ = 0.0 +end + +function sample_position_velocity(m::AbstractMetric, model::LampPostModel{T}) where {T} + x = SVector{4,T}(0, model.h, model.θ, model.ϕ) + gcomp = metric_components(m, SVector(x[2], x[3])) + v = inv(√(-gcomp[1])) * SVector{4,T}(1, 0, 0, 0) + x, v +end diff --git a/src/transfer-functions/transfer-functions-2d.jl b/src/transfer-functions/transfer-functions-2d.jl index 6615afea..7d0d7c6b 100644 --- a/src/transfer-functions/transfer-functions-2d.jl +++ b/src/transfer-functions/transfer-functions-2d.jl @@ -144,7 +144,7 @@ function lagtransfer( verbose, ) - ce = emissivity_profile( + ce = tracecorona( m, d, model; diff --git a/test/transfer-functions/test-2d.jl b/test/transfer-functions/test-2d.jl index 4f041f46..868bf186 100644 --- a/test/transfer-functions/test-2d.jl +++ b/test/transfer-functions/test-2d.jl @@ -33,7 +33,7 @@ fluxsum = sum(filter(!isnan, f)) # test semi-analytic method -ep = @time Gradus.emissivity_profile( +prof = @time Gradus.emissivity_profile( m, d, model; @@ -41,7 +41,6 @@ ep = @time Gradus.emissivity_profile( callback = domain_upper_hemisphere(), sampler = EvenSampler(domain = BothHemispheres(), generator = GoldenSpiralGenerator()), ) -prof = RadialDiscProfile(ep) radii = Gradus.Grids._inverse_grid(Gradus.isco(m), 100.0, 5) d = GeometricThinDisc(0.0, 500.0, π / 2) From 29bf67dbda09349fd13d35e5b41c4932e3534dad Mon Sep 17 00:00:00 2001 From: fjebaker Date: Mon, 24 Jul 2023 15:44:03 +0100 Subject: [PATCH 2/9] feat: optional normalisation for emissivity profile plots --- src/plotting-recipes.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/plotting-recipes.jl b/src/plotting-recipes.jl index f33817db..cef56c61 100644 --- a/src/plotting-recipes.jl +++ b/src/plotting-recipes.jl @@ -127,7 +127,7 @@ end (a, b) end -@recipe function f(p::RadialDiscProfile) +@recipe function f(p::RadialDiscProfile; normalize = identity) legend --> false xlabel --> "r (rg)" ylabel --> "ε (arb.)" @@ -135,8 +135,7 @@ end yscale --> :log10 y = p.f.ε.u[2:end-1] - y = y ./ sum(y) - y = y ./ minimum(y) + @. y = normalize(y) p.f.ε.t[2:end-1], y end From 6411a6c06245bcd4de82ea1966705d3eda988569 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Mon, 24 Jul 2023 15:46:22 +0100 Subject: [PATCH 3/9] feat: emissivity for point-like corona with symmetry --- src/corona/disc-profiles.jl | 28 ++++++++----- src/corona/emissivity.jl | 75 ++++++++++++++++++++++++++++++++++ src/corona/models/lamp-post.jl | 8 ++++ 3 files changed, 101 insertions(+), 10 deletions(-) diff --git a/src/corona/disc-profiles.jl b/src/corona/disc-profiles.jl index c2c32230..d520861c 100644 --- a/src/corona/disc-profiles.jl +++ b/src/corona/disc-profiles.jl @@ -46,6 +46,22 @@ function RadialDiscProfile( RadialDiscProfile(m, radii, times, gs, intensities; kwargs...) end +function RadialDiscProfile(rs::AbstractArray, ts::AbstractArray, εs::AbstractArray) + # create interpolations + t = DataInterpolations.LinearInterpolation(ts, rs) + ε = DataInterpolations.LinearInterpolation(εs, rs) + # wrap geodesic point wrappers + RadialDiscProfile(gp -> ε(gp.x[2]), gp -> t(gp.x[2]) + gp.x[1]) +end + +function RadialDiscProfile(rdp::RadialDiscProfile) + Base.depwarn( + "This function is deprecated. Note that `emissivity_profile` now returns a `RadialDiscProfile`.", + :RadialDiscProfile, + ) + rdp +end + """ The relativistic correction is calculated via @@ -87,23 +103,15 @@ function RadialDiscProfile( R = bins[i] r = i == 1 ? 0 : bins[i-1] - x = SVector(R, π / 2) - gcomp = metric_components(m, x) - # 2π comes from integrating over dϕ dr = R - r - A = 2π * dr * √(gcomp[2] * gcomp[4]) + A = dr * _proper_area(m, SVector(0, R, π / 2, 0)) # I now stores emissivity I[i] = source_to_disc_emissivity(m, I[i], A, x, eratios(R)) end - # create interpolations - t = DataInterpolations.LinearInterpolation(ts, bins) - ε = DataInterpolations.LinearInterpolation(I, bins) - - # wrap geodesic point wrappers - RadialDiscProfile(gp -> ε(gp.x[2]), gp -> t(gp.x[2]) + gp.x[1]) + RadialDiscProfile(bins, ts, I) end emitted_flux(profile::RadialDiscProfile, gps) = map(profile.f, gps) diff --git a/src/corona/emissivity.jl b/src/corona/emissivity.jl index 88627b13..a72ef927 100644 --- a/src/corona/emissivity.jl +++ b/src/corona/emissivity.jl @@ -16,6 +16,9 @@ function source_to_disc_emissivity(m::AbstractStaticAxisSymmetric, 𝒩, A, x, g 𝒩 / (g^Γ * A * γ) end +function source_to_disc_emissivity(δ, g, A, γ; Γ = 2) + sin(δ) / (g^Γ * A * γ) +end struct CoronaGeodesics{T,M,G,C,P,V} trace::T @@ -119,5 +122,77 @@ function emissivity_profile( RadialDiscProfile(tracecorona(m, d, model; kwargs...); grid = grid, N = N) end +function _proper_area(m, x::SVector{4}) + gcomp = Gradus.metric_components(m, SVector(x[2], x[3])) + det_g = √(gcomp[2] * gcomp[4]) + 2π * det_g +end + +function polar_angle_to_velfunc(m::AbstractMetric, x, v, δs) + function _polar_angle_velfunc(i) + sky_angles_to_velocity(m, x, v, δs[i], 0.0) + end +end + +function _point_source_symmetric_emissivity_profile( + m::AbstractMetric, + d::AbstractAccretionGeometry, + model::AbstractCoronaModel; + N = 100, + λ_max = 10_000.0, + callback = domain_upper_hemisphere(), + kwargs..., +) + δs = deg2rad.(range(0.1, 179.9, N)) + # we assume a point source + x, v = Gradus.sample_position_velocity(m, model) + velfunc = polar_angle_to_velfunc(m, x, v, δs) + gps = tracegeodesics( + m, + x, + velfunc, + d, + λ_max; + save_on = false, + ensemble = EnsembleEndpointThreads(), + callback = callback, + trajectories = length(δs), + kwargs..., + ) + + # filter only those that intersected, and sort radially + I = [i.status == StatusCodes.IntersectedWithGeometry for i in gps] + points = gps[I] + δs = δs[I] + J = sortperm(points, by = i -> i.x[2]) + points = points[J] + δs = δs[J] + + r, ε = _point_source_emissivity(m, d, v, δs, points) + t = [i.x[1] for i in @views(points[1:end-1])] + + RadialDiscProfile(r, t, ε) +end + +function _point_source_emissivity( + m::AbstractMetric, + d::AbstractAccretionGeometry, + source_velocity, + δs, + points, +) + _points = @views points[1:end-1] + # get the time and radial position + r = [i.x[2] for i in points] + # create a view + gs = Gradus.energy_ratio.(m, _points, (source_velocity,)) + A = [_proper_area(m, i.x) for i in _points] + γ = [Gradus.lorentz_factor(m, d, i.x) for i in _points] + Δr = diff(r) + @. A = A * Δr + r = r[1:end-1] + ε = source_to_disc_emissivity.(@views(δs[1:end-1]), gs, A, γ) + r, ε +end export emissivity_profile, tracecorona diff --git a/src/corona/models/lamp-post.jl b/src/corona/models/lamp-post.jl index 82b40195..f0fd8dd2 100644 --- a/src/corona/models/lamp-post.jl +++ b/src/corona/models/lamp-post.jl @@ -11,3 +11,11 @@ function sample_position_velocity(m::AbstractMetric, model::LampPostModel{T}) wh v = inv(√(-gcomp[1])) * SVector{4,T}(1, 0, 0, 0) x, v end + +# can exploit point source symmetry for certain disc models +emissivity_profile( + m::AbstractMetric, + d::AbstractAccretionDisc, + model::LampPostModel; + kwargs..., +) = _point_source_symmetric_emissivity_profile(m, d, model; kwargs...) From 45296509da9c36be026b906da07fd565af78aa17 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Mon, 24 Jul 2023 16:00:01 +0100 Subject: [PATCH 4/9] fix: dispatch angular method based on presence of sampler kwarg --- src/corona/disc-profiles.jl | 4 ++-- src/corona/emissivity.jl | 24 ++++++++++++++++++------ src/corona/models/lamp-post.jl | 1 + 3 files changed, 21 insertions(+), 8 deletions(-) diff --git a/src/corona/disc-profiles.jl b/src/corona/disc-profiles.jl index d520861c..b756f876 100644 --- a/src/corona/disc-profiles.jl +++ b/src/corona/disc-profiles.jl @@ -105,8 +105,8 @@ function RadialDiscProfile( # 2π comes from integrating over dϕ dr = R - r - A = dr * _proper_area(m, SVector(0, R, π / 2, 0)) - + x = SVector(0, R, π / 2, 0) + A = dr * _proper_area(m, x) # I now stores emissivity I[i] = source_to_disc_emissivity(m, I[i], A, x, eratios(R)) end diff --git a/src/corona/emissivity.jl b/src/corona/emissivity.jl index a72ef927..c718ddd1 100644 --- a/src/corona/emissivity.jl +++ b/src/corona/emissivity.jl @@ -9,9 +9,9 @@ count `𝒩`, central position `x`, and redshift `g`. Evaluates ``` """ function source_to_disc_emissivity(m::AbstractStaticAxisSymmetric, 𝒩, A, x, g; Γ = 2) - v = CircularOrbits.fourvelocity(m, x) + v = CircularOrbits.fourvelocity(m, SVector(x[2], x[3])) # account for relativistic effects in area due to lorentz shift - γ = lorentz_factor(m, SVector(0, x[1], x[2], 0), v) + γ = lorentz_factor(m, x, v) # divide by area to get number density 𝒩 / (g^Γ * A * γ) end @@ -111,7 +111,15 @@ This function assumes axis symmetry, and therefore always interpolates the emiss as a function of the radial coordinate on the disc. If non-symmetric profiles are desired, consider using [`VoronoiDiscProfile`](@ref). """ +emissivity_profile( + m::AbstractMetric, + d::AbstractAccretionGeometry, + model::AbstractCoronaModel; + sampler = nothing, + kwargs..., +) = emissivity_profile(sampler, m, d, model; kwargs...) function emissivity_profile( + sampler::AbstractDirectionSampler, m::AbstractMetric, d::AbstractAccretionGeometry, model::AbstractCoronaModel; @@ -119,7 +127,11 @@ function emissivity_profile( N = 100, kwargs..., ) - RadialDiscProfile(tracecorona(m, d, model; kwargs...); grid = grid, N = N) + RadialDiscProfile( + tracecorona(m, d, model; sampler = sampler, kwargs...); + grid = grid, + N = N, + ) end function _proper_area(m, x::SVector{4}) @@ -138,14 +150,14 @@ function _point_source_symmetric_emissivity_profile( m::AbstractMetric, d::AbstractAccretionGeometry, model::AbstractCoronaModel; - N = 100, + n_samples = 100, λ_max = 10_000.0, callback = domain_upper_hemisphere(), kwargs..., ) - δs = deg2rad.(range(0.1, 179.9, N)) + δs = deg2rad.(range(0.1, 179.9, n_samples)) # we assume a point source - x, v = Gradus.sample_position_velocity(m, model) + x, v = sample_position_velocity(m, model) velfunc = polar_angle_to_velfunc(m, x, v, δs) gps = tracegeodesics( m, diff --git a/src/corona/models/lamp-post.jl b/src/corona/models/lamp-post.jl index f0fd8dd2..c26b6546 100644 --- a/src/corona/models/lamp-post.jl +++ b/src/corona/models/lamp-post.jl @@ -14,6 +14,7 @@ end # can exploit point source symmetry for certain disc models emissivity_profile( + ::Nothing, m::AbstractMetric, d::AbstractAccretionDisc, model::LampPostModel; From 6bacf36df2cf24e49ca3b2298901a2690ee405f3 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Mon, 24 Jul 2023 16:00:10 +0100 Subject: [PATCH 5/9] test: fix for api changes --- test/transfer-functions/test-2d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/transfer-functions/test-2d.jl b/test/transfer-functions/test-2d.jl index 868bf186..1133d518 100644 --- a/test/transfer-functions/test-2d.jl +++ b/test/transfer-functions/test-2d.jl @@ -33,12 +33,11 @@ fluxsum = sum(filter(!isnan, f)) # test semi-analytic method -prof = @time Gradus.emissivity_profile( +prof = emissivity_profile( m, d, model; n_samples = 5_000, - callback = domain_upper_hemisphere(), sampler = EvenSampler(domain = BothHemispheres(), generator = GoldenSpiralGenerator()), ) From 5c10ab54c9d3541e47ac58229acadbfd88793cd8 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Tue, 25 Jul 2023 08:40:59 +0100 Subject: [PATCH 6/9] chore: refactor code and disc profile structures Seperated the disc profile structures into their own logical units, so that all relevant code to a particular structure (e.g. Voronoi) is in the same place. Added a new dispatch to `emissivity_profile` that raises a descriptive error if an attempt to calculate a profile without a sampler and valid implementation is invoked. Rename `sky-geometry.jl` to `samplers.jl`, as it better illustrates what that file does. --- src/Gradus.jl | 2 +- src/corona/corona-models.jl | 36 ++ src/corona/disc-profiles.jl | 354 +------------------- src/corona/emissivity.jl | 95 +----- src/corona/profiles/radial.jl | 140 ++++++++ src/corona/profiles/voronoi.jl | 248 ++++++++++++++ src/corona/{sky-geometry.jl => samplers.jl} | 0 7 files changed, 445 insertions(+), 430 deletions(-) create mode 100644 src/corona/profiles/radial.jl create mode 100644 src/corona/profiles/voronoi.jl rename src/corona/{sky-geometry.jl => samplers.jl} (100%) diff --git a/src/Gradus.jl b/src/Gradus.jl index 8d852e49..057d7945 100644 --- a/src/Gradus.jl +++ b/src/Gradus.jl @@ -218,7 +218,7 @@ include("transfer-functions/types.jl") include("transfer-functions/cunningham-transfer-functions.jl") include("transfer-functions/integration.jl") -include("corona/sky-geometry.jl") +include("corona/samplers.jl") include("corona/corona-models.jl") include("corona/disc-profiles.jl") # needs the types from disc profiles so defer include diff --git a/src/corona/corona-models.jl b/src/corona/corona-models.jl index cb0d64d4..88fe6fd6 100644 --- a/src/corona/corona-models.jl +++ b/src/corona/corona-models.jl @@ -93,6 +93,42 @@ function tracegeodesics( tracegeodesics(m, xs, vs, args...; kwargs...) end +struct CoronaGeodesics{T,M,G,C,P,V} + trace::T + metric::M + geometry::G + model::C + geodesic_points::P + source_velocity::V +end + +function tracecorona( + m::AbstractMetric, + g::AbstractAccretionGeometry, + model::AbstractCoronaModel; + λ_max = 10_000, + n_samples = 1024, + sampler = EvenSampler(domain = BothHemispheres(), generator = RandomGenerator()), + trace = TraceGeodesic(), + callback = domain_upper_hemisphere(), + kwargs..., +) + xs, vs, source_vels = sample_position_direction_velocity(m, model, sampler, n_samples) + gps = tracegeodesics( + m, + xs, + vs, + g, + λ_max; + trace = trace, + save_on = false, + ensemble = EnsembleEndpointThreads(), + callback = callback, + kwargs..., + ) + mask = [i.status == StatusCodes.IntersectedWithGeometry for i in gps] + CoronaGeodesics(trace, m, g, model, gps[mask], source_vels[mask]) +end include("models/lamp-post.jl") diff --git a/src/corona/disc-profiles.jl b/src/corona/disc-profiles.jl index b756f876..a4031e47 100644 --- a/src/corona/disc-profiles.jl +++ b/src/corona/disc-profiles.jl @@ -14,355 +14,5 @@ function delay_flux(profile::AbstractDiscProfile, gps) (delay(profile, gps), emitted_flux(profile, gps)) end -# implementations - -struct RadialDiscProfile{F,R} <: AbstractDiscProfile - # geodesic point to flux - f::F - # geodesic point to time - t::R -end - -function RadialDiscProfile( - m::AbstractMetric, - model::AbstractCoronaModel, - gps::AbstractVector{<:GeodesicPoint}, - source_vels::AbstractVector, - intensities = nothing; - kwargs..., -) - radii = map(i -> i.x[2], gps) - # ensure sorted: let the user sort so that everything is sure to be - # in order - if !issorted(radii) - error( - "geodesic points (and therefore source velocities) must be sorted by radii: use `sortperm(gps; by = i -> i.x[2])` to get the sorting permutation for both", - ) - end - - times = map(i -> i.x[1], gps) - gs = energy_ratio.(m, gps, source_vels) - - RadialDiscProfile(m, radii, times, gs, intensities; kwargs...) -end - -function RadialDiscProfile(rs::AbstractArray, ts::AbstractArray, εs::AbstractArray) - # create interpolations - t = DataInterpolations.LinearInterpolation(ts, rs) - ε = DataInterpolations.LinearInterpolation(εs, rs) - # wrap geodesic point wrappers - RadialDiscProfile(gp -> ε(gp.x[2]), gp -> t(gp.x[2]) + gp.x[1]) -end - -function RadialDiscProfile(rdp::RadialDiscProfile) - Base.depwarn( - "This function is deprecated. Note that `emissivity_profile` now returns a `RadialDiscProfile`.", - :RadialDiscProfile, - ) - rdp -end - -""" -The relativistic correction is calculated via - -```math -\\tilde{A} = A \\sqrt{g_{\\mu,\\nu}(x)} -``` -""" -function RadialDiscProfile( - m::AbstractMetric, - radii::AbstractArray{T}, - times, - gs, - intensity, - ; - grid = GeometricGrid(), - N = 100, -) where {T} - bins = collect(grid(extrema(radii)..., N)) - - ibucket = Buckets.IndexBucket(Int, size(radii)) - bucket!(ibucket, Buckets.Simple(), radii, bins) - groups = Buckets.unpack_bucket(ibucket) - - gs = [@views(mean(gs[g])) for g in groups] - ts = [@views(mean(times[g])) for g in groups] - I = if isnothing(intensity) - # count number of photons in each radial bin - @. convert(T, length(groups)) - else - # sum the intensity over the bin - [@views(sum(intensity[g])) for g in groups] - end - - # interpolate the energy ratio over the disc - bin_domain = bins[1:end-1] - eratios = DataInterpolations.LinearInterpolation(gs, bin_domain) - - for i in eachindex(I) - R = bins[i] - r = i == 1 ? 0 : bins[i-1] - - # 2π comes from integrating over dϕ - dr = R - r - x = SVector(0, R, π / 2, 0) - A = dr * _proper_area(m, x) - # I now stores emissivity - I[i] = source_to_disc_emissivity(m, I[i], A, x, eratios(R)) - end - - RadialDiscProfile(bins, ts, I) -end - -emitted_flux(profile::RadialDiscProfile, gps) = map(profile.f, gps) -delay(profile::RadialDiscProfile, gps) = map(profile.t, gps) -get_emissivity(prof::RadialDiscProfile) = (prof.f.ε.t, prof.f.ε.u) - -struct VoronoiDiscProfile{D,V,G} <: AbstractDiscProfile - disc::D - polys::Vector{Vector{V}} - generators::Vector{V} - geodesic_points::Vector{G} - - function VoronoiDiscProfile( - d::D, - polys::Vector{Vector{V}}, - gen::Vector{V}, - gps::Vector{G}, - ) where {D<:AbstractAccretionDisc,V<:AbstractArray,G} - if !isapprox(d.inclination, π / 2) - return error( - "Currently only supported for discs in the equitorial plane (θ=π/2).", - ) - end - new{D,V,G}(d, polys, gen, gps) - end -end - -function emitted_flux(profile::VoronoiDiscProfile, gps) - error("TODO") -end -function delay(profile::VoronoiDiscProfile, gps) - error("TODO") -end - -# tuple so the calculation may be combined if desired -@inline function delay_flux(profile::VoronoiDiscProfile, gps) - (delay(profile, gps), emitted_flux(profile, gps)) -end - -function Base.show(io::IO, vdp::VoronoiDiscProfile{D}) where {D} - write(io, "VoronoiDiscProfile for $D with $(length(vdp.generators)) generators") -end - -function VoronoiDiscProfile( - m::AbstractMetric{T}, - d::AbstractAccretionDisc{T}, - endpoints::AbstractVector{<:GeodesicPoint{T}}; - padding = 1, -) where {T} - dim = d.outer_radius + padding - rect = VoronoiCells.Rectangle( - GeometryBasics.Point2{T}(-dim, -dim), - GeometryBasics.Point2{T}(dim, dim), - ) - - generators = to_cartesian.(endpoints) - - generators_points = convert.(GeometryBasics.Point2{T}, generators) - - tess = VoronoiCells.voronoicells(generators_points, rect) - polys = GeometryBasics.Polygon.(tess.Cells) - - # cut the polygons down to shape - cutchart!(polys, m, d) - - polys_vecs = unpack_polys(polys) - - VoronoiDiscProfile(d, polys_vecs, generators, endpoints) -end - -function VoronoiDiscProfile( - m::AbstractMetric, - d::AbstractAccretionDisc, - sols::AbstractArray{S}, -) where {S<:SciMLBase.AbstractODESolution} - VoronoiDiscProfile(m, d, map(sol -> unpack_solution(m, sol), sols)) -end - -function VoronoiDiscProfile( - m::AbstractMetric, - d::AbstractAccretionDisc, - simsols::SciMLBase.EnsembleSolution, -) - VoronoiDiscProfile( - m, - d, - filter(i -> i.prob.p.status == StatusCodes.IntersectedWithGeometry, simsols.u), - ) -end - -@noinline function findindex(vdp::VoronoiDiscProfile, p) - for (i, poly) in enumerate(vdp.polys) - # check if we're at all close - let p1 = poly[1] - # todo: what's going on here? - d = @inbounds (p1[1] - p[1])^2 + (p1[2] - p[2])^2 - - if inpolygon(poly, p) - return i - end - end - end - -1 -end - -@inline function findindex(vdp::VoronoiDiscProfile, gp::GeodesicPoint) - p = to_cartesian(gp) - findindex(vdp, p) -end - -function findindex(vdp::VoronoiDiscProfile, gps::AbstractArray{<:GeodesicPoint}; kwargs...) - ret = fill(-1, size(gps)) - Threads.@threads for i in eachindex(gps) - gp = gps[i] - ret[i] = findindex(vdp, gp) - end - ret -end - -getareas(vdp::VoronoiDiscProfile) = getarea.(vdp.polys) - -function getproperarea(poly::AbstractArray, m::AbstractMetric) - A = getarea(poly) - c = getbarycenter(poly) - # get value of metric at the radius of the polygon's barycenter, and in the equitorial plane - m_params = metric_components(m, (sqrt(c[1]^2 + c[2]^2), π / 2)) - # need radial and azimuthal components of the metric - sqrt(m_params[2] * m_params[4]) * A -end - -getproperarea(vdp::VoronoiDiscProfile, m::AbstractMetric) = - map(p -> getproperarea(p, m), vdp.polys) - -function unpack_polys( - polys::AbstractVector{GeometryBasics.Polygon{2,T,GeometryBasics.Point2{T},L,V}}, -) where {T,L,V} - map(polys) do poly - map(SVector{2,T}, getpoints(poly)) - end -end - -""" - getbarycenter(poly) - -Calculate the barycentric point of `poly`. Given that `poly` is an array of 2D points, this -method calculates - -```math -\\vec{c} = \\frac{1}{N} \\sum_i^N \\vec{p}_i, -``` - -where ``\\vec{p}_i`` are the vectors to point ``i`` of the polygon. -""" -function getbarycenter(poly) - xs = sum(i -> first(i), poly) - ys = sum(i -> last(i), poly) - n = length(poly) - SVector{2}(xs / n, ys / n) -end - -""" - cutchart!(polys, m::AbstractMetric{T}, d::AbstractAccretionDisc{T}) - -Trims each polygon in `polys` to match the chart of the disc. Walks the edges of each -polygon until some intersection ``A`` is found. Continues walking until a second intersection -``B`` is found, and then interpolates an arc between ``A`` and ``B``. -""" -function cutchart!(polys, m::AbstractMetric{T}, d::AbstractAccretionDisc{T}) where {T} - rs = inner_radius(m) - inside_radius = rs < d.inner_radius ? d.inner_radius : rs - outside_radius = d.outer_radius - for i in eachindex(polys) - poly = polys[i] - poly = _cut_polygon(inside_radius, poly; inner = true) - poly = _cut_polygon(outside_radius, poly; inner = false) - polys[i] = poly - end - # filter!(!isnothing, polys) -end - -function _cut_polygon(radius, poly; inner = true) - clines = getcycliclines(poly) - (i1, t1) = _circ_path_intersect(radius, clines) - if t1 > 0 - (i2, t2) = _circ_path_intersect(radius, @view(clines[i1+1:end])) - if inner - _circ_cut(radius, clines, i1, t1, i2 + i1, t2; outer = false) - else - _circ_cut(radius, clines, i1, t1, i2 + i1, t2) - end - else - poly - end -end - -function _circ_cut(radius, clines, i1, t1, i2, t2; outer = true) - # get intersection points - A = clines[i1][1] .+ t1 .* (clines[i1][2] - clines[i1][1]) - B = clines[i2][1] .+ t2 .* (clines[i2][2] - clines[i2][1]) - - @assert length(clines) > 2 - - # check if we started in or outside of the circle - dist = √sum(clines[1][1] .^ 2) - if (outer && (dist ≤ radius)) || (!outer && (dist ≥ radius)) - # inside circle - points = first.(@view(clines[1:i1])) - push!(points, A) - push!(points, B) - points = vcat(points, first.(@view(clines[i2+1:end]))) - - return GeometryBasics.Polygon(points) - else - # outside circle - points = first.(@view(clines[i1+1:i2])) - push!(points, B) - push!(points, A) - - return GeometryBasics.Polygon(points) - end -end - -function _circ_path_intersect(radius, lines::AbstractArray{<:GeometryBasics.Line}) - for (i, line) in enumerate(lines) - t = _circ_line_intersect(radius, line) - if t > 0 - return i, t - end - end - -1, -1 -end - -@fastmath function _circ_line_intersect(radius, line::GeometryBasics.Line) - let A = line.points[1], B = line.points[2] - d = B .- A - d2 = sum(d .* d) - A2 = sum(A .* A) - da = sum(d .* A) - Δ = da^2 - d2 * (A2 - radius^2) - - if Δ > 0 - tp = (-da + √Δ) / d2 - tn = (-da - √Δ) / d2 - - if (1 ≥ tp ≥ 0) - return tp - end - if (1 ≥ tn ≥ 0) - return tn - end - end - end - -1 -end +include("profiles/radial.jl") +include("profiles/voronoi.jl") diff --git a/src/corona/emissivity.jl b/src/corona/emissivity.jl index c718ddd1..161376a1 100644 --- a/src/corona/emissivity.jl +++ b/src/corona/emissivity.jl @@ -20,82 +20,6 @@ function source_to_disc_emissivity(δ, g, A, γ; Γ = 2) sin(δ) / (g^Γ * A * γ) end -struct CoronaGeodesics{T,M,G,C,P,V} - trace::T - metric::M - geometry::G - model::C - geodesic_points::P - source_velocity::V -end - -function tracecorona( - m::AbstractMetric, - g::AbstractAccretionGeometry, - model::AbstractCoronaModel; - λ_max = 10_000, - n_samples = 1024, - sampler = EvenSampler(domain = BothHemispheres(), generator = RandomGenerator()), - trace = TraceGeodesic(), - callback = domain_upper_hemisphere(), - kwargs..., -) - xs, vs, source_vels = sample_position_direction_velocity(m, model, sampler, n_samples) - gps = tracegeodesics( - m, - xs, - vs, - g, - λ_max; - trace = trace, - save_on = false, - ensemble = EnsembleEndpointThreads(), - callback = callback, - kwargs..., - ) - mask = [i.status == StatusCodes.IntersectedWithGeometry for i in gps] - CoronaGeodesics(trace, m, g, model, gps[mask], source_vels[mask]) -end - -function RadialDiscProfile(ce::CoronaGeodesics; kwargs...) - J = sortperm(ce.geodesic_points; by = i -> i.x[2]) - @views RadialDiscProfile( - ce.metric, - ce.model, - ce.geodesic_points[J], - ce.source_velocity[J], - ) -end - -function RadialDiscProfile(ce::CoronaGeodesics{<:TraceRadiativeTransfer}; kwargs...) - J = sortperm(ce.geodesic_points; by = i -> i.x[2]) - @views RadialDiscProfile( - ce.metric, - ce.model, - ce.geodesic_points[J], - ce.source_velocity[J], - [i.aux[1] for i in ce.geodesic_points[J]], - ) -end - -function RadialDiscProfile(ε, ce::CoronaGeodesics; kwargs...) - J = sortperm(ce.geodesic_points; by = i -> i.x[2]) - radii = @views [i.x[2] for i in ce.geodesic_points[J]] - times = @views [i.x[1] for i in ce.geodesic_points[J]] - - t = DataInterpolations.LinearInterpolation(times, radii) - - function _emissivity_wrapper(gp) - ε(gp.x[2]) - end - - function _delay_wrapper(gp) - t(gp.x[2]) + gp.x[1] - end - - RadialDiscProfile(_emissivity_wrapper, _delay_wrapper) -end - """ function emissivity_profile( m::AbstractMetric, @@ -133,6 +57,21 @@ function emissivity_profile( N = N, ) end +function emissivity_profile( + ::Nothing, + m::AbstractMetric, + d::AbstractAccretionGeometry, + model::AbstractCoronaModel; + kwargs..., +) + error( + """Not yet implemented for $(Base.typename(typeof(model)).name) with $(Base.typename(typeof(d)).name). + This dispatch is reserved for more sophisticated (rather that just pure Monte-Carlo sampling) + of the emissivity profile, and may not be applicable for all coronal models. Pass the `sampler` kwarg + with an `AbstractDirectionSampler` to use the general strategy. + """, + ) +end function _proper_area(m, x::SVector{4}) gcomp = Gradus.metric_components(m, SVector(x[2], x[3])) @@ -153,9 +92,11 @@ function _point_source_symmetric_emissivity_profile( n_samples = 100, λ_max = 10_000.0, callback = domain_upper_hemisphere(), + δmin = 0.1, + δmax = 179.9, kwargs..., ) - δs = deg2rad.(range(0.1, 179.9, n_samples)) + δs = deg2rad.(range(δmin, δmax, n_samples)) # we assume a point source x, v = sample_position_velocity(m, model) velfunc = polar_angle_to_velfunc(m, x, v, δs) diff --git a/src/corona/profiles/radial.jl b/src/corona/profiles/radial.jl new file mode 100644 index 00000000..489b5f5c --- /dev/null +++ b/src/corona/profiles/radial.jl @@ -0,0 +1,140 @@ +struct RadialDiscProfile{F,R} <: AbstractDiscProfile + # geodesic point to flux + f::F + # geodesic point to time + t::R +end + +function RadialDiscProfile( + m::AbstractMetric, + model::AbstractCoronaModel, + gps::AbstractVector{<:GeodesicPoint}, + source_vels::AbstractVector, + intensities = nothing; + kwargs..., +) + radii = map(i -> i.x[2], gps) + # ensure sorted: let the user sort so that everything is sure to be + # in order + if !issorted(radii) + error( + "geodesic points (and therefore source velocities) must be sorted by radii: use `sortperm(gps; by = i -> i.x[2])` to get the sorting permutation for both", + ) + end + + times = map(i -> i.x[1], gps) + gs = energy_ratio.(m, gps, source_vels) + + RadialDiscProfile(m, radii, times, gs, intensities; kwargs...) +end + +function RadialDiscProfile(rs::AbstractArray, ts::AbstractArray, εs::AbstractArray) + # create interpolations + t = DataInterpolations.LinearInterpolation(ts, rs) + ε = DataInterpolations.LinearInterpolation(εs, rs) + # wrap geodesic point wrappers + RadialDiscProfile(gp -> ε(gp.x[2]), gp -> t(gp.x[2]) + gp.x[1]) +end + +function RadialDiscProfile(rdp::RadialDiscProfile) + Base.depwarn( + "This function is deprecated. Note that `emissivity_profile` now returns a `RadialDiscProfile`.", + :RadialDiscProfile, + ) + rdp +end + +""" +The relativistic correction is calculated via + +```math +\\tilde{A} = A \\sqrt{g_{\\mu,\\nu}(x)} +``` +""" +function RadialDiscProfile( + m::AbstractMetric, + radii::AbstractArray{T}, + times, + gs, + intensity, + ; + grid = GeometricGrid(), + N = 100, +) where {T} + bins = collect(grid(extrema(radii)..., N)) + + ibucket = Buckets.IndexBucket(Int, size(radii)) + bucket!(ibucket, Buckets.Simple(), radii, bins) + groups = Buckets.unpack_bucket(ibucket) + + gs = [@views(mean(gs[g])) for g in groups] + ts = [@views(mean(times[g])) for g in groups] + I = if isnothing(intensity) + # count number of photons in each radial bin + @. convert(T, length(groups)) + else + # sum the intensity over the bin + [@views(sum(intensity[g])) for g in groups] + end + + # interpolate the energy ratio over the disc + bin_domain = bins[1:end-1] + eratios = DataInterpolations.LinearInterpolation(gs, bin_domain) + + for i in eachindex(I) + R = bins[i] + r = i == 1 ? 0 : bins[i-1] + + # 2π comes from integrating over dϕ + dr = R - r + x = SVector(0, R, π / 2, 0) + A = dr * _proper_area(m, x) + # I now stores emissivity + I[i] = source_to_disc_emissivity(m, I[i], A, x, eratios(R)) + end + + RadialDiscProfile(bins, ts, I) +end + +function RadialDiscProfile(ce::CoronaGeodesics; kwargs...) + J = sortperm(ce.geodesic_points; by = i -> i.x[2]) + @views RadialDiscProfile( + ce.metric, + ce.model, + ce.geodesic_points[J], + ce.source_velocity[J], + ) +end + +function RadialDiscProfile(ce::CoronaGeodesics{<:TraceRadiativeTransfer}; kwargs...) + J = sortperm(ce.geodesic_points; by = i -> i.x[2]) + @views RadialDiscProfile( + ce.metric, + ce.model, + ce.geodesic_points[J], + ce.source_velocity[J], + [i.aux[1] for i in ce.geodesic_points[J]], + ) +end + +function RadialDiscProfile(ε, ce::CoronaGeodesics; kwargs...) + J = sortperm(ce.geodesic_points; by = i -> i.x[2]) + radii = @views [i.x[2] for i in ce.geodesic_points[J]] + times = @views [i.x[1] for i in ce.geodesic_points[J]] + + t = DataInterpolations.LinearInterpolation(times, radii) + + function _emissivity_wrapper(gp) + ε(gp.x[2]) + end + + function _delay_wrapper(gp) + t(gp.x[2]) + gp.x[1] + end + + RadialDiscProfile(_emissivity_wrapper, _delay_wrapper) +end + +emitted_flux(profile::RadialDiscProfile, gps) = map(profile.f, gps) +delay(profile::RadialDiscProfile, gps) = map(profile.t, gps) +get_emissivity(prof::RadialDiscProfile) = (prof.f.ε.t, prof.f.ε.u) diff --git a/src/corona/profiles/voronoi.jl b/src/corona/profiles/voronoi.jl new file mode 100644 index 00000000..e2931612 --- /dev/null +++ b/src/corona/profiles/voronoi.jl @@ -0,0 +1,248 @@ +struct VoronoiDiscProfile{D,V,G} <: AbstractDiscProfile + disc::D + polys::Vector{Vector{V}} + generators::Vector{V} + geodesic_points::Vector{G} + + function VoronoiDiscProfile( + d::D, + polys::Vector{Vector{V}}, + gen::Vector{V}, + gps::Vector{G}, + ) where {D<:AbstractAccretionDisc,V<:AbstractArray,G} + if !isapprox(d.inclination, π / 2) + return error( + "Currently only supported for discs in the equitorial plane (θ=π/2).", + ) + end + new{D,V,G}(d, polys, gen, gps) + end +end + +function emitted_flux(profile::VoronoiDiscProfile, gps) + error("TODO") +end +function delay(profile::VoronoiDiscProfile, gps) + error("TODO") +end + +# tuple so the calculation may be combined if desired +@inline function delay_flux(profile::VoronoiDiscProfile, gps) + (delay(profile, gps), emitted_flux(profile, gps)) +end + +function Base.show(io::IO, vdp::VoronoiDiscProfile{D}) where {D} + write(io, "VoronoiDiscProfile for $D with $(length(vdp.generators)) generators") +end + +function VoronoiDiscProfile( + m::AbstractMetric{T}, + d::AbstractAccretionDisc{T}, + endpoints::AbstractVector{<:GeodesicPoint{T}}; + padding = 1, +) where {T} + dim = d.outer_radius + padding + rect = VoronoiCells.Rectangle( + GeometryBasics.Point2{T}(-dim, -dim), + GeometryBasics.Point2{T}(dim, dim), + ) + + generators = to_cartesian.(endpoints) + + generators_points = convert.(GeometryBasics.Point2{T}, generators) + + tess = VoronoiCells.voronoicells(generators_points, rect) + polys = GeometryBasics.Polygon.(tess.Cells) + + # cut the polygons down to shape + cutchart!(polys, m, d) + + polys_vecs = unpack_polys(polys) + + VoronoiDiscProfile(d, polys_vecs, generators, endpoints) +end + +function VoronoiDiscProfile( + m::AbstractMetric, + d::AbstractAccretionDisc, + sols::AbstractArray{S}, +) where {S<:SciMLBase.AbstractODESolution} + VoronoiDiscProfile(m, d, map(sol -> unpack_solution(m, sol), sols)) +end + +function VoronoiDiscProfile( + m::AbstractMetric, + d::AbstractAccretionDisc, + simsols::SciMLBase.EnsembleSolution, +) + VoronoiDiscProfile( + m, + d, + filter(i -> i.prob.p.status == StatusCodes.IntersectedWithGeometry, simsols.u), + ) +end + +@noinline function findindex(vdp::VoronoiDiscProfile, p) + for (i, poly) in enumerate(vdp.polys) + # check if we're at all close + let p1 = poly[1] + # todo: what's going on here? + d = @inbounds (p1[1] - p[1])^2 + (p1[2] - p[2])^2 + + if inpolygon(poly, p) + return i + end + end + end + -1 +end + +@inline function findindex(vdp::VoronoiDiscProfile, gp::GeodesicPoint) + p = to_cartesian(gp) + findindex(vdp, p) +end + +function findindex(vdp::VoronoiDiscProfile, gps::AbstractArray{<:GeodesicPoint}; kwargs...) + ret = fill(-1, size(gps)) + Threads.@threads for i in eachindex(gps) + gp = gps[i] + ret[i] = findindex(vdp, gp) + end + ret +end + +getareas(vdp::VoronoiDiscProfile) = getarea.(vdp.polys) + +function getproperarea(poly::AbstractArray, m::AbstractMetric) + A = getarea(poly) + c = getbarycenter(poly) + # get value of metric at the radius of the polygon's barycenter, and in the equitorial plane + m_params = metric_components(m, (sqrt(c[1]^2 + c[2]^2), π / 2)) + # need radial and azimuthal components of the metric + sqrt(m_params[2] * m_params[4]) * A +end + +getproperarea(vdp::VoronoiDiscProfile, m::AbstractMetric) = + map(p -> getproperarea(p, m), vdp.polys) + +function unpack_polys( + polys::AbstractVector{GeometryBasics.Polygon{2,T,GeometryBasics.Point2{T},L,V}}, +) where {T,L,V} + map(polys) do poly + map(SVector{2,T}, getpoints(poly)) + end +end + +""" + getbarycenter(poly) + +Calculate the barycentric point of `poly`. Given that `poly` is an array of 2D points, this +method calculates + +```math +\\vec{c} = \\frac{1}{N} \\sum_i^N \\vec{p}_i, +``` + +where ``\\vec{p}_i`` are the vectors to point ``i`` of the polygon. +""" +function getbarycenter(poly) + xs = sum(i -> first(i), poly) + ys = sum(i -> last(i), poly) + n = length(poly) + SVector{2}(xs / n, ys / n) +end + +""" + cutchart!(polys, m::AbstractMetric{T}, d::AbstractAccretionDisc{T}) + +Trims each polygon in `polys` to match the chart of the disc. Walks the edges of each +polygon until some intersection ``A`` is found. Continues walking until a second intersection +``B`` is found, and then interpolates an arc between ``A`` and ``B``. +""" +function cutchart!(polys, m::AbstractMetric{T}, d::AbstractAccretionDisc{T}) where {T} + rs = inner_radius(m) + inside_radius = rs < d.inner_radius ? d.inner_radius : rs + outside_radius = d.outer_radius + for i in eachindex(polys) + poly = polys[i] + poly = _cut_polygon(inside_radius, poly; inner = true) + poly = _cut_polygon(outside_radius, poly; inner = false) + polys[i] = poly + end + # filter!(!isnothing, polys) +end + +function _cut_polygon(radius, poly; inner = true) + clines = getcycliclines(poly) + (i1, t1) = _circ_path_intersect(radius, clines) + if t1 > 0 + (i2, t2) = _circ_path_intersect(radius, @view(clines[i1+1:end])) + if inner + _circ_cut(radius, clines, i1, t1, i2 + i1, t2; outer = false) + else + _circ_cut(radius, clines, i1, t1, i2 + i1, t2) + end + else + poly + end +end + +function _circ_cut(radius, clines, i1, t1, i2, t2; outer = true) + # get intersection points + A = clines[i1][1] .+ t1 .* (clines[i1][2] - clines[i1][1]) + B = clines[i2][1] .+ t2 .* (clines[i2][2] - clines[i2][1]) + + @assert length(clines) > 2 + + # check if we started in or outside of the circle + dist = √sum(clines[1][1] .^ 2) + if (outer && (dist ≤ radius)) || (!outer && (dist ≥ radius)) + # inside circle + points = first.(@view(clines[1:i1])) + push!(points, A) + push!(points, B) + points = vcat(points, first.(@view(clines[i2+1:end]))) + + return GeometryBasics.Polygon(points) + else + # outside circle + points = first.(@view(clines[i1+1:i2])) + push!(points, B) + push!(points, A) + + return GeometryBasics.Polygon(points) + end +end + +function _circ_path_intersect(radius, lines::AbstractArray{<:GeometryBasics.Line}) + for (i, line) in enumerate(lines) + t = _circ_line_intersect(radius, line) + if t > 0 + return i, t + end + end + -1, -1 +end + +@fastmath function _circ_line_intersect(radius, line::GeometryBasics.Line) + let A = line.points[1], B = line.points[2] + d = B .- A + d2 = sum(d .* d) + A2 = sum(A .* A) + da = sum(d .* A) + Δ = da^2 - d2 * (A2 - radius^2) + + if Δ > 0 + tp = (-da + √Δ) / d2 + tn = (-da - √Δ) / d2 + + if (1 ≥ tp ≥ 0) + return tp + end + if (1 ≥ tn ≥ 0) + return tn + end + end + end + -1 +end diff --git a/src/corona/sky-geometry.jl b/src/corona/samplers.jl similarity index 100% rename from src/corona/sky-geometry.jl rename to src/corona/samplers.jl From 19a24a9d192f5e01551aa926f0b9fa137f9d9e68 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Tue, 25 Jul 2023 08:54:16 +0100 Subject: [PATCH 7/9] fix: forward kwargs in RadialDiscProfile --- src/corona/profiles/radial.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/corona/profiles/radial.jl b/src/corona/profiles/radial.jl index 489b5f5c..cf763394 100644 --- a/src/corona/profiles/radial.jl +++ b/src/corona/profiles/radial.jl @@ -102,7 +102,8 @@ function RadialDiscProfile(ce::CoronaGeodesics; kwargs...) ce.metric, ce.model, ce.geodesic_points[J], - ce.source_velocity[J], + ce.source_velocity[J]; + kwargs..., ) end @@ -113,7 +114,8 @@ function RadialDiscProfile(ce::CoronaGeodesics{<:TraceRadiativeTransfer}; kwargs ce.model, ce.geodesic_points[J], ce.source_velocity[J], - [i.aux[1] for i in ce.geodesic_points[J]], + [i.aux[1] for i in ce.geodesic_points[J]]; + kwargs..., ) end @@ -132,7 +134,7 @@ function RadialDiscProfile(ε, ce::CoronaGeodesics; kwargs...) t(gp.x[2]) + gp.x[1] end - RadialDiscProfile(_emissivity_wrapper, _delay_wrapper) + RadialDiscProfile(_emissivity_wrapper, _delay_wrapper; kwargs...) end emitted_flux(profile::RadialDiscProfile, gps) = map(profile.f, gps) From 8213fa7c13aec5f6fd5f6b0e43c1273fa94b9c5e Mon Sep 17 00:00:00 2001 From: fjebaker Date: Tue, 25 Jul 2023 11:45:15 +0100 Subject: [PATCH 8/9] test: added numeric tests for emissivity --- test/runtests.jl | 1 + test/unit/emissivity.jl | 46 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+) create mode 100644 test/unit/emissivity.jl diff --git a/test/runtests.jl b/test/runtests.jl index 570b898b..1c0417cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,6 +21,7 @@ end @time @testset "corona" verbose = true begin include("unit/coronal-beaming.jl") + include("unit/emissivity.jl") end @time @testset "integration" verbose = true begin diff --git a/test/unit/emissivity.jl b/test/unit/emissivity.jl new file mode 100644 index 00000000..43bff8bd --- /dev/null +++ b/test/unit/emissivity.jl @@ -0,0 +1,46 @@ +using Test +using Gradus + +m = KerrMetric(1.0, 0.998) +model = LampPostModel(h = 10.0) +d = GeometricThinDisc(0.0, 500.0, π / 2) + +# using point source angular method +profile = emissivity_profile(m, d, model; n_samples = 20) + +r, em = get_emissivity(profile) +@test em ≈ [ + 0.03530163295788805, + 0.01618161075710721, + 0.010034247429248367, + 0.006231511729195148, + 0.0035296670918803694, + 0.0016859014503211888, + 0.0005946901244498561, + 0.0001052324320767813, +] atol = 1e-5 + +# using generic monte-carlo sampling +profile = emissivity_profile( + m, + d, + model; + n_samples = 1000, + sampler = EvenSampler(BothHemispheres(), GoldenSpiralGenerator()), + N = 10, +) + +r, em = get_emissivity(profile) + +@test em ≈ [ + 1.4536471470344825, + 3.2159234810191997, + 1.9624518676973357, + 0.8654137966466704, + 0.23973494196571282, + 0.046232848740113686, + 0.00789393186371828, + 0.0012974087513837867, + 0.0001872022563608107, + 3.825853506250807e-5, +] atol = 1e-5 From 9c3833d7c668902d5228a09d977791d59adb1af0 Mon Sep 17 00:00:00 2001 From: fjebaker Date: Tue, 25 Jul 2023 11:45:58 +0100 Subject: [PATCH 9/9] chore: update buckets dep --- Project.toml | 2 +- src/corona/profiles/radial.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 2ca4566f..b980b8dc 100644 --- a/Project.toml +++ b/Project.toml @@ -29,5 +29,5 @@ Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" VoronoiCells = "e3e34ffb-84e9-5012-9490-92c94d0c60a4" [compat] -Buckets = ">=0.1.10" +Buckets = ">=0.1.11" julia = "1.9" diff --git a/src/corona/profiles/radial.jl b/src/corona/profiles/radial.jl index cf763394..54c49922 100644 --- a/src/corona/profiles/radial.jl +++ b/src/corona/profiles/radial.jl @@ -63,7 +63,7 @@ function RadialDiscProfile( ) where {T} bins = collect(grid(extrema(radii)..., N)) - ibucket = Buckets.IndexBucket(Int, size(radii)) + ibucket = Buckets.IndexBucket(Int, size(radii), length(bins)) bucket!(ibucket, Buckets.Simple(), radii, bins) groups = Buckets.unpack_bucket(ibucket)