Merge pull request #180 from phajy/main

Add extrapolation to orbit-solving interpolation

src/corona/profiles/analytic.jl

@@ -16,10 +16,22 @@ function AnalyticRadialDiscProfile(emissivity, cg::CoronaGeodesics)
     AnalyticRadialDiscProfile(emissivity, t)
 end
 
-emissivity_at(prof::AnalyticRadialDiscProfile, r::Number) = prof.ε(r)
+function emissivity_at(prof::AnalyticRadialDiscProfile, r::Number)
+    r_bounded = _enforce_interpolation_bounds(r, prof)
+    prof.ε(r)
+end
 emissivity_at(prof::AnalyticRadialDiscProfile, gp::AbstractGeodesicPoint) =
     emissivity_at(prof, _equatorial_project(gp.x))
 
-coordtime_at(prof::AnalyticRadialDiscProfile, r::Number) = prof.t(r)
+function coordtime_at(prof::AnalyticRadialDiscProfile, r::Number)
+    r_bounded = _enforce_interpolation_bounds(r, prof)
+    prof.t(r_bounded)
+end
 coordtime_at(prof::AnalyticRadialDiscProfile, gp::AbstractGeodesicPoint) =
     coordtime_at(prof, _equatorial_project(gp.x)) + gp.x[1]
+
+function _enforce_interpolation_bounds(r::Number, prof::AnalyticRadialDiscProfile)
+    r_min = first(prof.t.t)
+    r_max = last(prof.t.t)
+    _enforce_interpolation_bounds(r, r_min, r_max)
+end

src/corona/profiles/radial.jl

@@ -6,11 +6,23 @@ struct RadialDiscProfile{T,I} <: AbstractDiscProfile
     interp_t::I
 end
 
-emissivity_at(prof::RadialDiscProfile, r::Number) = prof.interp_ε(r)
+@inline function _enforce_interpolation_bounds(r::Number, prof::RadialDiscProfile)
+    r_min = first(prof.radii)
+    r_max = last(prof.radii)
+    return _enforce_interpolation_bounds(r, r_min, r_max)
+end
+
+function emissivity_at(prof::RadialDiscProfile, r::Number)
+    r_bounded = _enforce_interpolation_bounds(r, prof)
+    prof.interp_ε(r_bounded)
+end
 emissivity_at(prof::RadialDiscProfile, gp::AbstractGeodesicPoint) =
     emissivity_at(prof, _equatorial_project(gp.x))
 
-coordtime_at(prof::RadialDiscProfile, r::Number) = prof.interp_t(r)
+function coordtime_at(prof::RadialDiscProfile, r::Number)
+    r_bounded = _enforce_interpolation_bounds(r, prof)
+    prof.interp_t(r_bounded)
+end
 coordtime_at(prof::RadialDiscProfile, gp::AbstractGeodesicPoint) =
     coordtime_at(prof, _equatorial_project(gp.x)) + gp.x[1]
 

src/orbits/orbit-solving.jl

@@ -112,13 +112,12 @@ struct PlungingInterpolation{M,_interp_type}
         vr = sol[6, :][I]
         vϕ = sol[8, :][I]
 
-        rinterp = DataInterpolations.LinearInterpolation(vt, r)
-
-        new{M,typeof(rinterp)}(
+        r_interp = _make_interpolation(r, vt)
+        new{M,typeof(r_interp)}(
             m,
-            rinterp,
-            DataInterpolations.LinearInterpolation(vr, r),
-            DataInterpolations.LinearInterpolation(vϕ, r),
+            r_interp,
+            _make_interpolation(r, vr),
+            _make_interpolation(r, vϕ),
         )
     end
 end
@@ -128,9 +127,10 @@ function Base.show(io::IO, ::PlungingInterpolation{M}) where {M}
 end
 
 function (pintrp::PlungingInterpolation)(r)
-    vt = pintrp.t(r)
-    vr = pintrp.r(r)
-    vϕ = pintrp.ϕ(r)
+    r_bounded = _enforce_interpolation_bounds(r, pintrp)
+    vt = pintrp.t(r_bounded)
+    vr = pintrp.r(r_bounded)
+    vϕ = pintrp.ϕ(r_bounded)
     SVector(vt, vr, 0, vϕ)
 end
 
@@ -166,5 +166,9 @@ function interpolate_plunging_velocities(
     PlungingInterpolation(m, sol)
 end
 
+function _enforce_interpolation_bounds(r::Number, pintrp::PlungingInterpolation)
+    _enforce_interpolation_bounds(r, first(pintrp.r.t), last(pintrp.r.t))
+end
+
 export solve_equatorial_circular_orbit,
     trace_equatorial_circular_orbit, PlungingInterpolation, interpolate_plunging_velocities

src/transfer-functions/transfer-functions-2d.jl

@@ -174,21 +174,21 @@ function lagtransfer(
 end
 
 function binflux(tf::LagTransferFunction; kwargs...)
-    profile = AnalyticRadialDiscProfile(r -> r^-3, tf.emissivity_profile)
+    profile = AnalyticRadialDiscProfile(r -> r^-3, tf.coronal_geodesics)
     binflux(tf, profile; kwargs...)
 end
 
 function binflux(
     tf::LagTransferFunction,
     profile::AbstractDiscProfile;
-    redshift = ConstPointFunctions.redshift(tf.emissivity_profile.metric, tf.x),
+    redshift = ConstPointFunctions.redshift(tf.coronal_geodesics.metric, tf.x),
     E₀ = 6.4,
     t0 = tf.x[2],
     kwargs...,
 )
     t = coordtime_at(profile, tf.observer_to_disc)
     ε = emissivity_at(profile, tf.observer_to_disc)
-    g = redshift.(tf.emissivity_profile.metric, tf.observer_to_disc, tf.max_t)
+    g = redshift.(tf.coronal_geodesics.metric, tf.observer_to_disc, tf.max_t)
     # calculate flux
     f = @. g^3 * ε * tf.image_plane_areas
     # normalize

src/transfer-functions/types.jl

@@ -32,17 +32,17 @@ struct LagTransferFunction{T,X,E,P}
     max_t::T
     x::X
     image_plane_areas::Vector{T}
-    emissivity_profile::E
+    coronal_geodesics::E
     observer_to_disc::Vector{P}
 end
 
 function Base.show(io::IO, ::MIME"text/plain", tf::LagTransferFunction)
-    text = """LagTransferFunction for $(typeof(tf.emissivity_profile.metric)) 
+    text = """LagTransferFunction for $(typeof(tf.coronal_geodesics.metric)) 
       . observer position      
           $(tf.x)
-      . model                         : $(typeof(tf.emissivity_profile.model))
+      . model                         : $(typeof(tf.coronal_geodesics.model))
       . observer to disc photon count : $(length(tf.observer_to_disc))
-      . source to disc photon count   : $(length(tf.emissivity_profile.geodesic_points))
+      . source to disc photon count   : $(length(tf.coronal_geodesics.geodesic_points))
       Total memory: $(Base.format_bytes(Base.summarysize(tf)))
     """
     print(io, text)

src/utils.jl

@@ -36,6 +36,13 @@ function _make_interpolation(x, y)
     DataInterpolations.LinearInterpolation(y, x)
 end
 
+@inline function _enforce_interpolation_bounds(r::Number, r_min::Number, r_max::Number)
+    if (r < r_min) || (r > r_max)
+        @warn "Interpolation out of bounds $r ∉ [$(r_min),  $(r_max)]. Additional geodesic samples may be required."
+    end
+    clamp(r, r_min, r_max)
+end
+
 @inline function _linear_interpolate(y1, y2, θ)
     (1 - θ) * y1 + θ * y2
 end

test/transfer-functions/test-2d.jl

@@ -23,7 +23,7 @@ tf = lagtransfer(
 
 # check number of intersection points
 @test length(tf.observer_to_disc) == 335
-@test length(tf.emissivity_profile.geodesic_points) == 58
+@test length(tf.coronal_geodesics.geodesic_points) == 58
 
 # ensure binning works as expected
 t, E, f = binflux(tf, N_t = 100, N_E = 100)