From ab8c441425481361983960e16df84765bf5d75d6 Mon Sep 17 00:00:00 2001 From: LuEdRaMo <73906617+LuEdRaMo@users.noreply.github.com> Date: Thu, 10 Oct 2024 18:16:01 -0600 Subject: [PATCH] Introduce `ODProblem` (#87) * Add errormodel.jl * Add odproblem.jl * Update process_radec.jl * Update parameters.jl * Update orbit determination functions * Update orbit determination tests * Update propagation tests * Update bplane tests * Update scripts/orbitdetermination.jl * Minor additions * New PropagationBuffer constructor * Update scripts/mcmov.jl * Update examples/recipes.jl * Bump minor version * Update scripts/adam.jl * Minor fix --- Project.toml | 2 +- examples/recipes.jl | 16 +- scripts/adam.jl | 52 +-- scripts/mcmov.jl | 24 +- scripts/orbitdetermination.jl | 8 +- src/NEOs.jl | 9 +- src/observations/errormodel.jl | 421 +++++++++++++++++++ src/observations/observations.jl | 1 + src/observations/process_radec.jl | 370 ++++------------ src/observations/tracklet.jl | 21 +- src/orbitdetermination/gaussinitcond.jl | 67 ++- src/orbitdetermination/neosolution.jl | 23 +- src/orbitdetermination/odproblem.jl | 91 ++++ src/orbitdetermination/orbitdetermination.jl | 126 +++--- src/orbitdetermination/tooshortarc.jl | 89 ++-- src/postprocessing/outlier_rejection.jl | 46 +- src/propagation/parameters.jl | 190 ++++----- test/bplane.jl | 9 +- test/orbitdetermination.jl | 69 +-- test/propagation.jl | 20 +- 20 files changed, 962 insertions(+), 692 deletions(-) create mode 100644 src/observations/errormodel.jl create mode 100644 src/orbitdetermination/odproblem.jl diff --git a/Project.toml b/Project.toml index beb4d61b..3e6f3072 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NEOs" uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df" authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"] -version = "0.9.1" +version = "0.10.0" [deps] Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" diff --git a/examples/recipes.jl b/examples/recipes.jl index f82cebf6..256eb772 100644 --- a/examples/recipes.jl +++ b/examples/recipes.jl @@ -2,23 +2,24 @@ # ext/NEOsRecipesBaseExt.jl using NEOs, Plots -# Download optical astrometry of asteroid 2023 DW -radec = fetch_radec_mpc("2023 DW") +# Download optical astrometry of asteroid 2016 TU93 +radec = fetch_radec_mpc("2016 TU93") # Orbit determination params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007) -sol = orbitdetermination(radec, params) +od = ODProblem(newtonian!, radec) +sol = orbitdetermination(od, params) # Plot optical residuals plot(sol.res, label = "", xlabel = "αcos(δ) [arcsec]", ylabel = "δ [arcsec]", - title = "2023 DW optical O-C residuals") + title = "2016 TU93 optical O-C residuals") # Plot asteroid's orbit (in barycentric cartesian coordinates) # N: number of points between t0 and tf (100 by default) # projection: which coordinates to plot, options are: :x, :y, :z, :xy, # :xz, :yz and :xyz (default) t0, tf = sol.bwd.t0 + sol.bwd.t[end], sol.fwd.t0 + sol.fwd.t[end] -plot(sol, t0, tf, N = 10_000, projection = :xyz, label = "2023 DW", +plot(sol, t0, tf, N = 10_000, projection = :xyz, label = "2016 TU93", xlabel = "x [au]", ylabel = "y [au]", zlabel = "z [au]", color = :blue) # We can also visualize the orbits of the Sun and the Earth @@ -33,6 +34,5 @@ plot!(params.eph_ea, t0, tf, N = 10_000, projection = :xyz, using NEOs: AdmissibleRegion A = AdmissibleRegion(sol.tracklets[1], params) plot(A, N = 10_000, ρscale = :log, label = "", xlabel = "log10(ρ) [au]", - ylabel = "v_ρ [au/day]", title = "2023 DW first tracklet AR", - c = :red, linewidth = 2) - + ylabel = "v_ρ [au/day]", title = "2016 TU93 first tracklet AR", + c = :red, linewidth = 2) \ No newline at end of file diff --git a/scripts/adam.jl b/scripts/adam.jl index 4d45dff1..d8f1a580 100644 --- a/scripts/adam.jl +++ b/scripts/adam.jl @@ -43,21 +43,21 @@ end @everywhere begin using NEOs, Dates, TaylorSeries, PlanetaryEphemeris, JLD2 using NEOs: RadecMPC, AdmissibleRegion, PropagationBuffer, OpticalResidual, - attr2bary, propres!, boundary_projection, reduce_tracklets, arboundary + attr2bary, propres!, boundary_projection, reduce_tracklets, arboundary, + indices - function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, + function adam(od::ODProblem{D, T}, i::Int, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, params::NEOParameters{T}; scale::Symbol = :linear, η::T = 25.0, - μ::T = 0.75, ν::T = 0.9, ϵ::T = 1e-8, Qtol::T = 0.001, adamorder::Int = 2, - dynamics::D = newtonian!) where {T <: Real, D} - # Initial time of integration [julian days] + μ::T = 0.75, ν::T = 0.9, ϵ::T = 1e-8, adamorder::Int = 2) where {D, T <: Real} + # Initial time of integration [julian days TDB] jd0 = dtutc2jdtdb(A.date) - # Maximum number of iterations - maxiter = params.adamiter + # Unfold maximum number of iterations and relative tolerance + maxiter, Qtol = params.adamiter, params.adamQtol # Allocate memory aes = Matrix{T}(undef, 6, maxiter+1) Qs = fill(T(Inf), maxiter+1) # Initial attributable elements - aes[:, 1] .= [A.α, A.δ, A.v_α, A.v_δ, ρ, v_ρ] + aes[:, 1] .= A.α, A.δ, A.v_α, A.v_δ, ρ, v_ρ # Scaling factors scalings = Vector{T}(undef, 6) scalings[1:4] .= abs.(aes[1:4, 1]) ./ 1e6 @@ -67,46 +67,44 @@ end scalings[5] = (log10(A.ρ_domain[2]) - log10(A.ρ_domain[1])) / 1_000 end scalings[6] = (A.v_ρ_domain[2] - A.v_ρ_domain[1]) / 1_000 - # Jet transport variables - set_variables(Float64, "dx"; order = adamorder, numvars = 6) + # Jet transport variables and initial condition dae = [scalings[i] * TaylorN(i, order = adamorder) for i in 1:6] + AE = aes[:, 1] .+ dae + # Subset of radec + idxs = indices(od.tracklets[i]) # Propagation buffer - t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end])) - tlim = (t0 - params.bwdoffset, tf + params.fwdoffset) - buffer = PropagationBuffer(dynamics, jd0, tlim, aes[:, 1] .+ dae, params) + buffer = PropagationBuffer(od, jd0, idxs[1], idxs[end], AE, params) # Vector of O-C residuals - res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec)) + res = [zero(OpticalResidual{T, TaylorN{T}}) for _ in eachindex(idxs)] # Origin - x0 = zeros(T, 6) - x1 = zeros(T, 6) + x0, x1 = zeros(T, 6), zeros(T, 6) # Gradient of objective function wrt (ρ, v_ρ) g_t = Vector{T}(undef, 2) # First momentum - m = zeros(T, 2) - _m_ = zeros(T, 2) + m, _m_ = zeros(T, 2), zeros(T, 2) # Second momentum - n = zeros(T, 2) - _n_ = zeros(T, 2) + n, _n_ = zeros(T, 2), zeros(T, 2) # Gradient descent for t in 1:maxiter # Current attributable elements (plain) ae = aes[:, t] # Attributable elements (JT) if scale == :linear - AE = ae + dae + AE .= ae + dae elseif scale == :log - AE = [ae[1] + dae[1], ae[2] + dae[2], ae[3] + dae[3], - ae[4] + dae[4], 10^(log10(ae[5]) + dae[5]), ae[6] + dae[6]] + AE .= [ae[1] + dae[1], ae[2] + dae[2], ae[3] + dae[3], + ae[4] + dae[4], 10^(log10(ae[5]) + dae[5]), ae[6] + dae[6]] end # Barycentric state vector q = attr2bary(A, AE, params) # Propagation and residuals # TO DO: `ρ::TaylorN` is too slow for `adam` due to evaluations # within the dynamical model - propres!(res, radec, jd0 - ae[5]/c_au_per_day, q, params; buffer, dynamics) + propres!(res, od, jd0 - ae[5]/c_au_per_day, q, params; buffer, idxs) iszero(length(res)) && break # Least squares fit fit = tryls(res, x0, 5, 1:4) + !fit.success && break x1 .= fit.x # Current Q Q = nms(res) @@ -177,11 +175,15 @@ function main() points_per_worker = [points[i] for i in idxs] # Evaluate `NEOs.adam` in each point result = pmap(points_per_worker) do points + # Orbit determination problem + od = ODProblem(newtonian!, radec) + # Pre-allocate output aes = Vector{Matrix{Float64}}(undef, length(points)) Qs = Vector{Vector{Float64}}(undef, length(points)) + # Main loop for i in eachindex(points) x, v_ρ = points[i] - aes[i], Qs[i] = adam(radec, A, 10^x, v_ρ, params; scale = :log, + aes[i], Qs[i] = adam(od, 1, A, 10^x, v_ρ, params; scale = :log, adamorder = varorder) j = findlast(!isinf, Qs[i]) if isnothing(j) diff --git a/scripts/mcmov.jl b/scripts/mcmov.jl index a0631476..a67d9b7f 100644 --- a/scripts/mcmov.jl +++ b/scripts/mcmov.jl @@ -42,17 +42,20 @@ end @everywhere begin using NEOs, Dates, TaylorSeries, PlanetaryEphemeris, JLD2 using NEOs: AdmissibleRegion, OpticalResidual, RadecMPC, PropagationBuffer, - reduce_tracklets, argoldensearch, scaled_variables, attr2bary, propres! + reduce_tracklets, argoldensearch, scaled_variables, attr2bary, + propres!, nobs function mcmov(points::Vector{Tuple{T, T}}, A::AdmissibleRegion{T}, radec::Vector{RadecMPC{T}}, params::NEOParameters{T}, varorder::Int = 2) where {T <: Real} + # Orbit determination problem + od = ODProblem(newtonian!, radec) # JT variables dae = scaled_variables("dx", ones(T, 6), order = varorder) - # Attributable elements + # Attributable elements (plain) ae = Vector{T}(undef, 6) - AE = Vector{TaylorN{T}}(undef, 6) ae[1:4] .= A.α, A.δ, A.v_α, A.v_δ + ae[5:6] .= points[1] # Scaling factors scalings = Vector{T}(undef, 6) scalings[1:4] .= abs.(ae[1:4]) ./ 1e6 @@ -60,13 +63,14 @@ end _, y_min = argoldensearch(A, 10^x_min, 10^x_max, :min, :outer, 1e-20) _, y_max = argoldensearch(A, 10^x_min, 10^x_max, :max, :outer, 1e-20) scalings[5:6] .= (x_max - x_min) / 100, (y_max - y_min) / 100 + # Attributable elements (jet transport) + AE = ae .+ scalings .* dae + # TDB epoch of admissible region + _jd0_ = dtutc2jdtdb(A.date) # Propagation buffer - t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end])) - tlim = (t0 - params.bwdoffset, tf + params.fwdoffset) - buffer = PropagationBuffer(newtonian!, dtutc2jdtdb(A.date), tlim, - ae .+ scalings .* dae, params) + buffer = PropagationBuffer(od, _jd0_, 1, nobs(od), AE, params) # Vector of residuals - res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec)) + res = [zero(OpticalResidual{T, TaylorN{T}}) for _ in eachindex(radec)] # Origin x0 = zeros(T, 6) # Manifold of variations @@ -82,9 +86,9 @@ end # Barycentric initial conditions (JT) q = attr2bary(A, AE, params) # Reference epoch in julian days (corrrected for light-time) - jd0 = dtutc2jdtdb(A.date) - constant_term(AE[5]) / c_au_per_day + jd0 = _jd0_ - constant_term(AE[5]) / c_au_per_day # Propagation and residuals - propres!(res, radec, jd0, q, params; buffer) + propres!(res, od, jd0, q, params; buffer) # Least squares fit fit = newtonls(res, x0, 10, 1:4) # Objective function diff --git a/scripts/orbitdetermination.jl b/scripts/orbitdetermination.jl index fa5705eb..69bdf8b0 100644 --- a/scripts/orbitdetermination.jl +++ b/scripts/orbitdetermination.jl @@ -78,8 +78,8 @@ end # Get at most 3 tracklets for orbit determination flag, radec = radecfilter(radec) !flag && return false - # Dynamical function - dynamics = newtonian! + # Orbit determination problem + od = ODProblem(newtonian!, radec) # Parameters params = NEOParameters( coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007, @@ -90,7 +90,7 @@ end # Start of computation init_time = now() # Initial orbit determination - sol = orbitdetermination(radec, params; dynamics, initcond) + sol = orbitdetermination(od, params; initcond) # Termination condition if length(sol.res) == length(radec) && nrms(sol) < 1.5 # Time of computation @@ -108,7 +108,7 @@ end jtlsiter = 20, newtoniter = 5 ) # Initial orbit determination - _sol_ = orbitdetermination(radec, params; dynamics, initcond) + _sol_ = orbitdetermination(od, params; initcond) # Time of computation Δ = (now() - init_time).value # Choose best orbit diff --git a/src/NEOs.jl b/src/NEOs.jl index 2026b2f2..cc0e6974 100644 --- a/src/NEOs.jl +++ b/src/NEOs.jl @@ -12,6 +12,7 @@ using Dates, InteractiveUtils, LazyArtifacts, LinearAlgebra, Printf, JSON, TaylorSeries, SatelliteToolboxTransformations, TaylorIntegration, SPICE, JLD2, Scratch using AutoHashEquals.Compat +using Base: RefValue using Dates: epochms2datetime using Downloads: download using DelimitedFiles: readdlm @@ -58,14 +59,14 @@ export loadpeeph, bwdfwdeph # Observer position in ECEF and ECI frames export obsposECEF, obsposvelECI # Optical astrometry processing -export compute_radec, select_debiasing_table, debiasing, w8sveres17, residuals, unfold, - relax_factor, outlier +export UniformWeights, Veres17, Farnocchia15, Eggl20, compute_radec, unfold, + residuals, outlier # Radar astrometry processing export compute_delay, radar_astrometry # Asteroid dynamical models export RNp1BP_pN_A_J23E_J2S_ng_eph_threads!, RNp1BP_pN_A_J23E_J2S_eph_threads!, newtonian! # Propagation -export NEOParameters, propagate, propagate_lyap, propagate_root +export propagate, propagate_lyap, propagate_root # Osculating export pv2kep, yarkp2adot # Least squares @@ -77,7 +78,7 @@ export tsaiod # Gauss method export gauss_method, gaussiod # Orbit determination -export curvature, issinglearc, orbitdetermination +export NEOParameters, ODProblem, curvature, issinglearc, orbitdetermination # B plane export valsecchi_circle, bopik, mtp # Outlier rejection diff --git a/src/observations/errormodel.jl b/src/observations/errormodel.jl new file mode 100644 index 00000000..3f95dbfc --- /dev/null +++ b/src/observations/errormodel.jl @@ -0,0 +1,421 @@ +@doc raw""" + AbstractErrorModel{T <: Real} + +Supertype for the optical astrometry error models API. +""" +abstract type AbstractErrorModel{T <: Real} end + +# Optical astrometry weighting schemes API +# All weighting schemes `W{T}` must: +# 1. be mutable subtypes of `AbstractWeightingScheme{T}`, +# 2. have a `w8s::Vector{T}` field, +# 3. implement a `W(::AbstractVector{RadecMPC{T}})` constructor, +# 4. override `getid(::W)`, +# 5. override `update!(::W{T}, ::AbstractVector{RadecMPC{T}})`. + +@doc raw""" + AbstractWeightingScheme{T} <: AbstractErrorModel{T} + +Supertype for the optical astrometry weighting schemes API. +""" +abstract type AbstractWeightingScheme{T} <: AbstractErrorModel{T} end + +# Print method for AbstractWeightingScheme +show(io::IO, w::AbstractWeightingScheme) = print(io, getid(w), + " weighting scheme with ", length(w.w8s), " observations") + +@doc raw""" + UniformWeights{T} <: AbstractWeightingScheme{T} + +Uniform optical astrometry weighting scheme. +""" +mutable struct UniformWeights{T} <: AbstractWeightingScheme{T} + w8s::Vector{T} + # Default constructor + UniformWeights(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} = + new{T}(ones(T, length(radec))) +end + +# Override getid +getid(w::UniformWeights) = "Uniform" + +# Override update! +function update!(w::UniformWeights{T}, + radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + w.w8s = ones(T, length(radec)) + return nothing +end + +@doc raw""" + Veres17{T} <: AbstractWeightingScheme{T} + +Veres et al. (2017) optical astrometry weighting scheme. + +!!! reference + See https://doi.org/10.1016/j.icarus.2017.05.021. +""" +mutable struct Veres17{T} <: AbstractWeightingScheme{T} + w8s::Vector{T} + # Default constructor + Veres17(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} = + new{T}(w8sveres17(radec)) +end + +# Override getid +getid(w::Veres17) = "Veres et al. (2017)" + +# Override update! +function update!(w::Veres17{T}, + radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + w.w8s = w8sveres17(radec) + return nothing +end + +@doc raw""" + σsveres17(obs::RadecMPC) + +Return the statistical uncertainty of `obs` according to Veres et al. (2017). + +!!! reference + See https://doi.org/10.1016/j.icarus.2017.05.021. +""" +function σsveres17(obs::RadecMPC) + + obscode = obs.observatory.code + dt_utc_obs = obs.date + catalogue = obs.catalogue.code + + # Unit weight (arcseconds) + w = 1.0 + # Table 2: epoch-dependent astrometric residuals + if obscode == "703" + return Date(dt_utc_obs) < Date(2014,1,1) ? w : 0.8w + elseif obscode ∈ ("691", "291") # Spacewatch, Kitt Peak, Arizona + return Date(dt_utc_obs) < Date(2003,1,1) ? 0.6w : 0.5w + elseif obscode == "644" + return Date(dt_utc_obs) < Date(2003,9,1) ? 0.6w : 0.4w + # Table 3: most active CCD asteroid observers + elseif obscode ∈ ("704", "C51", "J75") + return w + elseif obscode == "G96" + return 0.5w + elseif obscode ∈ ("F51", "F52") # Pan-STARRS 1 & 2, Haleakala, Hawaii + return 0.2w + elseif obscode ∈ ("G45", "608") + return 0.6w + elseif obscode == "699" + return 0.8w + elseif obscode ∈ ("D29", "E12") + return 0.75w + # Table 4: + elseif obscode ∈ ("645", "673", "H01") + return 0.3w + elseif obscode ∈ ("J04", "K92", "K93", "Q63", "Q64", "V37", "W85", "W86", "W87", + "K91", "E10", "F65") # Tenerife + Las Cumbres + return 0.4w + elseif obscode ∈ ("689", "950", "W84") + return 0.5w + # Applies only to program code assigned to M. Micheli + #elseif obscode ∈ ("G83", "309") + # if catalogue ∈ ("q", "t") # "q"=>"UCAC-4", "t"=>"PPMXL" + # return 0.3w + # elseif catalogue ∈ ("U", "V") # Gaia-DR1, Gaia-DR2 + # return 0.2w + # end + elseif obscode ∈ ("Y28",) + if catalogue ∈ ("t", "U", "V") + return 0.3w + else + return w + end + elseif obscode ∈ ("568",) + if catalogue ∈ ("o", "s") # "o"=>"USNO-B1.0", "s"=>"USNO-B2.0" + return 0.5w + elseif catalogue ∈ ("U", "V") # Gaia DR1, DR2 + return 0.1w + elseif catalogue ∈ ("t",) #"t"=>"PPMXL" + return 0.2w + else + return w + end + elseif obscode ∈ ("T09", "T12", "T14") && catalogue ∈ ("U", "V") # Gaia DR1, DR2 + return 0.1w + elseif catalogue == "" + return 1.5w + elseif catalogue != "" + return w + else + return w + end +end + +@doc raw""" + rexveres17(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + +Return the relax factor for each element of `radec` according to +Veres et al. (2017), which mitigates unresolved systematic errors +in observations taken on the same night by the same observatory. + +!!! reference + See https://doi.org/10.1016/j.icarus.2017.05.021. +""" +function rexveres17(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + # Convert to DataFrame + df = DataFrame(radec) + # Group by observatory and TimeOfDay + df.TimeOfDay = TimeOfDay.(radec) + gdf = groupby(df, [:observatory, :TimeOfDay]) + # Number of observations per tracklet + cdf = combine(gdf, nrow) + # Count observations in each group + Nv = cdf[gdf.groups, :nrow] + # Relaxation factor + return map(N -> N > 4.0 ? sqrt(N/4.0) : 1.0, Nv) +end + +@doc raw""" + w8sveres17(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + +Return the statistical weight of each element of `radec` according +to Veres et al. (2017). + +!!! reference + See https://doi.org/10.1016/j.icarus.2017.05.021. +""" +function w8sveres17(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + σs = σsveres17.(radec) + rex = rexveres17(radec) + return 1 ./ (rex .* σs) .^ 2 +end + +# Optical astrometry debiasing schemes API +# All debiasing schemes `D{T}` must: +# 1. be mutable subtypes of `AbstractDebiasingScheme{T}`, +# 2. have a `bias::Vector{Tuple{T, T}}` field, +# 3. implement a `D(::AbstractVector{RadecMPC{T}})` constructor, +# 4. override `getid(::D)`, +# 5. override `update!(::D{T}, ::AbstractVector{RadecMPC{T}})`. + +@doc raw""" + AbstractDebiasingScheme{T} <: AbstractErrorModel{T} + +Supertype for the optical astrometry debiasing schemes API. +""" +abstract type AbstractDebiasingScheme{T} <: AbstractErrorModel{T} end + +# Print method for AbstractDebiasingScheme +show(io::IO, d::AbstractDebiasingScheme) = print(io, getid(d), + " debiasing scheme with ", length(d.bias), " observations") + +@doc raw""" + Farnocchia15{T} <: AbstractDebiasingScheme{T} + +Farnocchia et al. (2015) optical astrometry debiasing scheme. + +!!! reference + See https://doi.org/10.1016/j.icarus.2014.07.033. +""" +mutable struct Farnocchia15{T} <: AbstractDebiasingScheme{T} + bias::Vector{Tuple{T, T}} + catcodes::Vector{String} + truth::String + resol::Resolution + table::Matrix{T} + # Default constructor + function Farnocchia15(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + catcodes, truth, resol, table = select_debiasing_table("2014") + bias = debiasing.(radec, Ref(catcodes), Ref(truth), + Ref(resol), Ref(table)) + return new{T}(bias, catcodes, truth, resol, table) + end +end + +# Override getid +getid(d::Farnocchia15) = "Farnocchia et al. (2015)" + +# Override update! +function update!(d::Farnocchia15{T}, + radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + d.bias = debiasing.(radec, Ref(d.catcodes), Ref(d.truth), + Ref(d.resol), Ref(d.table)) + return nothing +end + +@doc raw""" + Eggl20{T} <: AbstractDebiasingScheme{T} + +Eggl et al. (2020) optical astrometry debiasing scheme. + +!!! reference + See https://doi.org/10.1016/j.icarus.2019.113596. +""" +mutable struct Eggl20{T} <: AbstractDebiasingScheme{T} + bias::Vector{Tuple{T, T}} + hires::Bool + catcodes::Vector{String} + truth::String + resol::Resolution + table::Matrix{T} + # Default constructor + function Eggl20(radec::AbstractVector{RadecMPC{T}}, + hires::Bool = false) where {T <: Real} + id = hires ? "hires2018" : "2018" + catcodes, truth, resol, table = select_debiasing_table(id) + bias = debiasing.(radec, Ref(catcodes), Ref(truth), + Ref(resol), Ref(table)) + return new{T}(bias, hires, catcodes, truth, resol, table) + end +end + +# Override getid +getid(d::Eggl20) = string("Eggl et al. (2020)", d.hires ? + " [high resolution]" : "") + +# Override update! +function update!(d::Eggl20{T}, + radec::AbstractVector{RadecMPC{T}}) where {T <: Real} + d.bias = debiasing.(radec, Ref(d.catcodes), Ref(d.truth), + Ref(d.resol), Ref(d.table)) + return nothing +end + +@doc raw""" + select_debiasing_table(id::String) + +Return the catalogue codes, truth catalogue, resolution and bias matrix of the +corresponding debiasing scheme. The possible values for `id` are: +- `2014` corresponds to https://doi.org/10.1016/j.icarus.2014.07.033, +- `2018` corresponds to https://doi.org/10.1016/j.icarus.2019.113596 (default), +- `hires2018` corresponds to https://doi.org/10.1016/j.icarus.2019.113596 + (high resolution). +""" +function select_debiasing_table(id::String = "2018") + # Debiasing tables are loaded "lazily" via Julia artifacts, + # according to rules in Artifacts.toml + if id == "2018" + debias_path = artifact"debias_2018" + catcodes = mpc_catalogue_codes_2018 + # The healpix tesselation resolution of the bias map from + # https://doi.org/10.1016/j.icarus.2019.113596 + NSIDE = 64 + # In 2018 debias table Gaia DR2 catalogue is regarded as the truth + truth = "V" + elseif id == "hires2018" + debias_path = artifact"debias_hires2018" + catcodes = mpc_catalogue_codes_2018 + # The healpix tesselation resolution of the high-resolution bias map from + # https://doi.org/10.1016/j.icarus.2019.113596 + NSIDE = 256 + # In 2018 debias table Gaia DR2 catalogue is regarded as the truth + truth = "V" + elseif id == "2014" + debias_path = artifact"debias_2014" + catcodes = mpc_catalogue_codes_2014 + # The healpix tesselation resolution of the bias map from + # https://doi.org/10.1016/j.icarus.2014.07.033 + NSIDE = 64 + # In 2014 debias table PPMXL catalogue is regarded as the truth + truth = "t" + else + @error "Unknown bias map: $(id). Possible values are `2014`, `2018` \ + and `hires2018`." + end + + # Debias table file + bias_file = joinpath(debias_path, "bias.dat") + # Read bias matrix + table = readdlm(bias_file, comment_char = '!', comments = true) + # Initialize healpix Resolution variable + resol = Resolution(NSIDE) + # Compatibility between bias matrix and resolution + @assert size(table) == (resol.numOfPixels, 4length(catcodes)) "\ + Bias table file $bias_file dimensions do not match expected parameter \ + NSIDE = $NSIDE and/or number of catalogs in table." + + return catcodes, truth, resol, table +end + +@doc raw""" + debiasing(obs::RadecMPC{T}, catcodes::Vector{String}, truth::String, + resol::Resolution, table::Matrix{T}) where {T <: Real} + +Return total debiasing correction [arcsec] in both right ascension and declination. + +## Arguments + +- `obs::RadecMPC{T}`: optical observation. +- `catcodes::Vector{String}`: catalogues present in debiasing table. +- `truth::String`: truth catalogue of debiasing table. +- `resol::Resolution`: resolution. +- `table::Matrix{T}`: debiasing table. +""" +function debiasing(obs::RadecMPC{T}, catcodes::Vector{String}, + truth::String, resol::Resolution, table::Matrix{T}) where {T <: Real} + + # Catalogue code + catcode = obs.catalogue.code + + # If star catalogue is not present in debiasing table, + # then set corrections equal to zero + if (catcode ∉ catcodes) && catcode != "Y" + # Catalogue exists in CATALOGUES_MPC[] + if !isunknown(obs.catalogue) + # Truth catalogue is not present in debiasing table + # but it does not send a warning + if catcode != truth + @warn "Catalogue $(obs.catalogue.name) not found in debiasing table. \ + Setting debiasing corrections equal to zero." + end + # Unknown catalogue + elseif catcode == "" + @warn "Catalog information not available in observation record. \ + Setting debiasing corrections equal to zero." + # Catalogue code is not empty but it does not match an MPC catalogue code either + else + @warn "Catalog code $catcode does not correspond to MPC catalogue code. \ + Setting debiasing corrections equal to zero." + end + α_corr, δ_corr = zero(T), zero(T) + # If star catalogue is present in debiasing table, then compute corrections + else + # Get pixel tile index, assuming iso-latitude rings indexing, which is the + # formatting in `tiles.dat`. Substracting 1 from the returned value of + # `ang2pixRing` corresponds to 0-based indexing, as in `tiles.dat`; not + # substracting 1 from the returned value of `ang2pixRing` corresponds to + # 1-based indexing, as in Julia. Since we use pix_ind to get the corresponding + # row number in `bias.dat`, it's not necessary to substract 1. + pix_ind = ang2pixRing(resol, π/2 - obs.δ, obs.α) + + # Handle edge case: in new MPC catalogue nomenclature, "UCAC-5" -> "Y"; + # but in debias tables "UCAC-5" -> "W" + if catcode == "Y" + cat_ind = findfirst(x -> x == "W", catcodes) + else + cat_ind = findfirst(x -> x == catcode, catcodes) + end + + # Read dRA, pmRA, dDEC, pmDEC data from bias.dat + # dRA: position correction in RA * cos(DEC) at epoch J2000.0 [arcsec] + # dDEC: position correction in DEC at epoch J2000.0 [arcsec] + # pmRA: proper motion correction in RA*cos(DEC) [mas/yr] + # pmDEC: proper motion correction in DEC [mas/yr] + dRA, dDEC, pmRA, pmDEC = table[pix_ind, 4*cat_ind-3:4*cat_ind] + # Seconds since J2000 (TDB) + et_secs_i = dtutc2et(obs.date) + # Seconds sinde J2000 (TT) + tt_secs_i = et_secs_i - ttmtdb(et_secs_i/daysec) + # Years since J2000 + yrs_J2000_tt = tt_secs_i/(daysec*yr) + # Total debiasing correction in right ascension (arcsec) + α_corr = dRA + yrs_J2000_tt*pmRA/1_000 + # Total debiasing correction in declination (arcsec) + δ_corr = dDEC + yrs_J2000_tt*pmDEC/1_000 + end + + return α_corr, δ_corr +end + +debiasing(obs::RadecMPC{T}, catcodes::RefValue{Vector{String}}, + truth::RefValue{String}, resol::RefValue{Resolution}, + table::RefValue{Matrix{T}}) where {T <: Real} = + debiasing(obs, catcodes[], truth[], resol[], table[]) \ No newline at end of file diff --git a/src/observations/observations.jl b/src/observations/observations.jl index 04019928..b205cfdc 100644 --- a/src/observations/observations.jl +++ b/src/observations/observations.jl @@ -8,4 +8,5 @@ include("units.jl") include("jpl_eph.jl") include("topocentric.jl") include("tracklet.jl") +include("errormodel.jl") include("process_radec.jl") diff --git a/src/observations/process_radec.jl b/src/observations/process_radec.jl index 4ecd63bf..dd334c94 100644 --- a/src/observations/process_radec.jl +++ b/src/observations/process_radec.jl @@ -3,13 +3,12 @@ An astrometric optical observed minus computed residual. -# Fields +## Fields - `ξ_α::U`: right ascension residual. - `ξ_δ::U`: declination residual. - `w_α::T`: right ascension weight. - `w_δ::T`: declination weight. -- `relax_factor::T`: relaxation factor. - `outlier::Bool`: whether the residual is an outlier or not. """ @auto_hash_equals struct OpticalResidual{T <: Real, U <: Number} @@ -17,47 +16,52 @@ An astrometric optical observed minus computed residual. ξ_δ::U w_α::T w_δ::T - relax_factor::T outlier::Bool # Inner constructor - function OpticalResidual{T, U}(ξ_α::U, ξ_δ::U, w_α::T, w_δ::T, relax_factor::T = one(T), - outlier::Bool = false) where {T <: Real, U <: Number} - new{T, U}(ξ_α, ξ_δ, w_α, w_δ, relax_factor, outlier) + function OpticalResidual{T, U}(ξ_α::U, ξ_δ::U, w_α::T, w_δ::T, + outlier::Bool = false) where {T <: Real, U <: Number} + new{T, U}(ξ_α, ξ_δ, w_α, w_δ, outlier) end end + # Outer constructor -function OpticalResidual(ξ_α::U, ξ_δ::U, w_α::T, w_δ::T, relax_factor::T = one(T), - outlier::Bool = false) where {T <: Real, U <: Number} - return OpticalResidual{T, U}(ξ_α, ξ_δ, w_α, w_δ, relax_factor, outlier) -end +OpticalResidual(ξ_α::U, ξ_δ::U, w_α::T, w_δ::T, + outlier::Bool = false) where {T <: Real, U <: Number} = + OpticalResidual{T, U}(ξ_α, ξ_δ, w_α, w_δ, outlier) + +# Definition of zero OpticalResidual +zero(::Type{OpticalResidual{T, U}}) where {T <: Real, U <: Number} = + OpticalResidual{T, U}(zero(U), zero(U), zero(T), zero(T), false) # Evaluate methods function evaluate(res::OpticalResidual{T, TaylorN{T}}, x::Vector{T}) where {T <: Real} - return OpticalResidual(res.ξ_α(x), res.ξ_δ(x), res.w_α, res.w_δ, res.relax_factor, res.outlier) + OpticalResidual(res.ξ_α(x), res.ξ_δ(x), res.w_α, res.w_δ, res.outlier) end (res::OpticalResidual{T, TaylorN{T}})(x::Vector{T}) where {T <: Real} = evaluate(res, x) -function evaluate(res::AbstractVector{OpticalResidual{T, TaylorN{T}}}, x::Vector{T}) where {T <: Real} +function evaluate(res::AbstractVector{OpticalResidual{T, TaylorN{T}}}, + x::Vector{T}) where {T <: Real} res_new = Vector{OpticalResidual{T, T}}(undef, length(res)) for i in eachindex(res) res_new[i] = evaluate(res[i], x) end return res_new end -(res::Vector{OpticalResidual{T, TaylorN{T}}})(x::Vector{T}) where {T <: Real} = evaluate(res, x) +(res::Vector{OpticalResidual{T, TaylorN{T}}})(x::Vector{T}) where {T <: Real} = + evaluate(res, x) # Print method for OpticalResidual # Examples: # α: -138.79801 δ: -89.80025 # α: -134.79450 δ: -91.42509 (outlier) -function show(io::IO, x::OpticalResidual{T, U}) where {T <: Real, U <: Number} +function show(io::IO, x::OpticalResidual) outlier_flag = outlier(x) ? " (outlier)" : "" print(io, "α: ", @sprintf("%+.5f", cte(x.ξ_α)), " δ: ", @sprintf("%+.5f", cte(x.ξ_δ)), outlier_flag) end @doc raw""" - unfold(ξs::Vector{OpticalResidual{T, U}}) where {T <: Real, U <: Number} + unfold(ξs::AbstractVector{OpticalResidual{T, U}}) where {T <: Real, U <: Number} Concatenate right ascension and declination residuals for an orbit fit. """ @@ -75,10 +79,10 @@ function unfold(ξs::AbstractVector{OpticalResidual{T, U}}) where {T <: Real, U if !ξs[i].outlier # Right ascension res[k] = ξs[i].ξ_α - w[k] = ξs[i].w_α / ξs[i].relax_factor + w[k] = ξs[i].w_α # Declination res[k+L] = ξs[i].ξ_δ - w[k+L] = ξs[i].w_δ / ξs[i].relax_factor + w[k+L] = ξs[i].w_δ # Update global counter k += 1 end @@ -88,12 +92,11 @@ function unfold(ξs::AbstractVector{OpticalResidual{T, U}}) where {T <: Real, U end # Functions to get specific fields of a OpticalResidual object -ra(res::OpticalResidual{T, U}) where {T <: Real, U <: Number} = res.ξ_α -dec(res::OpticalResidual{T, U}) where {T <: Real, U <: Number} = res.ξ_δ -weight_ra(res::OpticalResidual{T, U}) where {T <: Real, U <: Number} = res.w_α -weight_dec(res::OpticalResidual{T, U}) where {T <: Real, U <: Number} = res.w_δ -relax_factor(res::OpticalResidual{T, U}) where {T <: Real, U <: Number} = res.relax_factor -outlier(res::OpticalResidual{T, U}) where {T <: Real, U <: Number} = res.outlier +ra(res::OpticalResidual) = res.ξ_α +dec(res::OpticalResidual) = res.ξ_δ +wra(res::OpticalResidual) = res.w_α +wdec(res::OpticalResidual) = res.w_δ +outlier(res::OpticalResidual) = res.outlier euclid3D(x::Vector{T}) where {T <: Real} = sqrt(dot3D(x, x)) function euclid3D(x::Vector{TaylorN{T}}) where {T <: Real} @@ -127,25 +130,25 @@ end dot3D(x::Vector{T}, y::Vector{TaylorN{T}}) where {T <: Real} = dot3D(y, x) @doc raw""" - compute_radec(observatory::ObservatoryMPC{T}, t_r_utc::DateTime; kwargs...) where {T <: AbstractFloat} - compute_radec(obs::RadecMPC{T}; kwargs...) where {T <: AbstractFloat} - compute_radec(obs::Vector{RadecMPC{T}}; kwargs...) where {T <: AbstractFloat} + compute_radec(observatory::ObservatoryMPC{T}, t_r_utc::DateTime; kwargs...) + compute_radec(obs::RadecMPC{T}; kwargs...) + compute_radec(obs::Vector{RadecMPC{T}}; kwargs...) where {T <: Real} -Compute astrometric right ascension and declination (both in arcsec) for a set of MPC-formatted -observations. Corrections due to Earth orientation, LOD, polar motion are considered in -computations. +Compute astrometric right ascension and declination [arcsec] for a set of MPC-formatted +observations. Corrections due to Earth orientation, LOD and polar motion are considered +in computations. -# Arguments +## Arguments - `observatory::ObservatoryMPC{T}`: observation site. - `t_r_utc::DateTime`: UTC time of astrometric observation. - `obs::Vector{RadecMPC{T}}`: optical observation(s). -# Keyword arguments +## Keyword arguments -- `niter::Int = 5`: number of light-time solution iterations. -- `xve::EarthEph = earthposvel`: Earth ephemeris. -- `xvs::SunEph = sunposvel`: Sun ephemeris. +- `niter::Int`: number of light-time solution iterations (default: `5`). +- `xve::EarthEph`: Earth ephemeris (default: `earthposvel`). +- `xvs::SunEph`: Sun ephemeris (default: `sunposvel`). - `xva::AstEph`: asteroid ephemeris. All ephemeris must take [et seconds since J2000] and return [barycentric position in km @@ -153,7 +156,7 @@ and velocity in km/sec]. """ function compute_radec(observatory::ObservatoryMPC{T}, t_r_utc::DateTime; niter::Int = 5, xvs::SunEph = sunposvel, xve::EarthEph = earthposvel, - xva::AstEph) where {T <: AbstractFloat, SunEph, EarthEph, AstEph} + xva::AstEph) where {T <: Real, SunEph, EarthEph, AstEph} # Transform receiving time from UTC to TDB seconds since J2000 et_r_secs = dtutc2et(t_r_utc) # Sun barycentric position and velocity at receive time @@ -282,11 +285,11 @@ function compute_radec(observatory::ObservatoryMPC{T}, t_r_utc::DateTime; niter: return α_as, δ_as # right ascension, declination both in arcsec end -function compute_radec(obs::RadecMPC{T}; kwargs...) where {T <: AbstractFloat} - return compute_radec(obs.observatory, obs.date; kwargs...) -end +compute_radec(obs::RadecMPC{T}; kwargs...) where {T <: Real} = + compute_radec(obs.observatory, obs.date; kwargs...) -function compute_radec(obs::Vector{RadecMPC{T}}; xva::AstEph, kwargs...) where {T <: AbstractFloat, AstEph} +function compute_radec(obs::Vector{RadecMPC{T}}; xva::AstEph, + kwargs...) where {T <: Real, AstEph} # Number of observations n_optical_obs = length(obs) @@ -312,225 +315,7 @@ function compute_radec(obs::Vector{RadecMPC{T}}; xva::AstEph, kwargs...) where { return vra, vdec # arcsec, arcsec end -@doc raw""" - select_debiasing_table(debias_table::String = "2018") - -Return the catalogue codes, truth catalogue, resolution and bias matrix of the corresponding -debias table. The possible values for `debias_table` are: -- `2014` corresponds to https://doi.org/10.1016/j.icarus.2014.07.033, -- `2018` corresponds to https://doi.org/10.1016/j.icarus.2019.113596, -- `hires2018` corresponds to https://doi.org/10.1016/j.icarus.2019.113596. -""" -function select_debiasing_table(debias_table::String = "2018") - # Debiasing tables are loaded "lazily" via Julia artifacts, according to rules in Artifacts.toml - if debias_table == "2018" - debias_path = artifact"debias_2018" - mpc_catalogue_codes_201X = mpc_catalogue_codes_2018 - # The healpix tesselation resolution of the bias map from https://doi.org/10.1016/j.icarus.2019.113596 - NSIDE = 64 - # In 2018 debias table Gaia DR2 catalogue is regarded as the truth - truth = "V" - elseif debias_table == "hires2018" - debias_path = artifact"debias_hires2018" - mpc_catalogue_codes_201X = mpc_catalogue_codes_2018 - # The healpix tesselation resolution of the high-resolution bias map from https://doi.org/10.1016/j.icarus.2019.113596 - NSIDE = 256 - # In 2018 debias table Gaia DR2 catalogue is regarded as the truth - truth = "V" - elseif debias_table == "2014" - debias_path = artifact"debias_2014" - mpc_catalogue_codes_201X = mpc_catalogue_codes_2014 - # The healpix tesselation resolution of the bias map from https://doi.org/10.1016/j.icarus.2014.07.033 - NSIDE= 64 - # In 2014 debias table PPMXL catalogue is regarded as the truth - truth = "t" - else - @error "Unknown bias map: $(debias_table). Possible values are `2014`, `2018` and `hires2018`." - end - - # Debias table file - bias_file = joinpath(debias_path, "bias.dat") - # Read bias matrix - bias_matrix = readdlm(bias_file, comment_char='!', comments=true) - # Initialize healpix Resolution variable - resol = Resolution(NSIDE) - # Compatibility between bias matrix and resolution - @assert size(bias_matrix) == (resol.numOfPixels, 4length(mpc_catalogue_codes_201X)) "Bias table file $bias_file dimensions do not match expected parameter NSIDE=$NSIDE and/or number of catalogs in table." - - return mpc_catalogue_codes_201X, truth, resol, bias_matrix -end - -@doc raw""" - debiasing(obs::RadecMPC{T}, mpc_catalogue_codes_201X::Vector{String}, truth::String, - resol::Resolution, bias_matrix::Matrix{T}) where {T <: AbstractFloat} - -Return total debiasing correction in right ascension and declination (both in arcsec). - -# Arguments - -- `obs::RadecMPC{T}`: optical observation. -- `mpc_catalogue_codes_201X::Vector{String}`: catalogues present in debiasing table. -- `truth::String`: truth catalogue of debiasing table. -- `resol::Resolution`: resolution. -- `bias_matrix::Matrix{T}`: debiasing table. -""" -function debiasing(obs::RadecMPC{T}, mpc_catalogue_codes_201X::Vector{String}, truth::String, - resol::Resolution, bias_matrix::Matrix{T}) where {T <: AbstractFloat} - - # Catalogue code - catcode = obs.catalogue.code - - # If star catalogue is not present in debiasing table, then set corrections equal to zero - if (catcode ∉ mpc_catalogue_codes_201X) && catcode != "Y" - # Catalogue exists in mpc_catalogues[] - if !isunknown(obs.catalogue) - # Truth catalogue is not present in debiasing table but it does not send a warning - if catcode != truth - @warn "Catalogue $(obs.catalogue.name) not found in debiasing table. Setting debiasing corrections equal to zero." - end - # Unknown catalogue - elseif catcode == "" - @warn "Catalog information not available in observation record. Setting debiasing corrections equal to zero." - # Catalogue code is not empty but it does not match an MPC catalogue code either - else - @warn "Catalog code $catcode does not correspond to MPC catalogue code. Setting debiasing corrections equal to zero." - end - α_corr = zero(T) - δ_corr = zero(T) - # If star catalogue is present in debiasing table, then compute corrections - else - # Get pixel tile index, assuming iso-latitude rings indexing, which is the formatting in `tiles.dat`. - # Substracting 1 from the returned value of `ang2pixRing` corresponds to 0-based indexing, as in `tiles.dat`; - # not substracting 1 from the returned value of `ang2pixRing` corresponds to 1-based indexing, as in Julia. - # Since we use pix_ind to get the corresponding row number in `bias.dat`, it's not necessary to substract 1. - pix_ind = ang2pixRing(resol, π/2 - obs.δ, obs.α) - - # Handle edge case: in new MPC catalogue nomenclature, "UCAC-5"->"Y"; but in debias tables "UCAC-5"->"W" - if catcode == "Y" - cat_ind = findfirst(x -> x == "W", mpc_catalogue_codes_201X) - else - cat_ind = findfirst(x -> x == catcode, mpc_catalogue_codes_201X) - end - - # Read dRA, pmRA, dDEC, pmDEC data from bias.dat - # dRA: position correction in RA * cos(DEC) at epoch J2000.0 [arcsec] - # dDEC: position correction in DEC at epoch J2000.0 [arcsec] - # pmRA: proper motion correction in RA*cos(DEC) [mas/yr] - # pmDEC: proper motion correction in DEC [mas/yr] - dRA, dDEC, pmRA, pmDEC = bias_matrix[pix_ind, 4*cat_ind-3:4*cat_ind] - # Seconds since J2000 (TDB) - et_secs_i = dtutc2et(obs.date) - # Seconds sinde J2000 (TT) - tt_secs_i = et_secs_i - ttmtdb(et_secs_i/daysec) - # Years since J2000 - yrs_J2000_tt = tt_secs_i/(daysec*yr) - # Total debiasing correction in right ascension (arcsec) - α_corr = dRA + yrs_J2000_tt*pmRA/1_000 - # Total debiasing correction in declination (arcsec) - δ_corr = dDEC + yrs_J2000_tt*pmDEC/1_000 - end - - return α_corr, δ_corr -end - -@doc raw""" - w8sveres17(obs::RadecMPC{T}) where {T <: AbstractFloat} - -Return the statistical weight from Veres et al. (2017) corresponding to `obs`. -""" -function w8sveres17(obs::RadecMPC{T}) where {T <: AbstractFloat} - - obscode = obs.observatory.code - dt_utc_obs = obs.date - catalogue = obs.catalogue.code - - # Unit weight (arcseconds) - w = 1.0 - # Table 2: epoch-dependent astrometric residuals - if obscode == "703" - return Date(dt_utc_obs) < Date(2014,1,1) ? w : 0.8w - elseif obscode ∈ ("691", "291") # Spacewatch, Kitt Peak, Arizona - return Date(dt_utc_obs) < Date(2003,1,1) ? 0.6w : 0.5w - elseif obscode == "644" - return Date(dt_utc_obs) < Date(2003,9,1) ? 0.6w : 0.4w - # Table 3: most active CCD asteroid observers - elseif obscode ∈ ("704", "C51", "J75") - return w - elseif obscode == "G96" - return 0.5w - elseif obscode ∈ ("F51", "F52") # Pan-STARRS 1 & 2, Haleakala, Hawaii - return 0.2w - elseif obscode ∈ ("G45", "608") - return 0.6w - elseif obscode == "699" - return 0.8w - elseif obscode ∈ ("D29", "E12") - return 0.75w - # Table 4: - elseif obscode ∈ ("645", "673", "H01") - return 0.3w - elseif obscode ∈ ("J04", "K92", "K93", "Q63", "Q64", "V37", "W85", "W86", "W87", - "K91", "E10", "F65") # Tenerife + Las Cumbres - return 0.4w - elseif obscode ∈ ("689", "950", "W84") - return 0.5w - #elseif obscode ∈ ("G83", "309") # Applies only to program code assigned to M. Micheli - # if catalogue ∈ ("q", "t") # "q"=>"UCAC-4", "t"=>"PPMXL" - # return 0.3w - # elseif catalogue ∈ ("U", "V") # Gaia-DR1, Gaia-DR2 - # return 0.2w - # end - elseif obscode ∈ ("Y28",) - if catalogue ∈ ("t", "U", "V") - return 0.3w - else - return w - end - elseif obscode ∈ ("568",) - if catalogue ∈ ("o", "s") # "o"=>"USNO-B1.0", "s"=>"USNO-B2.0" - return 0.5w - elseif catalogue ∈ ("U", "V") # Gaia DR1, DR2 - return 0.1w - elseif catalogue ∈ ("t",) #"t"=>"PPMXL" - return 0.2w - else - return w - end - elseif obscode ∈ ("T09", "T12", "T14") && catalogue ∈ ("U", "V") # Gaia DR1, DR2 - return 0.1w - elseif catalogue == "" - return 1.5w - elseif catalogue != "" - return w - else - return w - end -end - -@doc raw""" - relax_factor(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat} - -Return a relax factor for each element of `radec` quantifying the correlation between -observations taken on the same night by the same observatory. - -!!! reference - See https://doi.org/10.1016/j.icarus.2017.05.021. -""" -function relax_factor(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat} - # Convert to DataFrame - df = DataFrame(radec) - # Group by observatory and TimeOfDay - df.TimeOfDay = TimeOfDay.(radec) - gdf = groupby(df, [:observatory, :TimeOfDay]) - # Number of observations per tracklet - cdf = combine(gdf, nrow) - # Count observations in each group - Nv = cdf[gdf.groups, :nrow] - # Relaxation factor - return map(x -> x > 4.0 ? x/4.0 : 1.0, Nv) -end - -# Angle difference taking into account the discontinuity in [0, 2π] -> [0, 2π] +# Angle difference taking into account the discontinuity in [0, 2π) -> [0, 2π) # x and y must be in arcsec function anglediff(x::T, y::S) where {T, S <: Number} # Signed difference @@ -546,42 +331,37 @@ function anglediff(x::T, y::S) where {T, S <: Number} end @doc raw""" - residuals(radec::Vector{RadecMPC{T}}; kwargs...) where {T <: Real, AstEph} + residuals(radec::AbstractVector{RadecMPC{T}}, w8s::AbstractVector{T}, + bias::AbstractVector{Tuple{T, T}}; xva::AstEph, + kwargs...) where {AstEph, T <: Real} -Compute O-C residuals for optical astrometry. Corrections due to Earth orientation, LOD, -polar motion are computed by default. +Compute observed minus computed residuals for optical astrometry `radec`. +Corrections due to Earth orientation, LOD and polar motion are computed +by default. -See also [`compute_radec`](@ref), [`debiasing`](@ref), [`w8sveres17`](@ref) and -[`Healpix.ang2pixRing`](@ref). +See also [`OpticalResidual`](@ref) and [`compute_radec`](@ref). -# Arguments +## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of observations. +- `radec::AbstractVector{RadecMPC{T}}`: optical astrometry. +- `w8s::AbstractVector{T}`: statistical weights. +- `bias::AbstractVector{Tuple{T, T}}`: debiasing corrections. -# Keyword arguments +## Keyword arguments -- `niter::Int = 5`: number of light-time solution iterations. -- `debias_table::String = "2018"`: possible values are: - - `2014` corresponds to https://doi.org/10.1016/j.icarus.2014.07.033, - - `2018` corresponds to https://doi.org/10.1016/j.icarus.2019.113596, - - `hires2018` corresponds to https://doi.org/10.1016/j.icarus.2019.113596. -- `xve::EarthEph = earthposvel`: Earth ephemeris. -- `xvs::SunEph = sunposvel`: Sun ephemeris. +- `niter::Int`: number of light-time solution iterations (default: `5`). +- `xve::EarthEph`: Earth ephemeris (default: `earthposvel`). +- `xvs::SunEph`: Sun ephemeris (default: `sunposvel`). - `xva::AstEph`: asteroid ephemeris. All ephemeris must take [et seconds since J2000] and return [barycentric position in km and velocity in km/sec]. """ -function residuals(radec::Vector{RadecMPC{T}}; debias_table::String = "2018", xva::AstEph, - kwargs...) where {T <: Real, AstEph} - mpc_catalogue_codes_201X, truth, resol, bias_matrix = select_debiasing_table(debias_table) - return residuals(radec, mpc_catalogue_codes_201X, truth, resol, bias_matrix; xva, kwargs...) -end - -function residuals(radec::Vector{RadecMPC{T}}, mpc_catalogue_codes_201X::Vector{String}, - truth::String, resol::Resolution, bias_matrix::Matrix{T}; xva::AstEph, - kwargs...) where {T <: Real, AstEph} +function residuals(radec::AbstractVector{RadecMPC{T}}, w8s::AbstractVector{T}, + bias::AbstractVector{Tuple{T, T}}; xva::AstEph, kwargs...) where {AstEph, T <: Real} + # Check consistency between arrays + @assert length(radec) == length(w8s) == length(bias) # UTC time of first astrometric observation utc1 = date(radec[1]) # TDB seconds since J2000.0 for first astrometric observation @@ -592,38 +372,32 @@ function residuals(radec::Vector{RadecMPC{T}}, mpc_catalogue_codes_201X::Vector{ U = typeof(a1_et1) # Vector of residuals res = Vector{OpticalResidual{T, U}}(undef, length(radec)) - residuals!(res, radec, mpc_catalogue_codes_201X, truth, resol, - bias_matrix; xva, kwargs...) + residuals!(res, radec, w8s, bias; xva, kwargs...) return res end -function residuals!(res::Vector{OpticalResidual{T, U}}, radec::Vector{RadecMPC{T}}, - mpc_catalogue_codes_201X::Vector{String}, truth::String, resol::Resolution, - bias_matrix::Matrix{T}; xva::AstEph, kwargs...) where {T <: Real, U <: Number, AstEph} - - # Relax factors - rex = relax_factor(radec) +function residuals!(res::Vector{OpticalResidual{T, U}}, + radec::AbstractVector{RadecMPC{T}}, w8s::AbstractVector{T}, + bias::AbstractVector{Tuple{T, T}}; xva::AstEph, kwargs...) where {AstEph, + T <: Real, U <: Number} - tmap!(res, radec, rex) do obs, r + tmap!(res, radec, w8s, bias) do obs, w8, bias # Observed ra/dec α_obs = rad2arcsec(ra(obs)) # arcsec δ_obs = rad2arcsec(dec(obs)) # arcsec # Computed ra/dec α_comp, δ_comp = compute_radec(obs; xva = xva, kwargs...) # arcsec # Debiasing corrections - α_corr, δ_corr = debiasing(obs, mpc_catalogue_codes_201X, truth, resol, bias_matrix) - # Statistical weights from Veres et al. (2017) - w8 = w8sveres17(obs) - # O-C residual ra/dec + α_corr, δ_corr = bias + # Observed minus computed residual ra/dec # Note: ra is multiplied by a metric factor cos(dec) to match the format of # debiasing corrections return OpticalResidual{T, U}( anglediff(α_obs, α_comp) * cos(dec(obs)) - α_corr, δ_obs - δ_comp - δ_corr, - 1 / w8^2, - 1 / w8^2, - r, + w8, + w8, false ) end diff --git a/src/observations/tracklet.jl b/src/observations/tracklet.jl index b47635d1..31816164 100644 --- a/src/observations/tracklet.jl +++ b/src/observations/tracklet.jl @@ -39,14 +39,15 @@ A set of optical observations taken by the same observatory on the same night. end # Functions to get specific fields of a Tracklet object -date(x::Tracklet{T}) where {T <: AbstractFloat} = x.date -ra(x::Tracklet{T}) where {T <: AbstractFloat} = x.α -dec(x::Tracklet{T}) where {T <: AbstractFloat} = x.δ -observatory(x::Tracklet{T}) where {T <: AbstractFloat} = x.observatory -vra(x::Tracklet{T}) where {T <: AbstractFloat} = x.v_α -vdec(x::Tracklet{T}) where {T <: AbstractFloat} = x.v_δ -mag(x::Tracklet{T}) where {T <: AbstractFloat} = x.mag -indices(x::Tracklet{T}) where {T <: AbstractFloat} = x.indices +date(x::Tracklet) = x.date +ra(x::Tracklet) = x.α +dec(x::Tracklet) = x.δ +observatory(x::Tracklet) = x.observatory +vra(x::Tracklet) = x.v_α +vdec(x::Tracklet) = x.v_δ +mag(x::Tracklet) = x.mag +nobs(x::Tracklet) = x.nobs +indices(x::Tracklet) = x.indices # Print method for Tracklet{T} # Examples: @@ -192,7 +193,7 @@ Return: See section 5.1 of https://doi.org/10.1016/j.icarus.2007.11.033 """ function curvature(ts::AbstractVector{T}, αs::AbstractVector{T}, δs::AbstractVector{T}, - σs::AbstractVector{T}) where {T <: Real} + σs::AbstractVector{T}) where {T <: Real} @assert length(ts) == length(αs) == length(δs) == length(σs) ≥ 3 """ At least three observations needed for significant curvature computation.""" # Days of observation [relative to first observation] @@ -261,7 +262,7 @@ function curvature(radec::AbstractVector{RadecMPC{T}}) where {T <: Real} # Declination δs = dec.(radec) # Standard deviations [rad] - σs = arcsec2rad.(w8sveres17.(radec)) + σs = arcsec2rad.(σsveres17.(radec)) return curvature(ts, αs, δs, σs) end \ No newline at end of file diff --git a/src/orbitdetermination/gaussinitcond.jl b/src/orbitdetermination/gaussinitcond.jl index bed54f42..9308ecf4 100644 --- a/src/orbitdetermination/gaussinitcond.jl +++ b/src/orbitdetermination/gaussinitcond.jl @@ -5,7 +5,7 @@ A preliminary orbit obtained from Gauss method of orbit determination. See also [`gauss_method`](@ref). -# Fields +## Fields - `statevect::Vector{U}`: state vector at middle observation. - `ρ::Vector{U}`: topocentric ranges. @@ -32,19 +32,18 @@ See also [`gauss_method`](@ref). f_3::U g_3::U # Inner constructor - function GaussSolution{T, U}( - statevect::Vector{U}, ρ::Vector{U}, D::Matrix{U}, R_vec::Matrix{T}, - ρ_vec::Matrix{U}, τ_1::T, τ_3::T, f_1::U, g_1::U, f_3::U, g_3::U) where {T <: Real, U <: Number} + function GaussSolution{T, U}(statevect::Vector{U}, ρ::Vector{U}, D::Matrix{U}, + R_vec::Matrix{T}, ρ_vec::Matrix{U}, τ_1::T, τ_3::T, f_1::U, g_1::U, f_3::U, + g_3::U) where {T <: Real, U <: Number} return new{T, U}(statevect, ρ, D, R_vec, ρ_vec, τ_1, τ_3, f_1, g_1, f_3, g_3) end end # Outer constructor -function GaussSolution( - statevect::Vector{U}, ρ::Vector{U}, D::Matrix{U}, R_vec::Matrix{T}, - ρ_vec::Matrix{U}, τ_1::T, τ_3::T, f_1::U, g_1::U, f_3::U, g_3::U) where {T <: Real, U <: Number} - return GaussSolution{T, U}(statevect, ρ, D, R_vec, ρ_vec, τ_1, τ_3, f_1, g_1, f_3, g_3) -end +GaussSolution(statevect::Vector{U}, ρ::Vector{U}, D::Matrix{U}, + R_vec::Matrix{T}, ρ_vec::Matrix{U}, τ_1::T, τ_3::T, f_1::U, + g_1::U, f_3::U, g_3::U) where {T <: Real, U <: Number} = + GaussSolution{T, U}(statevect, ρ, D, R_vec, ρ_vec, τ_1, τ_3, f_1, g_1, f_3, g_3) # Print method for GaussSolution # Examples: @@ -110,7 +109,8 @@ function _format_Lagrange_equation(a::T, b::T, c::T) where {T <: Real} b_sgn = b ≥ 0 ? "+" : "-" c_sgn = c ≥ 0 ? "+" : "-" - return join(["r⁸ ", a_sgn, " ", abs(a), " r⁶ ", b_sgn, " ", abs(b), " r³ ", c_sgn, " ", abs(c), " = 0"]) + return join(["r⁸ ", a_sgn, " ", abs(a), " r⁶ ", b_sgn, " ", abs(b), " r³ ", c_sgn, + " ", abs(c), " = 0"]) end @doc raw""" @@ -190,11 +190,11 @@ end @doc raw""" gauss_method(obs::Vector{RadecMPC{T}}, params::NEOParameters{T}) where {T <: Real} gauss_method(observatories::Vector{ObservatoryMPC{T}}, dates::Vector{DateTime}, - α::Vector{U}, δ::Vector{U}, params::NEOParameters{T}) where {T <: Real, U <: Number} + α::Vector{U}, δ::Vector{U}, params::NEOParameters{T}) where {T <: Real, U <: Number} Core Gauss method of Initial Orbit determination (IOD). -# Arguments +## Arguments - `obs::Vector{RadecMPC{T}}`: three observations. - `observatories::Vector{ObservatoryMPC{T}}`: sites of observation. @@ -407,10 +407,12 @@ end # Update an initial condition `q`, obtained by Gauss' method, via ADAM # minimization over the middle tracklet's manifold of variations. -function _adam!(q::Vector{TaylorN{T}}, jd0::T, tracklet::Tracklet, - params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D} +function _adam!(od::ODProblem{D, T}, i::Int, q::Vector{TaylorN{T}}, jd0::T, + params::NEOParameters{T}) where {D, T <: Real} + # Middle tracklet + tracklet = od.tracklets[i] # Exploratory propagation - bwd, fwd, _ = propres(tracklet.radec, jd0, q(), params) + bwd, fwd, _ = propres(od, jd0, q(), params; idxs = indices(tracklet)) # Admissible region A = AdmissibleRegion(tracklet, params) # Epoch [days since J2000] @@ -426,41 +428,34 @@ function _adam!(q::Vector{TaylorN{T}}, jd0::T, tracklet::Tracklet, # Boundary projection ρ, v_ρ = boundary_projection(A, ρ, v_ρ) # ADAM - ae, _ = adam(tracklet.radec, A, ρ, v_ρ, params; scale = :log, dynamics = dynamics) + ae, _ = adam(od, i, A, ρ, v_ρ, params; scale = :log) # Epoch [julian days] (corrected for light-time) jd0 = dtutc2jdtdb(A.date) - ae[5] / c_au_per_day # Convert attributable elements to barycentric cartesian coordinates q0 = attr2bary(A, ae, params) # Jet Transport initial condition - q .= [q0[i] + (abs(q0[i]) / 10^5) * TaylorN(i, order = params.tsaorder) for i in 1:6] + q .= [q0[j] + (abs(q0[j]) / 10^5) * TaylorN(j, order = params.tsaorder) for j in 1:6] return jd0 end @doc raw""" - gaussiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, - params::NEOParameters{T}; kwargs...) where {T <: Real, D} + gaussiod(od::ODProblem{D, T}, params::NEOParameters{T}) where {D, T <: Real} -Fit a preliminary orbit to `radec` via jet transport Gauss method. +Fit a preliminary orbit to `od` via jet transport Gauss method. See also [`gauss_method`](@ref). ## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. -- `tracklets::Vector{Tracklet{T}},`: vector of tracklets. +- `od::ODProblem{D, T}`: an orbit determination problem. - `params::NEOParameters{T}`: see `Gauss' Method Parameters` of [`NEOParameters`](@ref). - -## Keyword arguments - -- `dynamics::D`: dynamical model (default: `newtonian!`). """ -function gaussiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, - params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D} +function gaussiod(od::ODProblem{D, T}, params::NEOParameters{T}) where {D, T <: Real} # Allocate memory for orbit sol = zero(NEOSolution{T, T}) # This function requires exactly 3 tracklets - length(tracklets) != 3 && return sol + length(od.tracklets) != 3 && return sol # Jet transport scaling factors scalings = fill(1e-6, 6) # Jet transport perturbation (ra/dec) @@ -471,7 +466,7 @@ function gaussiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, α = Vector{T}(undef, 3) δ = Vector{T}(undef, 3) # Find best triplet of observations - _gausstriplet!(observatories, dates, α, δ, tracklets) + _gausstriplet!(observatories, dates, α, δ, od.tracklets) # Julian day of middle observation _jd0_ = dtutc2jdtdb(dates[2]) # Gauss method solution @@ -487,18 +482,18 @@ function gaussiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, # Jet transport initial conditions q = solG[i].statevect .+ params.eph_su(jd0 - JD_J2000) # Jet Transport Least Squares - _sol_ = jtls(radec, tracklets, jd0, q, 0, params, true; dynamics) + _sol_ = jtls(od, jd0, q, 0, params, true) # Update solution - sol = updatesol(sol, _sol_, radec) + sol = updatesol(sol, _sol_, od.radec) # Termination condition nrms(sol) <= params.gaussQmax && return sol # ADAM help - tracklets[2].nobs < 2 && continue - jd0 = _adam!(q, jd0, tracklets[2], params; dynamics) + od.tracklets[2].nobs < 2 && continue + jd0 = _adam!(od, 2, q, jd0, params) # Jet Transport Least Squares - _sol_ = jtls(radec, tracklets, jd0, q, 0, params, true; dynamics) + _sol_ = jtls(od, jd0, q, 0, params, true) # Update solution - sol = updatesol(sol, _sol_, radec) + sol = updatesol(sol, _sol_, od.radec) # Termination condition nrms(sol) <= params.gaussQmax && return sol end diff --git a/src/orbitdetermination/neosolution.jl b/src/orbitdetermination/neosolution.jl index 1c80e5ad..26aec639 100644 --- a/src/orbitdetermination/neosolution.jl +++ b/src/orbitdetermination/neosolution.jl @@ -139,7 +139,7 @@ function zero(::Type{NEOSolution{T, U}}) where {T <: Real, U <: Number} x_fwd = Matrix{U}(undef, 0, 0) g_fwd = Vector{U}(undef, 0) res = Vector{OpticalResidual{T, U}}(undef, 0) - fit = LeastSquaresFit{T}(false, Vector{T}(undef, 0), Matrix{T}(undef, 0, 0), :newton) + fit = zero(LeastSquaresFit{T}) jacobian = Matrix{T}(undef, 0, 0) NEOSolution{T, U}( @@ -234,17 +234,16 @@ function jplcompare(des::String, sol::NEOSolution{T, U}) where {T <: Real, U <: end @doc raw""" - uncertaintyparameter(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, params::NEOParameters{T}; - dynamics::D = newtonian!) where {T <: Real, D} + uncertaintyparameter(od::ODProblem{D, T}, sol::NEOSolution{T, T}, + params::NEOParameters{T}) where {D, T <: Real} Return the Minor Planet Center Uncertainty Parameter. -# Arguments +## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. +- `od::ODProblem{D, T}`: an orbit determination problem. - `sol::NEOSolution{T, T}:` reference orbit. - `params::NEOParameters{T}`: see [`NEOParameters`](@ref). -- `dynamics::D`: dynamical model. !!! reference https://www.minorplanetcenter.net/iau/info/UValue.html @@ -252,10 +251,10 @@ Return the Minor Planet Center Uncertainty Parameter. !!! warning This function will change the (global) `TaylorSeries` variables. """ -function uncertaintyparameter(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, params::NEOParameters{T}; - dynamics::D = newtonian!) where {T <: Real, D} - # Reduce tracklets by polynomial regression - tracklets = reduce_tracklets(radec) +function uncertaintyparameter(od::ODProblem{D, T}, sol::NEOSolution{T, T}, + params::NEOParameters{T}) where {D, T <: Real} + # Check consistency between od and sol + @assert od.tracklets == sol.tracklets # Epoch [Julian days TDB] jd0 = sol.bwd.t0 + PE.J2000 # Barycentric initial conditions @@ -269,13 +268,13 @@ function uncertaintyparameter(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T} # Initial conditions q = q0 + dq # Propagation and residuals - bwd, fwd, res = propres(radec, jd0, q, params; dynamics) + bwd, fwd, res = propres(od, jd0, q, params) # Orbit fit fit = tryls(res, x0, params.newtoniter) # Residuals space to barycentric coordinates jacobian. J = Matrix(TS.jacobian(dq)) # Update solution - _sol_ = NEOSolution(tracklets, bwd, fwd, res, fit, J) + _sol_ = NEOSolution(od.tracklets, bwd, fwd, res, fit, J) # Osculating keplerian elements osc = pv2kep(_sol_() - params.eph_su(sol.bwd.t0); jd = jd0, frame = :ecliptic) # Eccentricity diff --git a/src/orbitdetermination/odproblem.jl b/src/orbitdetermination/odproblem.jl new file mode 100644 index 00000000..906dc008 --- /dev/null +++ b/src/orbitdetermination/odproblem.jl @@ -0,0 +1,91 @@ +@doc raw""" + ODProblem{D, T <: Real, WT <: AbstractWeightingScheme{T}, + DT <: AbstractDebiasingScheme{T}} + +An orbit determination problem. + +## Fields + +- `dynamics::D`: dynamical model. +- `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. +- `tracklets::Vector{Tracklet{T}}`: vector of tracklets. +- `w8s::WT`: weighting scheme. +- `bias::DT`: debiasing scheme. +""" +mutable struct ODProblem{D, T <: Real, WT <: AbstractWeightingScheme{T}, + DT <: AbstractDebiasingScheme{T}} + dynamics::D + radec::Vector{RadecMPC{T}} + tracklets::Vector{Tracklet{T}} + w8s::WT + bias::DT + @doc raw""" + ODProblem(dynamics::D, radec::Vector{RadecMPC{T}} [, + w8s::WT [, bias::DT]]) where {D, T <: Real, WT, DT} + + Return an orbit determination problem. + + ## Arguments + + - `dynamics::D`: dynamical model. + - `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. + - `w8s::WT`: weighting scheme (default: `Veres17`). + - `bias::DT`: debiasing scheme (default: `Eggl20`). + """ + function ODProblem(dynamics::D, radec::Vector{RadecMPC{T}}, + w8s::WT = Veres17, bias::DT = Eggl20) where {D, T <: Real, WT, DT} + # Reduce tracklets by polynomial regression + tracklets = reduce_tracklets(radec) + # Build weighting scheme + _w8s_ = w8s(radec) + # Build debiasing scheme + _bias_ = bias(radec) + # Consistency check + @assert length(radec) == sum(nobs, tracklets; init = 0) == + length(_w8s_.w8s) == length(_bias_.bias) + # Assemble orbit determination problem + new{D, T, typeof(_w8s_), typeof(_bias_)}(dynamics, radec, tracklets, + _w8s_, _bias_) + end +end + +const ODProblem{D, T} = ODProblem{D, T, WT, DT} where {D, T <: Real, + WT <: AbstractWeightingScheme{T}, DT <: AbstractDebiasingScheme{T}} + +# Print method for ODProblem +function show(io::IO, p::ODProblem) + t = " " + print(io, "Orbit determination problem:\n") + print(io, t, rpad("Dynamical model:", 21), p.dynamics, "\n") + print(io, t, rpad("Astrometry:", 21), length(p.radec), + " optical observations (", length(p.tracklets), " tracklets)\n") + print(io, t, rpad("Weighting scheme:", 21), getid(p.w8s), "\n") + print(io, t, rpad("Debiasing scheme:", 21), getid(p.bias), "\n") +end + +# Override update! +function update!(p::ODProblem{D, T}, radec::Vector{RadecMPC{T}}) where {D, T <: Real} + # Update ODProblem fields + p.radec = radec + p.tracklets = reduce_tracklets(radec) + update!(p.w8s, radec) + update!(p.bias, radec) + # Consistency check + @assert length(p.radec) == sum(nobs, p.tracklets; init = 0) == + length(p.w8s.w8s) == length(p.bias.bias) + return nothing +end + +# Override nobs +nobs(od::ODProblem) = length(od.radec) + +# Special PropagationBuffer constructor +function PropagationBuffer(od::ODProblem{D, T}, jd0::V, k0::Int, kf::Int, q0::Vector{U}, + params::NEOParameters{T}) where {D, T <: Real, U <: Number, V <: Number} + t0 = dtutc2days(date(od.radec[k0])) + tf = dtutc2days(date(od.radec[kf])) + tlim = (t0 - params.bwdoffset, tf + params.fwdoffset) + buffer = PropagationBuffer(od.dynamics, jd0, tlim, q0, params) + + return buffer +end \ No newline at end of file diff --git a/src/orbitdetermination/orbitdetermination.jl b/src/orbitdetermination/orbitdetermination.jl index f282f70b..fd07b2e6 100644 --- a/src/orbitdetermination/orbitdetermination.jl +++ b/src/orbitdetermination/orbitdetermination.jl @@ -1,10 +1,12 @@ include("osculating.jl") include("least_squares.jl") +include("odproblem.jl") include("neosolution.jl") include("admissibleregion.jl") # Times used within propres -function _proprestimes(radec::Vector{RadecMPC{T}}, jd0::U, params::NEOParameters{T}) where {T <: Real, U <: Number} +function _proprestimes(radec::AbstractVector{RadecMPC{T}}, jd0::U, + params::NEOParameters{T}) where {T <: Real, U <: Number} # Time of first (last) observation t0, tf = dtutc2jdtdb(date(radec[1])), dtutc2jdtdb(date(radec[end])) # TDB epoch (plain) @@ -17,26 +19,29 @@ function _proprestimes(radec::Vector{RadecMPC{T}}, jd0::U, params::NEOParameters end # Propagate an orbit and compute residuals -function propres(radec::Vector{RadecMPC{T}}, jd0::V, q0::Vector{U}, params::NEOParameters{T}; - buffer::Union{Nothing, PropagationBuffer{T, U, V}} = nothing, - dynamics::D = newtonian!) where {D, T <: Real, U <: Number, V <: Number} +function propres(od::ODProblem{D, T}, jd0::V, q0::Vector{U}, params::NEOParameters{T}; + buffer::Union{Nothing, PropagationBuffer{T, U, V}} = nothing, + idxs::AbstractVector{Int} = eachindex(od.radec)) where {D, T <: Real, U <: Number, + V <: Number} + # Subset of radec for propagation and residuals + radec = view(od.radec, idxs) # Times of first/last observation, epoch and years in backward/forward propagation t0, tf, _jd0_, nyears_bwd, nyears_fwd = _proprestimes(radec, jd0, params) # Propagation buffer if isnothing(buffer) - tlim = (t0 - JD_J2000 - params.bwdoffset, tf - JD_J2000 + params.fwdoffset) - buffer = PropagationBuffer(dynamics, jd0, tlim, q0, params) + buffer = PropagationBuffer(od, jd0, idxs[1], idxs[end], q0, params) end # Backward (forward) integration - bwd = _propagate(dynamics, jd0, nyears_bwd, q0, buffer, params) - fwd = _propagate(dynamics, jd0, nyears_fwd, q0, buffer, params) + bwd = _propagate(od.dynamics, jd0, nyears_bwd, q0, buffer, params) + fwd = _propagate(od.dynamics, jd0, nyears_fwd, q0, buffer, params) if !issuccessfulprop(bwd, t0 - _jd0_; tol = params.coeffstol) || !issuccessfulprop(fwd, tf - _jd0_; tol = params.coeffstol) return bwd, fwd, Vector{OpticalResidual{T, U}}(undef, 0) end # O-C residuals try - res = residuals(radec, params; + w8s, bias = view(od.w8s.w8s, idxs), view(od.bias.bias, idxs) + res = residuals(radec, w8s, bias; xvs = et -> auday2kmsec(params.eph_su(et/daysec)), xve = et -> auday2kmsec(params.eph_ea(et/daysec)), xva = et -> bwdfwdeph(et, bwd, fwd)) @@ -47,19 +52,22 @@ function propres(radec::Vector{RadecMPC{T}}, jd0::V, q0::Vector{U}, params::NEOP end # In-place method of propres -function propres!(res::Vector{OpticalResidual{T, U}}, radec::Vector{RadecMPC{T}}, jd0::V, q0::Vector{U}, - params::NEOParameters{T}; buffer::Union{Nothing, PropagationBuffer{T, U, V}} = nothing, - dynamics::D = newtonian!) where {D, T <: Real, U <: Number, V <: Number} +function propres!(res::Vector{OpticalResidual{T, U}}, od::ODProblem{D, T}, + jd0::V, q0::Vector{U}, params::NEOParameters{T}; + buffer::Union{Nothing, PropagationBuffer{T, U, V}} = nothing, + idxs::AbstractVector{Int} = eachindex(od.radec)) where {D, T <: Real, U <: Number, + V <: Number} + # Subset of radec for propagation and residuals + radec = view(od.radec, idxs) # Times of first/last observation, epoch and years in backward/forward propagation t0, tf, _jd0_, nyears_bwd, nyears_fwd = _proprestimes(radec, jd0, params) # Propagation buffer if isnothing(buffer) - tlim = (t0 - JD_J2000 - params.bwdoffset, tf - JD_J2000 + params.fwdoffset) - buffer = PropagationBuffer(dynamics, jd0, tlim, q0, params) + buffer = PropagationBuffer(od, jd0, idxs[1], idxs[end], q0, params) end # Backward (forward) integration - bwd = _propagate(dynamics, jd0, nyears_bwd, q0, buffer, params) - fwd = _propagate(dynamics, jd0, nyears_fwd, q0, buffer, params) + bwd = _propagate(od.dynamics, jd0, nyears_bwd, q0, buffer, params) + fwd = _propagate(od.dynamics, jd0, nyears_fwd, q0, buffer, params) if !issuccessfulprop(bwd, t0 - _jd0_; tol = params.coeffstol) || !issuccessfulprop(fwd, tf - _jd0_; tol = params.coeffstol) empty!(res) @@ -67,7 +75,8 @@ function propres!(res::Vector{OpticalResidual{T, U}}, radec::Vector{RadecMPC{T}} end # O-C residuals try - residuals!(res, radec, params; + w8s, bias = view(od.w8s.w8s, idxs), view(od.bias.bias, idxs) + residuals!(res, radec, w8s, bias; xvs = et -> auday2kmsec(params.eph_su(et/daysec)), xve = et -> auday2kmsec(params.eph_ea(et/daysec)), xva = et -> bwdfwdeph(et, bwd, fwd)) @@ -142,49 +151,39 @@ function addradec!(::Val{false}, rin::Vector{Int}, fit::LeastSquaresFit{T}, end @doc raw""" - jtls(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, jd0::V, - q::Vector{TayloN{T}}, i::Int, params::NEOParameters{T}, mode::Bool; - kwargs...) where {D, T <: Real, V <: Number} + jtls(od::ODProblem{D, T}, jd0::V, q::Vector{TayloN{T}}, i::Int, + params::NEOParameters{T} [, mode::Bool]) where {D, T <: Real, V <: Number} Compute an orbit via Jet Transport Least Squares. ## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. -- `tracklets::Vector{Tracklet{T}},`: vector of tracklets. +- `od::ODProblem{D, T}`: orbit determination problem. - `jd0::V`: reference epoch [Julian days TDB]. - `q::Vector{TaylorN{T}}`: jet transport initial condition. - `i::Int`: index of `tracklets` to start least squares fit. - `params::NEOParameters{T}`: see `Jet Transport Least Squares Parameters` of [`NEOParameters`](@ref). - `mode::Bool`: `addradec!` mode (default: `true`). - -## Keyword arguments - -- `dynamics::D`: dynamical model (default: `newtonian!`). """ -function jtls(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, - jd0::V, q::Vector{TaylorN{T}}, i::Int, params::NEOParameters{T}, - mode::Bool = true; dynamics::D = newtonian!) where {D, T <: Real, V <: Number} +function jtls(od::ODProblem{D, T}, jd0::V, q::Vector{TaylorN{T}}, i::Int, + params::NEOParameters{T}, mode::Bool = true) where {D, T <: Real, V <: Number} # Plain initial condition q0 = constant_term.(q) # JT tail dq = q - q0 # Vector of O-C residuals - res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec)) + res = [zero(OpticalResidual{T, TaylorN{T}}) for _ in eachindex(od.radec)] # Propagation buffer - t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end])) - tlim = (t0 - params.bwdoffset, tf + params.fwdoffset) - buffer = PropagationBuffer(dynamics, jd0, tlim, q, params) + buffer = PropagationBuffer(od, jd0, 1, nobs(od), q, params) # Origin x0 = zeros(T, 6) # Initial subset of radec for orbit fit - tin, tout, rin = _initialtracklets(tracklets, i) + tin, tout, rin = _initialtracklets(od.tracklets, i) # Residuals space to barycentric coordinates jacobian J = Matrix(TS.jacobian(dq)) # Best orbit - best_sol = zero(NEOSolution{T, T}) - best_Q = T(Inf) + best_sol, best_Q = zero(NEOSolution{T, T}), T(Inf) # Convergence flag flag = false # Jet transport least squares @@ -192,7 +191,7 @@ function jtls(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, # Initial conditions q = q0 + dq # Propagation & residuals - bwd, fwd = propres!(res, radec, jd0, q, params; buffer, dynamics) + bwd, fwd = propres!(res, od, jd0, q, params; buffer) iszero(length(res)) && break # Orbit fit fit = tryls(res[rin], x0, params.newtoniter) @@ -257,77 +256,68 @@ function issinglearc(radec::Vector{RadecMPC{T}}, arc::Day = Day(30)) where {T <: end @doc raw""" - orbitdetermination(radec::Vector{RadecMPC{T}}, params::NEOParameters{T}; - kwargs...) where {T <: Real, D1, D2} + orbitdetermination(od::ODProblem{D, T}, params::NEOParameters{T}; + kwargs...) where {D, I, T <: Real} Initial Orbit Determination (IOD) routine. ## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. +- `od::ODProblem{D, T}`: an orbit determination problem. - `params::NEOParameters{T}`: see [`NEOParameters`](@ref). ## Keyword arguments -- `dynamics::D1`: dynamical model (default: `newtonian!`). -- `initcond::D2`: naive initial conditions function; takes as input an - `AdmissibleRegion` and outputs a `Vector{Tuple{T, T, Symbol}}`, +- `initcond::I`: naive initial conditions function; takes as input an + `AdmissibleRegion{T}` and outputs a `Vector{Tuple{T, T, Symbol}}`, where each element has the form `(ρ, v_ρ, scale)` (default: `iodinitcond`). !!! warning This function will change the (global) `TaylorSeries` variables. """ -function orbitdetermination(radec::Vector{RadecMPC{T}}, params::NEOParameters{T}; - dynamics::D1 = newtonian!, initcond::D2 = iodinitcond) where {T <: Real, D1, D2} +function orbitdetermination(od::ODProblem{D, T}, params::NEOParameters{T}; + initcond::I = iodinitcond) where {D, I, T <: Real} # Allocate memory for orbit sol = zero(NEOSolution{T, T}) - # Eliminate observatories without coordinates - filter!(x -> hascoord(observatory(x)), radec) + # Cannot handle observatories without coordinates + all(x -> hascoord(observatory(x)), od.radec) || return sol # Cannot handle zero observations or multiple arcs - (isempty(radec) || !issinglearc(radec)) && return sol - # Reduce tracklets by polynomial regression - tracklets = reduce_tracklets(radec) + (isempty(od.radec) || !issinglearc(od.radec)) && return sol # Set jet transport variables varorder = max(params.tsaorder, params.gaussorder, params.jtlsorder) scaled_variables("dx", ones(T, 6); order = varorder) # Gauss method - _sol_ = gaussiod(radec, tracklets, params; dynamics) + _sol_ = gaussiod(od, params) # Update solution - sol = updatesol(sol, _sol_, radec) + sol = updatesol(sol, _sol_, od.radec) # Termination condition nrms(sol) <= params.gaussQmax && return sol # Too short arc - _sol_ = tsaiod(radec, tracklets, params; dynamics, initcond) + _sol_ = tsaiod(od, params; initcond) # Update solution - sol = updatesol(sol, _sol_, radec) + sol = updatesol(sol, _sol_, od.radec) return sol end @doc raw""" - orbitdetermination(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, - params::NEOParameters{T}; kwargs...) where {T <: Real, D} + orbitdetermination(od::ODProblem{D, T}, sol::NEOSolution{T, T}, + params::NEOParameters{T}) where {D, T <: Real} -Fit a least squares orbit to `radec` using `sol` as an initial condition. +Fit a least squares orbit to `od` using `sol` as an initial condition. ## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. +- `od::ODProblem{D, T}`: orbit determination problem. - `sol::NEOSolution{T, T}:` preliminary orbit. - `params::NEOParameters{T}`: see [`NEOParameters`](@ref). -## Keyword arguments - -- `dynamics::D`: dynamical model (default: `newtonian!`). - !!! warning This function will change the (global) `TaylorSeries` variables. """ -function orbitdetermination(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, - params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D} - # Reduce tracklets by polynomial regression - tracklets = reduce_tracklets(radec) +function orbitdetermination(od::ODProblem{D, T}, sol::NEOSolution{T, T}, + params::NEOParameters{T}) where {D, T <: Real} # Reference epoch [Julian days TDB] jd0 = sol.bwd.t0 + PE.J2000 # Plain barycentric initial condition @@ -339,6 +329,6 @@ function orbitdetermination(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, # Jet Transport initial condition q = q0 + dq # Jet Transport Least Squares - _, i = findmin(tracklet -> abs(dtutc2days(tracklet.date) - sol.bwd.t0), tracklets) - return jtls(radec, tracklets, jd0, q, i, params; dynamics) + _, i = findmin(t -> abs(dtutc2days(t.date) - sol.bwd.t0), od.tracklets) + return jtls(od, jd0, q, i, params) end \ No newline at end of file diff --git a/src/orbitdetermination/tooshortarc.jl b/src/orbitdetermination/tooshortarc.jl index 5dd140b8..cd444a3f 100644 --- a/src/orbitdetermination/tooshortarc.jl +++ b/src/orbitdetermination/tooshortarc.jl @@ -1,10 +1,10 @@ @doc raw""" - adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, - params::NEOParameters{T}; kwargs...) where {T <: Real, D} + adam(od::ODProblem{D, T}, i::Int, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, + params::NEOParameters{T}; kwargs...) where {D, T <: Real} Adaptative moment estimation (ADAM) minimizer of normalized mean square -residual of `radec` over the manifold of variations of `A`, starting -from `(ρ, v_ρ)`. +residual of `od.tracklets[i].radec` over the manifold of variations of `A`, +starting from `(ρ, v_ρ)`. ## Keyword arguments @@ -14,26 +14,22 @@ from `(ρ, v_ρ)`. - `ν::T`: second moment (default: `0.9`). - `ϵ::T`: numerical stability constant (default: `1e-8`). - `adamorder::Int`: jet transport order (default: `2`). -- `dynamics::D`: dynamical model (default: `newtonian!`). !!! reference See Algorithm 1 of https://doi.org/10.48550/arXiv.1412.6980. """ -function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, - params::NEOParameters{T}; scale::Symbol = :linear, η::T = 25.0, - μ::T = 0.75, ν::T = 0.9, ϵ::T = 1e-8, adamorder::Int = 2, - dynamics::D = newtonian!) where {T <: Real, D} - # Initial time of integration [julian days] +function adam(od::ODProblem{D, T}, i::Int, A::AdmissibleRegion{T}, ρ::T, v_ρ::T, + params::NEOParameters{T}; scale::Symbol = :linear, η::T = 25.0, + μ::T = 0.75, ν::T = 0.9, ϵ::T = 1e-8, adamorder::Int = 2) where {D, T <: Real} + # Initial time of integration [julian days TDB] jd0 = dtutc2jdtdb(A.date) - # Maximum number of iterations - maxiter = params.adamiter - # Target function relative tolerance - Qtol = params.adamQtol + # Unfold maximum number of iterations and relative tolerance + maxiter, Qtol = params.adamiter, params.adamQtol # Allocate memory aes = Matrix{T}(undef, 6, maxiter+1) Qs = fill(T(Inf), maxiter+1) # Initial attributable elements - aes[:, 1] .= [A.α, A.δ, A.v_α, A.v_δ, ρ, v_ρ] + aes[:, 1] .= A.α, A.δ, A.v_α, A.v_δ, ρ, v_ρ # Scaling factors scalings = Vector{T}(undef, 6) scalings[1:4] .= abs.(aes[1:4, 1]) ./ 1e6 @@ -43,34 +39,32 @@ function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T scalings[5] = (log10(A.ρ_domain[2]) - log10(A.ρ_domain[1])) / 1_000 end scalings[6] = (A.v_ρ_domain[2] - A.v_ρ_domain[1]) / 1_000 - # Jet transport variables + # Jet transport variables and initial condition dae = [scalings[i] * TaylorN(i, order = adamorder) for i in 1:6] + AE = aes[:, 1] .+ dae + # Subset of radec + idxs = indices(od.tracklets[i]) # Propagation buffer - t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end])) - tlim = (t0 - params.bwdoffset, tf + params.fwdoffset) - buffer = PropagationBuffer(dynamics, jd0, tlim, aes[:, 1] .+ dae, params) + buffer = PropagationBuffer(od, jd0, idxs[1], idxs[end], AE, params) # Vector of O-C residuals - res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec)) + res = [zero(OpticalResidual{T, TaylorN{T}}) for _ in eachindex(idxs)] # Origin - x0 = zeros(T, 6) - x1 = zeros(T, 6) + x0, x1 = zeros(T, 6), zeros(T, 6) # Gradient of objective function wrt (ρ, v_ρ) g_t = Vector{T}(undef, 2) # First momentum - m = zeros(T, 2) - _m_ = zeros(T, 2) + m, _m_ = zeros(T, 2), zeros(T, 2) # Second momentum - n = zeros(T, 2) - _n_ = zeros(T, 2) + n, _n_ = zeros(T, 2), zeros(T, 2) # Gradient descent for t in 1:maxiter # Current attributable elements (plain) ae = aes[:, t] # Attributable elements (JT) if scale == :linear - AE = ae + dae + AE .= ae + dae elseif scale == :log - AE = [ae[1] + dae[1], ae[2] + dae[2], ae[3] + dae[3], + AE .= [ae[1] + dae[1], ae[2] + dae[2], ae[3] + dae[3], ae[4] + dae[4], 10^(log10(ae[5]) + dae[5]), ae[6] + dae[6]] end # Barycentric state vector @@ -78,10 +72,11 @@ function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T # Propagation and residuals # TO DO: `ρ::TaylorN` is too slow for `adam` due to evaluations # within the dynamical model - propres!(res, radec, jd0 - ae[5]/c_au_per_day, q, params; buffer, dynamics) + propres!(res, od, jd0 - ae[5]/c_au_per_day, q, params; buffer, idxs) iszero(length(res)) && break # Least squares fit fit = tryls(res, x0, 5, 1:4) + !fit.success && break x1 .= fit.x # Current Q Q = nms(res) @@ -112,47 +107,41 @@ function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T end @doc raw""" - tsaiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, - params::NEOParameters{T}; kwargs...) where {T <: Real, D1, D2} + tsaiod(od::ODProblem{D, T}, params::NEOParameters{T}; + initcond::I = iodinitcond) where {D, I, T <: Real} - - -Compute an orbit by minimizing the normalized root mean square residual -over the manifold of variations. +Fit a preliminary orbit to `od` via jet transport minimization +of the normalized mean square residual over the manifold of variations. ## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of optical astrometry. -- `tracklets::Vector{Tracklet{T}},`: vector of tracklets. +- `od::ODProblem{D, T}`: an orbit determination problem. - `params::NEOParameters{T}`: see `Too Short Arc Parameters` of [`NEOParameters`](@ref). ## Keyword arguments -- `dynamics::D1`: dynamical model (default: `newtonian!`). -- `initcond::D2`: naive initial conditions function; takes as input an - `AdmissibleRegion` and outputs a `Vector{Tuple{T, T, Symbol}}`, +- `initcond::I`: naive initial conditions function; takes as input an + `AdmissibleRegion{T}` and outputs a `Vector{Tuple{T, T, Symbol}}`, where each element has the form `(ρ, v_ρ, scale)` (default: `iodinitcond`). """ -# Too short arc initial orbit determination -function tsaiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, - params::NEOParameters{T}; dynamics::D1 = newtonian!, - initcond::D2 = iodinitcond) where {T <: Real, D1, D2} +function tsaiod(od::ODProblem{D, T}, params::NEOParameters{T}; + initcond::I = iodinitcond) where {D, I, T <: Real} # Allocate memory for orbit sol = zero(NEOSolution{T, T}) # Iterate tracklets - for i in eachindex(tracklets) + for i in eachindex(od.tracklets) # ADAM requires a minimum of 2 observations - tracklets[i].nobs < 2 && continue + od.tracklets[i].nobs < 2 && continue # Admissible region - A = AdmissibleRegion(tracklets[i], params) + A = AdmissibleRegion(od.tracklets[i], params) # List of naive initial conditions I0 = initcond(A) # Iterate naive initial conditions for j in eachindex(I0) # ADAM minimization over manifold of variations ρ, v_ρ, scale = I0[j] - ae, Q = adam(tracklets[i].radec, A, ρ, v_ρ, params; scale, dynamics) + ae, Q = adam(od, i, A, ρ, v_ρ, params; scale) # ADAM failed to converge isinf(Q) && continue # Initial time of integration [julian days] @@ -165,9 +154,9 @@ function tsaiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, # Jet Transport initial condition q = [q0[k] + scalings[k] * TaylorN(k, order = params.tsaorder) for k in 1:6] # Jet Transport Least Squares - _sol_ = jtls(radec, tracklets, jd0, q, i, params, false; dynamics) + _sol_ = jtls(od, jd0, q, i, params, false) # Update solution - sol = updatesol(sol, _sol_, radec) + sol = updatesol(sol, _sol_, od.radec) # Termination condition nrms(sol) <= params.tsaQmax && return sol end diff --git a/src/postprocessing/outlier_rejection.jl b/src/postprocessing/outlier_rejection.jl index 37be8e91..187a29e2 100644 --- a/src/postprocessing/outlier_rejection.jl +++ b/src/postprocessing/outlier_rejection.jl @@ -5,26 +5,30 @@ include("bplane.jl") Return the contribution of `x` to the nrms. """ -residual_norm(x::OpticalResidual{T, T}) where {T <: Real} = x.w_α * x.ξ_α^2 / x.relax_factor + x.w_δ * x.ξ_δ^2 / x.relax_factor +residual_norm(x::OpticalResidual{T, T}) where {T <: Real} = + x.w_α * x.ξ_α^2 + x.w_δ * x.ξ_δ^2 @doc raw""" - outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, - params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D} + outlier_rejection(od::ODProblem{D, T}, sol::NEOSolution{T, T}, + params::NEOParameters{T}) where {D, T <: Real} -Refine an orbit, computed by [`orbitdetermination`](@ref), via propagation and/or outlier rejection. +Refine an orbit, computed by [`orbitdetermination`](@ref), +via propagation and/or outlier rejection. -# Arguments +## Arguments -- `radec::Vector{RadecMPC{T}}`: vector of observations. +- `od::ODProblem{D, T}`: an orbit determination problem. - `sol::NEOSolution{T, T}`: orbit to be refined. -- `params::NEOParameters{T}`: see `Outlier Rejection Parameters` of [`NEOParameters`](@ref). -- `dynamics::D`: dynamical model. +- `params::NEOParameters{T}`: see `Outlier Rejection Parameters` + of [`NEOParameters`](@ref). !!! warning This function will change the (global) `TaylorSeries` variables. """ -function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, - params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D} +function outlier_rejection(od::ODProblem{D, T}, sol::NEOSolution{T, T}, + params::NEOParameters{T}) where {D, T <: Real} + # Check consistency between od and sol + @assert od.tracklets == sol.tracklets # Julian day (TDB) to start propagation jd0 = sol.bwd.t0 + PE.J2000 # Initial conditions (T) @@ -36,7 +40,7 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, # Initial conditions (jet transport) q = q0 .+ dq # Propagation and residuals - bwd, fwd, res = propres(radec, jd0, q, params; dynamics) + bwd, fwd, res = propres(od, jd0, q, params) iszero(length(res)) && return zero(NEOSolution{T, T}) # Origin x0 = zeros(T, 6) @@ -52,14 +56,14 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, end # Number of observations - N_radec = length(radec) + N_radec = length(od.radec) # Maximum allowed outliers max_drop = ceil(Int, N_radec * params.max_per / 100) # Boolean mask (0: included in fit, 1: outlier) new_outliers = BitVector(zeros(Int, N_radec)) # Drop loop - for i in 1:max_drop + for _ in 1:max_drop # Contribution of each residual to nrms norms = residual_norm.(res(fit.x)) # Iterate norms from largest to smallest @@ -70,8 +74,8 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, # Drop residual new_outliers[j] = true # Update residuals - res = OpticalResidual.(ra.(res), dec.(res), weight_ra.(res), weight_dec.(res), - relax_factor.(res), new_outliers) + res = OpticalResidual.(ra.(res), dec.(res), wra.(res), wdec.(res), + new_outliers) # Update fit fit = tryls(res, x0, params.newtoniter) break @@ -103,8 +107,8 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, # Recover residual new_outliers[j] = false # Update residuals - res = OpticalResidual.(ra.(res), dec.(res), weight_ra.(res), weight_dec.(res), - relax_factor.(res), new_outliers) + res = OpticalResidual.(ra.(res), dec.(res), wra.(res), wdec.(res), + new_outliers) # Update fit fit = tryls(res, x0, params.newtoniter) end @@ -117,8 +121,8 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, # Reset boolean mask new_outliers[1:end] .= false # Update residuals - res = OpticalResidual.(ra.(res), dec.(res), weight_ra.(res), weight_dec.(res), - relax_factor.(res), new_outliers) + res = OpticalResidual.(ra.(res), dec.(res), wra.(res), wdec.(res), + new_outliers) # Update fit fit = tryls(res, x0, params.newtoniter) @@ -147,8 +151,8 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, # Outliers new_outliers[idxs] .= true # Update residuals - res = OpticalResidual.(ra.(res), dec.(res), weight_ra.(res), weight_dec.(res), - relax_factor.(res), new_outliers) + res = OpticalResidual.(ra.(res), dec.(res), wra.(res), wdec.(res), + new_outliers) # Update fit fit = tryls(res, x0, params.newtoniter) diff --git a/src/propagation/parameters.jl b/src/propagation/parameters.jl index 4665f4f9..55f1d961 100644 --- a/src/propagation/parameters.jl +++ b/src/propagation/parameters.jl @@ -1,92 +1,32 @@ -struct NEOParameters{T <: Real} - # Propagation Parameters - maxsteps::Int - μ_ast::Vector{T} - order::Int - abstol::T - parse_eqs::Bool - bwdoffset::T - fwdoffset::T - coeffstol::T - # Residuals Parameters - mpc_catalogue_codes_201X::Vector{String} - truth::String - resol::Resolution - bias_matrix::Matrix{T} - eph_su::TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}} - eph_ea::TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}} - # Gauss' Method Parameters - max_triplets::Int - gaussorder::Int - adamhelp::Bool - gaussQmax::T - # Too Short Arc Parameters - H_max::T - a_max::T - adamiter::Int - adamQtol::T - tsaorder::Int - tsaQmax::T - # Jet Transport Least Squares Parameters - jtlsiter::Int - newtoniter::Int - jtlsorder::Int - # Outlier Rejection Parameters - max_per::T - # Inner constructor is generated by default -end - -# Print method for NEOParameters -function show(io::IO, p::NEOParameters{T}) where {T <: Real} - params = [ - :maxsteps, :order, :abstol, :parse_eqs, :bwdoffset, - :fwdoffset, :coeffstol, :max_triplets, :gaussorder, - :adamhelp, :gaussQmax, :H_max, :a_max, :adamiter, - :adamQtol, :tsaorder, :tsaQmax, :jtlsiter, :newtoniter, - :jtlsorder, :max_per - ] - s = Vector{String}(undef, length(params)+1) - s[1] = "NEOParameters{$T}:\n" - for i in eachindex(params) - x = string(params[i], ":") - s[i+1] = string(" ", rpad(x, 15), getfield(p, params[i]), "\n") - end - - print(io, join(s)) -end - -# Outer constructors - @doc raw""" NEOParameters([params::NEOParameters{T};] kwargs...) where {T <: Real} Parameters for all orbit determination functions. -# Propagation Parameters +## Propagation Parameters - `maxsteps::Int`: maximum number of steps for the integration (default: `500`). -- `μ_ast::Vector{T}`: vector of gravitational parameters (default: `μ_ast343_DE430[1:end]`). +- `μ_ast::Vector{T}`: vector of gravitational parameters + (default: `μ_ast343_DE430[1:end]`). - `order::Int`: order of Taylor expansions wrt time (default: 25). -- `abstol::T`: absolute tolerance used to compute propagation timestep (default: `1e-20`). -- `parse_eqs::Bool`: whether to use the specialized method of `jetcoeffs` or not (default: `true`). -- `bwdoffset/fwdoffset::T`: days to propagate beyond first (bwd) / last (fwd) observation (default: `0.5`). +- `abstol::T`: absolute tolerance used to compute propagation timestep + (default: `1e-20`). +- `parse_eqs::Bool`: whether to use the specialized method of `jetcoeffs` or not + (default: `true`). +- `bwdoffset/fwdoffset::T`: days to propagate beyond first (bwd) / last (fwd) observation + (default: `0.5`). - `coeffstol::T`: maximum size of the coefficients (default: `10.0`). -# Residuals Parameters - -- `debias_table::String`: debiasing scheme (default: `"2018"`). Possible values are: - - `"2014"` corresponds to https://doi.org/10.1016/j.icarus.2014.07.033, - - `"2018"` corresponds to https://doi.org/10.1016/j.icarus.2019.113596, - - `"hires2018"` corresponds to https://doi.org/10.1016/j.icarus.2019.113596. +## Gauss Method Parameters -# Gauss' Method Parameters - -- `max_triplets::Int`: maximum number of triplets to check for a solution (default: `10`). +- `max_triplets::Int`: maximum number of triplets to check for a solution + (default: `10`). - `gaussorder::Int`: order of the jet transport perturbation (default: `5`). -- `adamhelp::Bool`: whether to refine Gauss' preliminary orbit via ADAM (default: `false`). +- `adamhelp::Bool`: whether to refine Gauss preliminary orbit via ADAM + (default: `false`). - `gaussQmax::T`: nrms threshold (default: `5.0`). -# Too Short Arc Parameters +## Too Short Arc Parameters - `H_max::T`: maximum absolute magnitude (default: `34.5`). - `a_max::T`: maximum semimajor axis (default: `100.0`). @@ -95,43 +35,53 @@ Parameters for all orbit determination functions. - `tsaorder::Int`: order of the jet transport perturbation (default: `6`). - `tsaQmax::T`: nrms threshold (default: `1.5`). -# Jet Transport Least Squares Parameters +## Jet Transport Least Squares Parameters -- `jtlsiter::Int`: maximum number of iterations for jet transport least squares (default: `5`). -- `newtoniter::Int`: number of iterations for differential corrections / Newton's method (default: `5`). +- `jtlsiter::Int`: maximum number of iterations for jet transport least squares + (default: `5`). +- `newtoniter::Int`: number of iterations for differential corrections / Newton's method + (default: `5`). - `jtlsorder::Int`: order of the jet transport perturbation (default: `5`). -# Outlier Rejection Parameters +## Outlier Rejection Parameters - `max_per::T`: maximum allowed rejection percentage (default: `18.0`). """ -function NEOParameters(; - maxsteps::Int = 500, μ_ast::Vector{T} = μ_ast343_DE430[1:end], order::Int = 25, - abstol::T = 1e-20, parse_eqs::Bool = true, bwdoffset::T = 0.5, fwdoffset::T = 0.5, - coeffstol::T = 10.0, debias_table::String = "2018", max_triplets::Int = 10, - gaussorder::Int = 5, adamhelp::Bool = false, gaussQmax::T = 5.0, H_max::T = 34.5, - a_max::T = 100.0, adamiter::Int = 200, adamQtol::T = 0.001, tsaorder::Int = 6, - tsaQmax::T = 1.5, jtlsiter::Int = 5, newtoniter::Int = 5, jtlsorder::Int = 5, - max_per::T = 18.0) where {T <: Real} - - # Unfold debiasing matrix - mpc_catalogue_codes_201X, truth, resol, bias_matrix = select_debiasing_table(debias_table) - # Sun (Earth) ephemeris - _su = selecteph(sseph, su) - eph_su = TaylorInterpolant(_su.t0, _su.t, collect(_su.x)) - _ea = selecteph(sseph, ea) - eph_ea = TaylorInterpolant(_ea.t0, _ea.t, collect(_ea.x)) - # Assemble NEOParameters - return NEOParameters{T}( - maxsteps, μ_ast, order, abstol, parse_eqs, bwdoffset, fwdoffset, - coeffstol, mpc_catalogue_codes_201X, truth, resol, bias_matrix, - eph_su, eph_ea, max_triplets, gaussorder, adamhelp, gaussQmax, H_max, - a_max, adamiter, adamQtol, tsaorder, tsaQmax, jtlsiter, newtoniter, - jtlsorder, max_per - ) +@kwdef struct NEOParameters{T <: Real} + # Propagation Parameters + maxsteps::Int = 500 + μ_ast::Vector{T} = μ_ast343_DE430[1:end] + order::Int = 25 + abstol::T = 1e-20 + parse_eqs::Bool = true + bwdoffset::T = 0.5 + fwdoffset::T = 0.5 + coeffstol::T = 10.0 + # Sun (earth) ephemeris + eph_su::TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}} = _loadephsu() + eph_ea::TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}} = _loadephea() + # Gauss' Method Parameters + max_triplets::Int = 10 + gaussorder::Int = 5 + adamhelp::Bool = false + gaussQmax::T = 5.0 + # Too Short Arc Parameters + H_max::T = 34.5 + a_max::T = 100.0 + adamiter::Int = 200 + adamQtol::T = 0.001 + tsaorder::Int = 6 + tsaQmax::T = 1.5 + # Jet Transport Least Squares Parameters + jtlsiter::Int = 5 + newtoniter::Int = 5 + jtlsorder::Int = 5 + # Outlier Rejection Parameters + max_per::T = 18.0 end -function NEOParameters(params::NEOParameters{T}; kwargs...) where {T <: AbstractFloat} +# Outer constructors +function NEOParameters(params::NEOParameters{T}; kwargs...) where {T <: Real} fields = fieldnames(NEOParameters{T}) vals = Vector{Any}(undef, length(fields)) for i in eachindex(vals) @@ -145,13 +95,31 @@ function NEOParameters(params::NEOParameters{T}; kwargs...) where {T <: Abstract return NEOParameters{T}(vals...) end -function residuals(radec::Vector{RadecMPC{T}}, params::NEOParameters{T}; kwargs...) where {T <: Real} - return residuals(radec, params.mpc_catalogue_codes_201X, params.truth, params.resol, - params.bias_matrix; kwargs...) +# Print method for NEOParameters +function show(io::IO, p::NEOParameters{T}) where {T <: Real} + avoid = [:μ_ast, :eph_su, :eph_ea] + params = fieldnames(NEOParameters{T}) + s = Vector{String}(undef, length(params)+1) + s[1] = "NEOParameters{$T}:\n" + for i in eachindex(params) + if params[i] in avoid + s[i+1] = "" + else + x = string(params[i], ":") + s[i+1] = string(" ", rpad(x, 15), getfield(p, params[i]), "\n") + end + end + + print(io, join(s)) end -function residuals!(res::Vector{OpticalResidual{T, U}}, radec::Vector{RadecMPC{T}}, - params::NEOParameters{T}; kwargs...) where {T <: Real, U <: Number} - return residuals!(res, radec, params.mpc_catalogue_codes_201X, params.truth, params.resol, - params.bias_matrix; kwargs...) +# Load Sun (Earth) ephemeris +function _loadephsu() + _su = selecteph(sseph, su) + return TaylorInterpolant(_su.t0, _su.t, collect(_su.x)) +end + +function _loadephea() + _ea = selecteph(sseph, ea) + return TaylorInterpolant(_ea.t0, _ea.t, collect(_ea.x)) end \ No newline at end of file diff --git a/test/bplane.jl b/test/bplane.jl index 4c1f7e54..dfb778bb 100644 --- a/test/bplane.jl +++ b/test/bplane.jl @@ -14,13 +14,16 @@ using NEOs: propres radec = fetch_radec_mpc("2018 LA") # Parameters params = NEOParameters(coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007) + # Orbit determination problem + od = ODProblem(newtonian!, radec[1:8]) # Initial Orbit Determination - sol = orbitdetermination(radec[1:8], params) + sol = orbitdetermination(od, params) # Update parameters params = NEOParameters(params; jtlsorder = 6) # Refine orbit - sol = orbitdetermination(radec, sol, params) + NEOs.update!(od, radec) + sol = orbitdetermination(od, sol, params) # Radial velocity with respect to the Earth. function rvelea(t, fwd, params) @@ -34,7 +37,7 @@ using NEOs: propres # Time of close approach params = NEOParameters(params; fwdoffset = 0.3) - bwd, fwd, res = propres(radec, epoch(sol) + J2000, sol(), params) + bwd, fwd, res = propres(od, epoch(sol) + J2000, sol(), params) t_CA = find_zeros(t -> rvelea(t, fwd, params), fwd.t0, fwd.t0 + fwd.t[end-1])[1] # Asteroid's geocentric state vector xae = fwd(t_CA) - params.eph_ea(t_CA) diff --git a/test/orbitdetermination.jl b/test/orbitdetermination.jl index f74e95da..4e5fe49a 100644 --- a/test/orbitdetermination.jl +++ b/test/orbitdetermination.jl @@ -29,11 +29,13 @@ end # Parameters params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007, parse_eqs = false) params = NEOParameters(params, parse_eqs = true) + # Orbit determination problem + od = ODProblem(newtonian!, subradec) # Initial Orbit Determination - sol = orbitdetermination(subradec, params) + sol = orbitdetermination(od, params) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Orbit solution @test isa(sol, NEOSolution{Float64, Float64}) @@ -74,7 +76,8 @@ end @test numberofdays(subradec) < 2.76 # Refine orbit - sol1 = orbitdetermination(subradec, sol, params) + NEOs.update!(od, subradec) + sol1 = orbitdetermination(od, sol, params) # Orbit solution @test isa(sol1, NEOSolution{Float64, Float64}) @@ -106,7 +109,7 @@ end # Compatibility with JPL @test all(abs.(sol1() - JPL) ./ sigmas(sol1) .< 0.31) # MPC Uncertainty Parameter - @test uncertaintyparameter(radec, sol1, params) == 6 + @test uncertaintyparameter(od, sol1, params) == 7 end @testset "Gauss Method (with ADAM)" begin @@ -121,11 +124,13 @@ end # Parameters params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007, adamhelp = true) + # Orbit determination problem + od = ODProblem(newtonian!, radec) # Initial Orbit Determination - sol = orbitdetermination(radec, params) + sol = orbitdetermination(od, params) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Orbit solution @test isa(sol, NEOSolution{Float64, Float64}) @@ -159,7 +164,7 @@ end -0.000263645106358707, 0.01837375504784015, 0.007208464525828968] @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 4.8e-2) # MPC Uncertainty Parameter - @test uncertaintyparameter(radec, sol, params) == 8 + @test uncertaintyparameter(od, sol, params) == 8 end @testset "Admissible region" begin @@ -177,7 +182,7 @@ end # Admissible region A = AdmissibleRegion(tracklet, params) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Zero AdmissibleRegion @test iszero(zero(AdmissibleRegion{Float64})) @@ -281,11 +286,13 @@ end # Parameters params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007) + # Orbit determination problem + od = ODProblem(newtonian!, radec) # Initial Orbit Determination - sol = orbitdetermination(radec, params) + sol = orbitdetermination(od, params) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Curvature C, Γ_C = curvature(radec) @@ -325,7 +332,7 @@ end -0.009512301266159554, -0.01532548565855646, -0.00809464581680694] @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.012) # MPC Uncertainty Parameter - @test uncertaintyparameter(radec, sol, params) == 11 + @test uncertaintyparameter(od, sol, params) == 11 end @testset "Outlier Rejection" begin @@ -339,11 +346,13 @@ end # Parameters params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007) + # Orbit determination problem + od = ODProblem(newtonian!, subradec) # Initial Orbit Determination - sol = orbitdetermination(subradec, params) + sol = orbitdetermination(od, params) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Orbit solution @test isa(sol, NEOSolution{Float64, Float64}) @@ -378,9 +387,10 @@ end @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 1.75) # Add remaining observations - sol1 = orbitdetermination(radec, sol, params) + NEOs.update!(od, radec) + sol1 = orbitdetermination(od, sol, params) # Outlier rejection - sol1 = outlier_rejection(radec, sol1, params) + sol1 = outlier_rejection(od, sol1, params) # Orbit solution @test isa(sol1, NEOSolution{Float64, Float64}) @@ -412,7 +422,7 @@ end # Compatibility with JPL @test all(abs.(sol1() - JPL) ./ sigmas(sol1) .< 1.2e-3) # MPC Uncertainty Parameter - @test uncertaintyparameter(radec, sol1, params) == 8 + @test uncertaintyparameter(od, sol1, params) == 8 end @testset "Interesting NEOs" begin @@ -430,11 +440,13 @@ end # Parameters params = NEOParameters(coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007) + # Orbit determination problem + od = ODProblem(newtonian!, radec) - # Orbit Determination - sol = orbitdetermination(radec, params) + # Initial Orbit Determination + sol = orbitdetermination(od, params) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Curvature C, Γ_C = curvature(radec) @@ -474,7 +486,7 @@ end -0.017557851117612377, -0.005781634223099801, -0.0020075106081869185] @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.3) # MPC Uncertainty Parameter - @test uncertaintyparameter(radec, sol, params) == 10 + @test uncertaintyparameter(od, sol, params) == 10 # 2008 TC3 entered the Earth's atmosphere around October 7, 2008, 02:46 UTC @@ -488,11 +500,13 @@ end # Parameters params = NEOParameters(coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007) + # Orbit determination problem + od = ODProblem(newtonian!, subradec) # Initial Orbit Determination - sol = orbitdetermination(subradec, params) + sol = orbitdetermination(od, params) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Orbit solution @test isa(sol, NEOSolution{Float64, Float64}) @@ -533,7 +547,8 @@ end @test numberofdays(subradec) < 0.70 # Refine orbit - sol1 = orbitdetermination(subradec, sol, params) + NEOs.update!(od, subradec) + sol1 = orbitdetermination(od, sol, params) # Orbit solution @test isa(sol1, NEOSolution{Float64, Float64}) @@ -569,7 +584,7 @@ end # TODO: understand better differences wrt JPL solutions # @test nrms(sol1) < nrms(sol) # MPC Uncertainty Parameter - @test uncertaintyparameter(subradec, sol1, params) == 5 + @test uncertaintyparameter(od, sol1, params) == 5 end @testset "scripts/orbitdetermination.jl" begin @@ -612,11 +627,13 @@ end adamiter = 500, adamQtol = 1e-5, tsaQmax = 2.0, jtlsiter = 20, newtoniter = 10 ) + # Orbit determination problem + od = ODProblem(newtonian!, radec) # Initial Orbit Determination - sol = orbitdetermination(radec, params; initcond = iodinitcond) + sol = orbitdetermination(od, params; initcond = iodinitcond) - # Values by Sep 25, 2024 + # Values by Oct 1, 2024 # Orbit solution @test isa(sol, NEOSolution{Float64, Float64}) diff --git a/test/propagation.jl b/test/propagation.jl index 36334b03..78eb3147 100644 --- a/test/propagation.jl +++ b/test/propagation.jl @@ -19,6 +19,8 @@ using InteractiveUtils: methodswith @test !isempty(methodswith(Val{RNp1BP_pN_A_J23E_J2S_eph_threads!}, TaylorIntegration.jetcoeffs!)) @test !isempty(methodswith(Val{RNp1BP_pN_A_J23E_J2S_eph_threads!}, TaylorIntegration._allocate_jetcoeffs!)) + @test !isempty(methodswith(Val{newtonian!}, TaylorIntegration.jetcoeffs!)) + @test !isempty(methodswith(Val{newtonian!}, TaylorIntegration._allocate_jetcoeffs!)) end using PlanetaryEphemeris: selecteph, ea, su, daysec, auday2kmsec @@ -93,13 +95,15 @@ using InteractiveUtils: methodswith @test norm(sol_bwd(sol_bwd.t0 + sol_bwd.t[end])-q_bwd_end, Inf) < 1e-12 # Read optical astrometry file - obs_radec_mpc_2023DW = read_radec_mpc(joinpath(pkgdir(NEOs), "test", "data", "RADEC_2023_DW.dat")) + # Make weigths and debiasing corrections + w8s = Veres17(obs_radec_mpc_2023DW).w8s + bias = Eggl20(obs_radec_mpc_2023DW).bias # Compute residuals _res_ = NEOs.residuals( obs_radec_mpc_2023DW, - params, + w8s, bias, xve = t -> auday2kmsec(eph_ea(t/daysec)), xvs = t -> auday2kmsec(eph_su(t/daysec)), xva = t -> auday2kmsec(sol(t/daysec)) @@ -132,7 +136,7 @@ using InteractiveUtils: methodswith # Compute residuals for orbit with perturbed initial conditions _res1_ = NEOs.residuals( obs_radec_mpc_2023DW, - params, + w8s, bias, xve = t -> auday2kmsec(eph_ea(t/daysec)), xvs = t -> auday2kmsec(eph_su(t/daysec)), xva = t -> auday2kmsec(sol1(t/daysec)) @@ -197,11 +201,14 @@ using InteractiveUtils: methodswith # Read optical astrometry file obs_radec_mpc_apophis = read_radec_mpc(joinpath(pkgdir(NEOs), "test", "data", "99942_Tholen_etal_2013.dat")) + # Make weights and debiasing corrections + w8s = Veres17(obs_radec_mpc_apophis).w8s + bias = Eggl20(obs_radec_mpc_apophis).bias # Compute optical astrometry residuals _res_radec_ = NEOs.residuals( obs_radec_mpc_apophis, - params, + w8s, bias, xve = t -> auday2kmsec(eph_ea(t/daysec)), xvs = t -> auday2kmsec(eph_su(t/daysec)), xva = t -> auday2kmsec(sol(t/daysec)) @@ -372,9 +379,12 @@ using InteractiveUtils: methodswith # Read optical astrometry file obs_radec_mpc_apophis = read_radec_mpc(joinpath(pkgdir(NEOs), "test", "data", "99942_Tholen_etal_2013.dat")) + # Make weights and debiasing corrections + w8s = Veres17(obs_radec_mpc_apophis).w8s + bias = Eggl20(obs_radec_mpc_apophis).bias # Compute optical astrometry residuals - _res_radec_ = NEOs.residuals(obs_radec_mpc_apophis, params; xvs, xve, xva) + _res_radec_ = NEOs.residuals(obs_radec_mpc_apophis, w8s, bias; xvs, xve, xva) res_radec, w_radec = NEOs.unfold(_res_radec_) nobsopt = round(Int, length(res_radec))