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))