diff --git a/Project.toml b/Project.toml
index 3c50d3cf..beb4d61b 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.0"
+version = "0.9.1"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
diff --git a/examples/recipes.jl b/examples/recipes.jl
index 7996d5db..f82cebf6 100644
--- a/examples/recipes.jl
+++ b/examples/recipes.jl
@@ -3,7 +3,7 @@
 using NEOs, Plots
 
 # Download optical astrometry of asteroid 2023 DW
-radec = fetch_radec_mpc("designation" => "2023 DW")
+radec = fetch_radec_mpc("2023 DW")
 
 # Orbit determination
 params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007)
diff --git a/pha/apophis.jl b/pha/apophis.jl
index a03b6b61..a34d67f3 100644
--- a/pha/apophis.jl
+++ b/pha/apophis.jl
@@ -41,7 +41,7 @@ function parse_commandline()
 
     @add_arg_table! s begin
         "--jd0"
-            help = "Initial epoch calendar date (TDB)"
+            help = "Initial epoch calendar date (UTC)"
             arg_type = DateTime
             default = DateTime(2020, 12, 17)
         "--varorder"
@@ -121,7 +121,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     q0 = vcat(q00, 0.0, 0.0) .+ dq
 
     # Initial date (in Julian days)
-    jd0 = datetime2julian(jd0_datetime)
+    jd0 = dtutc2jdtdb(jd0_datetime)
 
     print_header("Integrator warmup", 2)
     params = NEOParameters(;maxsteps=1, order, abstol, parse_eqs)
@@ -132,7 +132,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     print_header("Main integration", 2)
     tmax = nyears_bwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
-    println("• Final time of integration: ", julian2datetime(jd0 + tmax))
+    println("• Final time of integration: ", jdtdb2utc(jd0 + tmax))
 
     params = NEOParameters(;maxsteps, order, abstol, parse_eqs)
     sol_bwd = NEOs.propagate(dynamics, jd0, nyears_bwd, q0, params)
@@ -141,7 +141,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
 
     tmax = nyears_fwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
-    println("• Final time of integration: ", julian2datetime(jd0 + tmax))
+    println("• Final time of integration: ", jdtdb2utc(jd0 + tmax))
 
     sol_fwd, tvS, xvS, gvS = NEOs.propagate_root(dynamics, jd0, nyears_fwd, q0, params)
     jldsave("Apophis_fwd.jld2"; sol_fwd, tvS, xvS, gvS)
diff --git a/scripts/adam.jl b/scripts/adam.jl
index b7c2b9db..4d45dff1 100644
--- a/scripts/adam.jl
+++ b/scripts/adam.jl
@@ -50,7 +50,7 @@ end
         μ::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]
-        jd0 = datetime2julian(A.date)
+        jd0 = dtutc2jdtdb(A.date)
         # Maximum number of iterations
         maxiter = params.adamiter
         # Allocate memory
@@ -71,7 +71,7 @@ end
         set_variables(Float64, "dx"; order = adamorder, numvars = 6)
         dae = [scalings[i] * TaylorN(i, order = adamorder) for i in 1:6]
         # Propagation buffer
-        t0, tf = datetime2days(date(radec[1])), datetime2days(date(radec[end]))
+        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)
         # Vector of O-C residuals
diff --git a/scripts/mcmov.jl b/scripts/mcmov.jl
index 284c7681..a0631476 100644
--- a/scripts/mcmov.jl
+++ b/scripts/mcmov.jl
@@ -61,9 +61,9 @@ end
         _, 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
         # Propagation buffer
-        t0, tf = datetime2days(date(radec[1])), datetime2days(date(radec[end]))
+        t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end]))
         tlim = (t0 - params.bwdoffset, tf + params.fwdoffset)
-        buffer = PropagationBuffer(newtonian!, datetime2julian(A.date), tlim,
+        buffer = PropagationBuffer(newtonian!, dtutc2jdtdb(A.date), tlim,
                                    ae .+ scalings .* dae, params)
         # Vector of residuals
         res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec))
@@ -82,7 +82,7 @@ end
             # Barycentric initial conditions (JT)
             q = attr2bary(A, AE, params)
             # Reference epoch in julian days (corrrected for light-time)
-            jd0 = datetime2julian(A.date) - constant_term(AE[5]) / c_au_per_day
+            jd0 = dtutc2jdtdb(A.date) - constant_term(AE[5]) / c_au_per_day
             # Propagation and residuals
             propres!(res, radec, jd0, q, params; buffer)
             # Least squares fit
diff --git a/scripts/names.jl b/scripts/names.jl
index 4402e53e..d4032615 100644
--- a/scripts/names.jl
+++ b/scripts/names.jl
@@ -1,5 +1,5 @@
-using ArgParse, Downloads, Random, HORIZONS, NEOs
-using NEOs: numberofdays
+using ArgParse, Downloads, Random, NEOs, Dates
+using NEOs: RadecMPC, numberofdays, issatellite, reduce_tracklets
 
 # Potentially Hazardous Asteroids MPC list
 const PHAS_URL = "https://cgi.minorplanetcenter.net/cgi-bin/textversion.cgi?f=lists/PHAs.html"
@@ -20,37 +20,61 @@ function parse_commandline()
 
     @add_arg_table! s begin
         "--N", "-N"
-            help = "Number of NEOs"
+            help = "number of NEOs"
             arg_type = Int
-            default = 100
+            default = 0
         "--output", "-o"
-            help = "Output file"
+            help = "output file"
             arg_type = String
             default = "names.txt"
     end
 
     s.epilog = """
-        examples:\n
-        \n
-        julia --project names.jl -N 250 -o names.txt\n
+        note: if izero(N), then the script will save all compatible NEOs\n
         \n
+        examples:\n
+        # Save 1000 NEOs to names.txt\n
+        julia --project names.jl -N 1000 -o names.txt\n
     """
 
     return parse_args(s)
 end
 
+function isodcompatible(radec::Vector{RadecMPC{T}}) where {T <: Real}
+    # Eliminate observations before oficial discovery
+    firstobs = findfirst(r -> !isempty(r.discovery), radec)
+    isnothing(firstobs) && return false
+    radec = radec[firstobs:end]
+    # Filter out incompatible observations
+    filter!(radec) do r
+        hascoord(r.observatory) && !issatellite(r.observatory) &&
+        date(r) >= Date(2000)
+    end
+    length(radec) < 3 && return false
+    # There is at least one set of 3 tracklets with a < 15 days timespan
+    tracklets = reduce_tracklets(radec)
+    for i in 1:length(tracklets)-2
+        numberofdays(tracklets[i:i+2]) > 15.0 && continue
+        tracklets = tracklets[i:i+2]
+        radec = reduce(vcat, getfield.(tracklets, :radec))
+        sort!(radec)
+        break
+    end
+    return numberofdays(radec) <= 15.0
+end
+
 function main()
 
     # Parse arguments from commandline
     parsed_args = parse_commandline()
-
-    println("NEOs selector for orbitdetermination.jl")
-
     # Number of NEOs
     N::Int = parsed_args["N"]
-    println("• Number of NEOs: ", N)
     # Output file
-    outfile::String = parsed_args["output"]
+    output::String = parsed_args["output"]
+    # Print header
+    println("NEOs selector for orbitdetermination.jl")
+    println("• Number of NEOs: ", N)
+    println("• Output file: ", output)
 
     # Parse PHAs, Atens, Apollos and Amors
     provdesig = Vector{SubString{String}}(undef, 0)
@@ -71,34 +95,40 @@ function main()
     unique!(provdesig)
     # We can only process asteroids discovered between 2000 and 2024
     filter!(desig -> "2000" <= desig[1:4] <= "2024", provdesig)
-    # Shuffle the provisional designations list
-    shuffle!(provdesig)
+
+    if 0 < N < length(provdesig)
+        # Shuffle the provisional designations list
+        shuffle!(provdesig)
+    else
+        # Consider all compatible NEOs
+        N = length(provdesig)
+    end
 
     # Assemble a list with N NEOs
     names = Vector{String}(undef, N)
 
-    for i in eachindex(names)
-        while !isempty(provdesig)
-            neo = String(pop!(provdesig))
-            radec = fetch_radec_mpc("designation" => neo)
-            jplorbit = sbdb("des" => neo)["orbit"]
-            if (numberofdays(radec) <= 30.0) &&  (jplorbit["n_obs_used"] == length(radec)) &&
-                isnothing(jplorbit["n_del_obs_used"]) && isnothing(jplorbit["n_dop_obs_used"])
-                names[i] = neo
-                break
-            end
+    i = 1
+    while i <= N && !isempty(provdesig)
+        neo = String(popfirst!(provdesig))
+        radec = fetch_radec_mpc(neo)
+        if isodcompatible(radec)
+            names[i] = neo
+            i += 1
         end
     end
+    # Eliminate #undef elements
+    names = names[1:i-1]
+    # Sort provisional designations
     sort!(names)
 
     # Save names list
-    open(outfile, "w") do file
+    open(output, "w") do file
         for i in eachindex(names)
             write(file, names[i], i == N ? "" : "\n")
         end
     end
 
-    println("• Saved ", N, " names to ", outfile)
+    println("• Saved ", N, " names to ", output)
 
     nothing
 end
diff --git a/scripts/orbitdetermination.jl b/scripts/orbitdetermination.jl
index 5434bf29..fa5705eb 100644
--- a/scripts/orbitdetermination.jl
+++ b/scripts/orbitdetermination.jl
@@ -11,12 +11,12 @@ function parse_commandline()
     using Julia's interface for distributed computing."""
 
     @add_arg_table! s begin
-        "--names", "-n"
-            help = "File containing the names of the NEOs to be processed"
+        "--input", "-i"
+            help = "input names file"
             arg_type = String
             default = "names.txt"
         "--output", "-o"
-            help = "Output directory"
+            help = "output directory"
             arg_type = String
             default = pwd()
     end
@@ -25,7 +25,7 @@ function parse_commandline()
         example:\n
         \n
         # Run orbitdetermination.jl with 10 workers and 5 threads each\n
-        julia -p 10 -t 5 --project scripts/orbitdetermination.jl -n names.txt -o orbits/\n
+        julia -p 10 -t 5 --project scripts/orbitdetermination.jl -i names.txt -o orbits/\n
         \n
     """
 
@@ -34,174 +34,125 @@ end
 
 @everywhere begin
     using NEOs, Dates, JLD2
-    using NEOs: Tracklet, reduce_tracklets, isgauss
-
-    # Initial orbit determination routines
-
-    # First filter
-    function iod1(neo::String, outdir::String)
-        # Output file
-        filename = joinpath(outdir, replace(neo, " " => "_") * ".jld2")
-        # Download optical astrometry
-        radec = fetch_radec_mpc("designation" => neo)
-        if length(radec) < 3
-            jldsave(filename; tracklets = Vector{Tracklet}(undef, 0))
-            return false
+    using NEOs: AdmissibleRegion, RadecMPC, reduce_tracklets, numberofdays,
+        issatellite
+
+    function radecfilter(radec::Vector{RadecMPC{T}}) where {T <: Real}
+        # Eliminate observations before oficial discovery
+        firstobs = findfirst(r -> !isempty(r.discovery), radec)
+        isnothing(firstobs) && return false, radec
+        radec = radec[firstobs:end]
+        # Filter out incompatible observations
+        filter!(radec) do r
+            hascoord(r.observatory) && !issatellite(r.observatory) &&
+            date(r) >= Date(2000)
         end
-        # Parameters
-        params = NEOParameters(coeffstol = 10.0, bwdoffset = 0.007,
-                 fwdoffset = 0.007, jtlsiter = 10, adamhelp = false)
-        # Select at most three tracklets
+        length(radec) < 3 && return false, radec
+        # Find the first set of 3 tracklets with a < 15 days timespan
         tracklets = reduce_tracklets(radec)
-        if length(tracklets) > 3
-            tracklets = tracklets[1:3]
+        for i in 1:length(tracklets)-2
+            numberofdays(tracklets[i:i+2]) > 15.0 && continue
+            tracklets = tracklets[i:i+2]
             radec = reduce(vcat, getfield.(tracklets, :radec))
             sort!(radec)
+            break
         end
+        return numberofdays(radec) <= 15.0, radec
+    end
 
-        try
-            # Start of computation
-            init_time = now()
-            # Orbit determination
-            sol = orbitdetermination(radec, params)
-            # Time of computation
-            Δ = (now() - init_time).value
-            # Save orbit
-            if length(sol.res) != length(radec)
-                jldsave(filename; tracklets = tracklets, Δ = Δ)
-                return false
-            else
-                jldsave(filename; sol = sol, Δ = Δ)
-                return true
-            end
-        catch
-            # An error ocurred
-            jldsave(filename; tracklets = tracklets)
-            return false
-        end
+    # Default naive initial conditions for iod
+    function initcond(A::AdmissibleRegion{T}) where {T <: Real}
+        v_ρ = sum(A.v_ρ_domain) / 2
+        return [
+            (A.ρ_domain[1], v_ρ, :log),
+            (10^(sum(log10, A.ρ_domain) / 2), v_ρ, :log),
+            (sum(A.ρ_domain) / 2, v_ρ, :log),
+            (A.ρ_domain[2], v_ρ, :log),
+        ]
     end
 
-    # Second filter
-    function iod2(neo::String, outdir::String)
-        # Output from last filter
-        filename = joinpath(outdir, replace(neo, " " => "_") * ".jld2")
-        dict = JLD2.load(filename)
-        # Previous filter already computed an orbit
-        "sol" in keys(dict) && return true
-        # Check tracklets are non empty
-        tracklets = dict["tracklets"]
-        isempty(tracklets) && return false
-        # Computation time so far
-        if "Δ" in keys(dict)
-            Δ = dict["Δ"]
-        else
-            Δ = 0
-        end
-        # Optical astrometry
-        radec = reduce(vcat, getfield.(tracklets, :radec))
-        sort!(radec)
+    # Initial orbit determination routine
+    function iod(neo::String, filename::String)
+        # Download optical astrometry
+        radec = fetch_radec_mpc(neo)
+        # Get at most 3 tracklets for orbit determination
+        flag, radec = radecfilter(radec)
+        !flag && return false
+        # Dynamical function
+        dynamics = newtonian!
         # Parameters
-        params = NEOParameters(coeffstol = Inf, bwdoffset = 0.5,
-                 fwdoffset = 0.5, jtlsiter = 10, adamhelp = true)
-
-        try
-            # Start of computation
-            init_time = now()
-            # Orbit determination
-            sol = orbitdetermination(radec, params)
+        params = NEOParameters(
+            coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007,
+            gaussorder = 6, gaussQmax = 1.0,
+            adamiter = 500, adamQtol = 1e-5, tsaQmax = 2.0,
+            jtlsiter = 20, newtoniter = 10
+        )
+        # Start of computation
+        init_time = now()
+        # Initial orbit determination
+        sol = orbitdetermination(radec, params; dynamics, initcond)
+        # Termination condition
+        if length(sol.res) == length(radec) && nrms(sol) < 1.5
             # Time of computation
-            Δ += (now() - init_time).value
+            Δ = (now() - init_time).value
             # Save orbit
-            if length(sol.res) != length(radec)
-                jldsave(filename; tracklets = tracklets, Δ = Δ)
-                return false
-            else
-                jldsave(filename; sol = sol, Δ = Δ)
-                return true
-            end
-        catch
-            # An error ocurred
-            jldsave(filename; tracklets = tracklets)
-            return false
-        end
-    end
-
-    # Third filter
-    function iod3(neo::String, outdir::String)
-        # Output from last filter
-        filename = joinpath(outdir, replace(neo, " " => "_") * ".jld2")
-        dict = JLD2.load(filename)
-        # Previous filter already computed an orbit
-        "sol" in keys(dict) && return true
-        # Check tracklets are non empty
-        tracklets = dict["tracklets"]
-        isempty(tracklets) && return false
-        # Computation time so far
-        if "Δ" in keys(dict)
-            Δ = dict["Δ"]
-        else
-            Δ = 0
+            jldsave(filename; sol = sol, Δ = Δ)
+            # Sucess flag
+            return true
         end
-        # Optical astrometry
-        radec = reduce(vcat, getfield.(tracklets, :radec))
-        sort!(radec)
         # Parameters
-        params = NEOParameters(coeffstol = Inf, bwdoffset = 0.5,
-                 fwdoffset = 0.5, jtlsiter = 10, adamhelp = true)
-
-        try
-            # Start of computation
-            init_time = now()
-            # Orbit determination
-            if isgauss(tracklets)
-                sol = tooshortarc(radec, tracklets, params)
-            else
-                sol = gaussinitcond(radec, tracklets, params)
-            end
-            # Time of computation
-            Δ += (now() - init_time).value
-            # Save orbit
-            if length(sol.res) != length(radec)
-                jldsave(filename; tracklets = tracklets, Δ = Δ)
-                return false
-            else
-                jldsave(filename; sol = sol, Δ = Δ)
-                return true
-            end
-        catch
-            # An error ocurred
-            jldsave(filename; tracklets = tracklets)
-            return false
+        params = NEOParameters(params;
+            coeffstol = 10.0, bwdoffset = 0.007, fwdoffset = 0.007,
+            gaussorder = 6, gaussQmax = 1.0,
+            adamiter = 200, adamQtol = 0.01, tsaQmax = 2.0,
+            jtlsiter = 20, newtoniter = 5
+        )
+        # Initial orbit determination
+        _sol_ = orbitdetermination(radec, params; dynamics, initcond)
+        # Time of computation
+        Δ = (now() - init_time).value
+        # Choose best orbit
+        if length(sol.res) == length(_sol_.res) == length(radec)
+            sol = min(sol, _sol_)
+        elseif length(sol.res) == length(radec)
+            sol = sol
+        elseif length(_sol_.res) == length(radec)
+            sol = _sol_
+        else
+            sol = zero(NEOSolution{Float64, Float64})
         end
+        # Save orbit
+        iszero(sol) && return false
+        jldsave(filename; sol = sol, Δ = Δ)
+        return true
     end
 end
 
 function main()
     # Parse arguments from commandline
     parsed_args = parse_commandline()
-
-    println("Orbit determination for NEOs via jet transport")
-    println()
-
-    # NEOs to be processed
-    namesfile = parsed_args["names"]
-    neos = readlines(namesfile)
-    println(length(neos), " NEOs to be processed with ", nworkers(), " workers (",
-            Threads.nthreads(), " threads each)")
-    println()
+    # Input names file
+    input::String = parsed_args["input"]
     # Output directory
-    outdir = parsed_args["output"]
+    output::String = parsed_args["output"]
+    # Print header
+    println("Initial orbit determination for NEOs via jet transport")
+    println("• Input names file: ", input)
+    println("• Output directory: ", output)
+
+    # Parse NEOs' designations
+    neos = readlines(input)
+    println("• ", length(neos), " NEOs to be processed with ", nworkers(),
+        " workers (", Threads.nthreads(), " threads each)")
+    # Output files
+    filenames = map(neos) do neo
+        return joinpath(output, replace(neo, " " => "") * ".jld2")
+    end
 
     # Distributed orbit determination
+    mask = pmap(iod, neos, filenames; on_error = ex -> false)
+    println("• ", count(mask), " / ", length(neos), " successful NEOs")
 
-    # First filter
-    mask = pmap(neo -> iod1(neo, outdir), neos)
-    all(mask) && return nothing
-    # Second filter
-    mask = pmap(neo -> iod2(neo, outdir), neos)
-    all(mask) && return nothing
-    # Third filter
-    mask = pmap(neo -> iod3(neo, outdir), neos)
     return nothing
 end
 
diff --git a/src/NEOs.jl b/src/NEOs.jl
index 12456c6b..2026b2f2 100644
--- a/src/NEOs.jl
+++ b/src/NEOs.jl
@@ -2,7 +2,7 @@ module NEOs
 
 # __precompile__(false)
 
-import Base: show, string, isless, convert, zero, iszero, isnan, in
+import Base: show, string, isless, convert, zero, iszero, isnan, in, min
 import Tables: istable, rowaccess, rows, schema, Schema
 import SatelliteToolboxTransformations: sv_ecef_to_eci, sv_ecef_to_ecef, ecef_to_geocentric
 import JLD2: writeas
@@ -12,6 +12,7 @@ using Dates, InteractiveUtils, LazyArtifacts, LinearAlgebra, Printf, JSON,
       TaylorSeries, SatelliteToolboxTransformations, TaylorIntegration,
       SPICE, JLD2, Scratch
 using AutoHashEquals.Compat
+using Dates: epochms2datetime
 using Downloads: download
 using DelimitedFiles: readdlm
 using HTTP: get
@@ -47,8 +48,9 @@ export hasdelay, hasdoppler, delay, doppler, rcvr, xmit, read_radar_jpl, write_r
 # NEOCPObject
 export fetch_objects_neocp, get_radec_neocp, fetch_radec_neocp, get_orbits_neocp
 # Units
-export julian2etsecs, etsecs2julian, dtutc2et, dtutc2jdtdb, et_to_200X, days_to_200X,
-       datetime_to_200X, datetime2days, days2datetime, rad2arcsec, arcsec2rad, mas2rad
+export julian2etsecs, etsecs2julian, dtutc2et, dtutc2jdtdb, et2dtutc, jdtdb2dtutc,
+       et_to_200X, days_to_200X, dtutc_to_200X, dtutc2days, days2dtutc,
+       rad2arcsec, arcsec2rad, mas2rad
 # JPL ephemerides
 export loadjpleph, sunposvel, earthposvel, moonposvel, apophisposvel197, apophisposvel199
 # PE and NEOs ephemerides
@@ -71,11 +73,11 @@ export project, chi2, nms, nrms, diffcorr, newtonls, levenbergmarquardt, tryls
 # NEOSolution
 export epoch, sigmas, snr, jplcompare, uncertaintyparameter
 # Too Short Arc
-export tooshortarc
+export tsaiod
 # Gauss method
-export gauss_method, gaussinitcond
+export gauss_method, gaussiod
 # Orbit determination
-export curvature, issinglearc, isgauss, orbitdetermination
+export curvature, issinglearc, orbitdetermination
 # B plane
 export valsecchi_circle, bopik, mtp
 # Outlier rejection
diff --git a/src/constants.jl b/src/constants.jl
index a92d303c..e30a931c 100644
--- a/src/constants.jl
+++ b/src/constants.jl
@@ -13,6 +13,9 @@ const poteph::TaylorInterpolant{Float64, Float64, 2, Vector{Float64}, Matrix{Tay
 const ttmtdb::TaylorInterpolant{Float64, Float64, 1, Vector{Float64}, Vector{Taylor1{Float64}}} = TaylorInterpolant(sseph.t0, sseph.t, sseph.x[:,end])
 const SSEPHORDER::Int = get_order(sseph.x[1])
 
+# Milliseconds between rounding epoch and J2000
+const EPOCHMSJ2000::Int = (DateTime(2000, 1, 1, 12) - DateTime(0)).value
+
 # Earth orientation parameters (eop) 2000
 const eop_IAU2000A::EopIau2000A = fetch_iers_eop(Val(:IAU2000A))
 
@@ -138,13 +141,12 @@ const CATALOGUES_MPC_URL = "https://www.minorplanetcenter.net/iau/info/Catalogue
 # MPC observatories file url
 const OBSERVATORIES_MPC_URL = "https://www.minorplanetcenter.net/iau/lists/ObsCodes.html"
 
-# MPC database search url
-const search_mpc_url = "https://www.minorplanetcenter.net/db_search/show_object?utf8=%E2%9C%93&object_id="
-# MPC observations url
-const obs_mpc_url = "https://www.minorplanetcenter.net/tmp2/"
-
+# MPC Oservations API url
+const MPC_OBS_API_URL = "https://data.minorplanetcenter.net/api/get-obs"
 # NEO Confirmation Page File URL
 const NEOCP_FILE_URL = "https://www.minorplanetcenter.net/Extended_Files/neocp.json"
+# MPC NEOCP Oservations API url
+const MPC_NEOCP_OBS_API_URL = "https://data.minorplanetcenter.net/api/get-obs-neocp"
 # NEO Confirmation Page Show Orbits URL
 const NEOCP_SHOWORBS_URL = "https://cgi.minorplanetcenter.net/cgi-bin/showobsorbs.cgi"
 
diff --git a/src/observations/neocp.jl b/src/observations/neocp.jl
index 505012c1..084a607a 100644
--- a/src/observations/neocp.jl
+++ b/src/observations/neocp.jl
@@ -117,59 +117,65 @@ function fetch_objects_neocp()
 end
 
 @doc raw"""
-    get_radec_neocp(id::String, filename::String = id * ".txt")
+    fetch_radec_neocp(id::AbstractString)
 
-Download MPC optical astrometry of NEOCP [1] object `id` and save the output to `filename`.
+Download MPC optical astrometry of NEOCP object `id` and
+return the output as `Vector{RadecMPC{Float64}}`.
 
 !!! reference
-    [1] https://www.minorplanetcenter.net/iau/NEO/toconfirm_tabular.html.
+    MPC NEOCP Observations API documentation:
+
+    https://www.minorplanetcenter.net/mpcops/documentation/neocp-observations-api/.
 """
-function get_radec_neocp(id::String, filename::String = id * ".txt")
-    # Assemble URL
-    url = string(NEOCP_SHOWORBS_URL, "?Obj=", id, "&obs=y")
-    # HTTP response
-    resp = get(url)
-    # Parse response body as String
+function fetch_radec_neocp(id::AbstractString)
+    # HTTP parameters
+    params = JSON.json(Dict("trksubs" => [id], "output_format" => ["OBS80"]))
+    resp = get(MPC_NEOCP_OBS_API_URL, ("Content-Type" => "application/json",), params)
+    # Convert to String
     text = String(resp.body)
-    # Parse optical astrometry
-    radec = read_radec_mpc(text)
-    # Save result to filename
-    write_radec_mpc(radec, filename)
+    # Parse JSON
+    dict = JSON.parse(text)
+    # Parse observations
+    radec = RadecMPC.(eachmatch(RADEC_MPC_REGEX, dict[1]["OBS80"]))
+    # Eliminate unsuccessful matches and repeated entries
+    filter!(r -> isa(r, RadecMPC{Float64}), radec)
+    unique!(radec)
 
-    return filename
+    return radec
 end
 
 @doc raw"""
-    fetch_radec_neocp(id::String)
+    get_radec_neocp(id::AbstractString [, filename::AbstractString])
 
-Download MPC optical astrometry of NEOCP [1] object `id` and
-return the output as `Vector{RadecMPC{Float64}}`.
+Download MPC optical astrometry of NEOCP object `id` and save the output to `filename`
+(default: `id.txt`).
 
 !!! reference
-    [1] https://www.minorplanetcenter.net/iau/NEO/toconfirm_tabular.html.
+    MPC NEOCP Observations API documentation:
+
+    https://www.minorplanetcenter.net/mpcops/documentation/neocp-observations-api/.
 """
-function fetch_radec_neocp(id::String)
-    # Temporary file
-    f = tempname()
-    # Download optical astrometry
-    get_radec_neocp(id, f)
-    # Parse optical astrometry
-    radec = read_radec_mpc(f)
-    # Delete temporary file
-    rm(f)
+function get_radec_neocp(id::AbstractString, filename::AbstractString =
+    replace(id, " " => "_") * ".txt")
+    # Fetch optical astrometry
+    radec = fetch_radec_neocp(id)
+    # Write observations to file
+    write_radec_mpc(radec, filename)
 
-    return radec
+    return filename
 end
 
 @doc raw"""
-    get_orbits_neocp(id::String, filename::String = id * ".txt")
+    get_orbits_neocp(id::AbstractString [, filename::AbstractString])
 
-Download sample orbits of NEOCP [1] object `id` and save the output to `fiilename`.
+Download sample orbits of NEOCP [1] object `id` and save the output to `filename`
+(default: `id.txt`).
 
 !!! reference
     [1] https://www.minorplanetcenter.net/iau/NEO/toconfirm_tabular.html.
 """
-function get_orbits_neocp(id::String, filename::String = id * ".txt")
+function get_orbits_neocp(id::AbstractString, filename::AbstractString =
+    replace(id, " " => "_") * ".txt")
     # Assemble URL
     url = string(NEOCP_SHOWORBS_URL, "?Obj=", id, "&orb=y")
     # HTTP response
diff --git a/src/observations/process_radec.jl b/src/observations/process_radec.jl
index 4c9e1271..4ecd63bf 100644
--- a/src/observations/process_radec.jl
+++ b/src/observations/process_radec.jl
@@ -449,7 +449,7 @@ function w8sveres17(obs::RadecMPC{T}) where {T <: AbstractFloat}
     # Table 2: epoch-dependent astrometric residuals
     if obscode == "703"
         return Date(dt_utc_obs) < Date(2014,1,1) ? w : 0.8w
-    elseif obscode == "691"
+    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
@@ -458,7 +458,7 @@ function w8sveres17(obs::RadecMPC{T}) where {T <: AbstractFloat}
         return w
     elseif obscode == "G96"
         return 0.5w
-    elseif obscode == "F51"
+    elseif obscode ∈ ("F51", "F52") # Pan-STARRS 1 & 2, Haleakala, Hawaii
         return 0.2w
     elseif obscode ∈ ("G45", "608")
         return 0.6w
@@ -469,7 +469,8 @@ function w8sveres17(obs::RadecMPC{T}) where {T <: AbstractFloat}
     # 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
+    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
diff --git a/src/observations/radec_mpc.jl b/src/observations/radec_mpc.jl
index de5dbfb3..65ee8aeb 100644
--- a/src/observations/radec_mpc.jl
+++ b/src/observations/radec_mpc.jl
@@ -384,49 +384,52 @@ function write_radec_mpc(obs::Vector{RadecMPC{T}}, filename::String) where {T <:
 end
 
 @doc raw"""
-    get_radec_mpc(id::Pair{String, String}, filename::String = replace(id[2], " " => "_") * ".txt")
+    fetch_radec_mpc(id::AbstractString)
 
-Download MPC optical astrometry of NEO `id` and save the output to `filename`.
+Download MPC optical astrometry of NEO `id` and return the output
+as `Vector{RadecMPC{Float64}}`.
+
+!!! reference
+    MPC Observations API documentation:
+
+    https://minorplanetcenter.net/mpcops/documentation/observations-api/.
 """
-function get_radec_mpc(id::Pair{String, String}, filename::String = replace(id[2], " " => "_") * ".txt")
-    # HTTP query
-    query = ["table" => "observations", id]
-    resp = get("http://minorplanetcenter.net/search_db"; query = query)
-    # Converty to String
+function fetch_radec_mpc(id::AbstractString)
+    # HTTP parameters
+    params = JSON.json(Dict("desigs" => [id], "output_format" => ["OBS80"]))
+    resp = get(MPC_OBS_API_URL, ("Content-Type" => "application/json",), params)
+    # Convert to String
     text = String(resp.body)
     # Parse JSON
-    obs = JSON.parse(text)
-    # Find matches
-    matches = Vector{Union{RegexMatch, Nothing}}(undef, length(obs))
-    for i in eachindex(obs)
-        s = obs[i]["original_record"]
-        matches[i] = match(RADEC_MPC_REGEX, s)
-    end
-    filter!(!isnothing, matches)
-    # Parse RadecMPC
-    radec = RadecMPC.(matches)
-    # Write observations to file
-    write_radec_mpc(radec, filename)
+    dict = JSON.parse(text)
+    # Parse observations
+    radec = RadecMPC.(eachmatch(RADEC_MPC_REGEX, dict[1]["OBS80"]))
+    # Eliminate unsuccessful matches and repeated entries
+    filter!(r -> isa(r, RadecMPC{Float64}), radec)
+    unique!(radec)
 
-    return filename
+    return radec
 end
 
 @doc raw"""
-    fetch_radec_mpc(id::Pair{String, String})
+    get_radec_mpc(id::AbstractString [, filename::AbstractString])
+
+Download MPC optical astrometry of NEO `id` and save the output to `filename`
+(default: `id.txt`).
+
+!!! reference
+    MPC Observations API documentation:
 
-Download MPC optical astrometry of NEO `id` and return the output as `Vector{RadecMPC{Float64}}`.
+    https://minorplanetcenter.net/mpcops/documentation/observations-api/.
 """
-function fetch_radec_mpc(id::Pair{String, String})
-    # Temporary file
-    f = tempname()
-    # Download optical astrometry
-    get_radec_mpc(id, f)
-    # Parse optical astrometry
-    radec = read_radec_mpc(f)
-    # Delete temporary file
-    rm(f)
+function get_radec_mpc(id::AbstractString, filename::AbstractString =
+    replace(id, " " => "_") * ".txt")
+    # Fetch optical astrometry
+    radec = fetch_radec_mpc(id)
+    # Write observations to file
+    write_radec_mpc(radec, filename)
 
-    return radec
+    return filename
 end
 
 # Methods to convert a Vector{<:AbstractAstrometry} to a DataFrame
diff --git a/src/observations/tracklet.jl b/src/observations/tracklet.jl
index 0c00cf62..b47635d1 100644
--- a/src/observations/tracklet.jl
+++ b/src/observations/tracklet.jl
@@ -120,13 +120,13 @@ function Tracklet(radec::Vector{RadecMPC{T}}, df::AbstractDataFrame) where {T <:
     df = combine(gdf, [:α, :δ, :mag] .=> mean, renamecols = false)
 
     # Julian days of observation
-    df.t_julian = datetime2julian.(df.date)
+    df.t_julian = dtutc2jdtdb.(df.date)
     # Days of observation [relative to first observation]
     df.t_rel = df.t_julian .- df.t_julian[1]
     # Mean date [relative to first observation]
     t_mean = mean(df.t_rel)
     # Mean date [DateTime]
-    date = julian2datetime(df.t_julian[1] + t_mean)
+    date = jdtdb2dtutc(df.t_julian[1] + t_mean)
 
     # Points in top quarter
     N_top = count(x -> x > 3π/2, df.α)
@@ -255,7 +255,7 @@ end
 
 function curvature(radec::AbstractVector{RadecMPC{T}}) where {T <: Real}
     # Days of observation [julian days]
-    ts = datetime2julian.(date.(radec))
+    ts = dtutc2jdtdb.(date.(radec))
     # Right ascension
     αs = ra.(radec)
     # Declination
diff --git a/src/observations/units.jl b/src/observations/units.jl
index 84d75e81..344ba0ea 100644
--- a/src/observations/units.jl
+++ b/src/observations/units.jl
@@ -12,7 +12,7 @@ end
 @doc raw"""
     etsecs2julian(et::T) where {T <: Number}
 
-Convert `et` TDB seconds since J2000 to Julian date (TDB).
+Convert `et` TDB seconds since J2000 to Julian date in TDB time scale.
 
 See also [`julian2etsecs`](@ref).
 """
@@ -53,40 +53,64 @@ function dtutc2jdtdb(dtutc::DateTime)
     return etsecs2julian(et) # JDTDB
 end
 
+@doc raw"""
+    et2dtutc(et::T) where {T <: Number}
+
+Convert `et` TDB seconds past the J2000 epoch to a UTC `DateTime`.
+"""
+function et2dtutc(et::T) where {T <: Number}
+    # UTC seconds past J2000
+    utc = et - tdb_utc(et)
+    # UTC milliseconds past J2000
+    mill = Millisecond(round(Int, utc * 1_000)).value
+    # UTC DateTime
+    return epochms2datetime(EPOCHMSJ2000 + mill)
+end
+
+@doc raw"""
+    jdtdb2dtutc(jdtdb::T) where {T <: Number}
+
+Convert `jdtdb`, a Julian date in TDB time scale, to a UTC `DateTime`.
+"""
+function jdtdb2dtutc(jdtdb::T) where {T <: Number}
+    et = julian2etsecs(jdtdb) # TDB seconds since J2000.0
+    return et2dtutc(et) # UTC DateTime
+end
+
 @doc raw"""
     et_to_200X(et::T) where {T <: Number}
 
-Convert `et` ephemeris seconds since J2000 to years `200X`.
+Convert `et` TDB seconds since J2000 to TDB years `200X`.
 """
 et_to_200X(et::T) where {T <: Number} = 2000 + et/daysec/yr
 
 @doc raw"""
     days_to_200X(d::T) where {T <: Number}
 
-Convert `d` days since J2000 to years `200X`.
+Convert `d` TDB days since J2000 to TDB years `200X`.
 """
 days_to_200X(d::T) where {T <: Number} = 2000 + d/yr
 
 @doc raw"""
-    datetime_to_200X(x::DateTime)
+    dtutc_to_200X(x::DateTime)
 
-Convert `DateTime` `x` to years `200X`.
+Convert `x`, a UTC `DateTime`, to TDB years `200X`.
 """
-datetime_to_200X(x::DateTime) = et_to_200X(dtutc2et(x))
+dtutc_to_200X(x::DateTime) = et_to_200X(dtutc2et(x))
 
 @doc raw"""
-    datetime2days(x::DateTime)
+    dtutc2days(x::DateTime)
 
-Convert `DateTime` `x` to days since J2000.
+Convert `x`, a UTC `DateTime`, to TDB days since J2000.
 """
-datetime2days(x::DateTime) = datetime2julian(x) - JD_J2000
+dtutc2days(x::DateTime) = dtutc2jdtdb(x) - JD_J2000
 
 @doc raw"""
-    days2datetime(d::T) where {T <: Number}
+    days2dtutc(d::T) where {T <: Number}
 
-Convert `d` days since J2000 to `DateTime`.
+Convert `d` TDB days since J2000 to a UTC `DateTime`.
 """
-days2datetime(d::T) where {T <: Number} = julian2datetime(d + JD_J2000)
+days2dtutc(d::T) where {T <: Number} = jdtdb2dtutc(d + JD_J2000)
 
 @doc raw"""
     tdb_utc(et::T) where {T <: Number}
@@ -170,7 +194,7 @@ Convert radians to arcseconds.
 
 See also [`arcsec2rad`](@ref) and [`mas2rad`](@ref).
 """
-rad2arcsec(x::T) where {T <: Number} = 3600 * rad2deg(x) # rad2deg(rad) -> deg; 3600 * deg -> arcsec
+rad2arcsec(x::T) where {T <: Number} = 3600 * rad2deg(x)
 
 @doc raw"""
     arcsec2rad(x::T) where {T <: Number}
@@ -179,7 +203,7 @@ Convert arcseconds to radians.
 
 See also [`rad2arcsec`](@ref) and [`mas2rad`](@ref).
 """
-arcsec2rad(x::T) where {T <: Number} = deg2rad(x / 3600) # arcsec/3600 -> deg; deg2rad(deg) -> rad
+arcsec2rad(x::T) where {T <: Number} = deg2rad(x / 3600)
 
 @doc raw"""
     mas2rad(x::T) where {T <: Number}
@@ -188,4 +212,4 @@ Convert milli-arcseconds to radians.
 
 See also [`rad2arcsec`](@ref) and [`arcsec2rad`](@ref).
 """
-mas2rad(x::T) where {T <: Number} = arcsec2rad(x / 1000) # mas/1000 -> arcsec; arcsec2rad(arcsec) -> rad
\ No newline at end of file
+mas2rad(x::T) where {T <: Number} = arcsec2rad(x / 1000)
\ No newline at end of file
diff --git a/src/orbitdetermination/admissibleregion.jl b/src/orbitdetermination/admissibleregion.jl
index 94608580..e647ec76 100644
--- a/src/orbitdetermination/admissibleregion.jl
+++ b/src/orbitdetermination/admissibleregion.jl
@@ -82,7 +82,7 @@ function AdmissibleRegion(tracklet::Tracklet{T}, params::NEOParameters{T}) where
     # Topocentric unit vector and partials
     ρ, ρ_α, ρ_δ = topounitpdv(α, δ)
     # Time of observation [days since J2000]
-    t_days = datetime2days(t_datetime)
+    t_days = dtutc2days(t_datetime)
     # Time of observation [et seconds]
     t_et = dtutc2et(t_datetime)
     # Heliocentric position of the observer
@@ -595,10 +595,10 @@ Convert topocentric range/range-rate `[ρ, v_ρ]` to barycentric cartesian coord
 """
 function topo2bary(A::AdmissibleRegion{T}, ρ::U, v_ρ::U) where {T <: Real, U <: Number}
     # Barycentric position
-    r = A.q[1:3] + ρ * A.ρ_unit + sseph(su, datetime2days(A.date))[1:3]
+    r = A.q[1:3] + ρ * A.ρ_unit + sseph(su, dtutc2days(A.date))[1:3]
     # Barycentric velocity
     v = A.q[4:6] + v_ρ * A.ρ_unit + ρ * A.v_α * A.ρ_α + ρ * A.v_δ * A.ρ_δ
-        + sseph(su, datetime2days(A.date))[4:6]
+        + sseph(su, dtutc2days(A.date))[4:6]
     # Barycentric state vector
     return vcat(r, v)
 end
@@ -611,7 +611,7 @@ Convert barycentric cartesian coordinates `q0` to topocentric range/range-rate.
 """
 function bary2topo(A::AdmissibleRegion{T}, q0::Vector{U}) where {T <: Real, U <: Number}
     # Heliocentric state vector
-    r = q0 - sseph(su, datetime2days(A.date))
+    r = q0 - sseph(su, dtutc2days(A.date))
     # Topocentric range
     ρ = euclid3D(r - A.q)
     # Topocentric range rate
@@ -632,7 +632,7 @@ function attr2bary(A::AdmissibleRegion{T}, a::Vector{U},
     # Unfold
     α, δ, v_α, v_δ, ρ, v_ρ = a
     # Light-time correction to epoch
-    t = datetime2days(A.date) - ρ/c_au_per_day
+    t = dtutc2days(A.date) - ρ/c_au_per_day
     # TO DO: `et::TaylorN` is too slow for `adam` due to
     # SatelliteToolboxTransformations overloads in src/observations/topocentric.jl
     et = dtutc2et(A.date) - cte(cte(ρ))/c_au_per_sec
diff --git a/src/orbitdetermination/gaussinitcond.jl b/src/orbitdetermination/gaussinitcond.jl
index 0d8d9c08..bed54f42 100644
--- a/src/orbitdetermination/gaussinitcond.jl
+++ b/src/orbitdetermination/gaussinitcond.jl
@@ -339,67 +339,34 @@ function gauss_method(observatories::Vector{ObservatoryMPC{T}}, dates::Vector{Da
 end
 
 @doc raw"""
-    numberofdays(dates::Vector{DateTime})
-    numberofdays(dates::Vector{RadecMPC{T}}) where {T <: Real}
-    numberofdays(dates::Vector{Tracklet{T}}) where {T <: Real}
+    numberofdays(x)
 
-Return the timespan of `dates` in days. The function assumes `dates` is sorted.
+Return the timespan of `x::AbstractVector{T}` in days, where
+`T` can be `DateTime`, `RadecMPC` or `Tracklet`.
 """
-numberofdays(dates::Vector{DateTime}) = (dates[end] - dates[1]).value / 86_400_000
-
-function numberofdays(dates::Vector{RadecMPC{T}}) where {T <: Real}
-    return (dates[end].date - dates[1].date).value / 86_400_000
+function numberofdays(dates::AbstractVector{DateTime})
+    t0, tf = extrema(dates)
+    return (tf - t0).value / 86_400_000
 end
 
-function numberofdays(dates::Vector{Tracklet{T}}) where {T <: Real}
-    return (dates[end].radec[end].date - dates[1].radec[1].date).value / 86_400_000
+function numberofdays(radec::AbstractVector{RadecMPC{T}}) where {T <: Real}
+    t0, tf = extrema(date, radec)
+    return (tf - t0).value / 86_400_000
 end
 
-# Return the `maxtriplets` best combinations of three `tracklets`
-# for Gauss' method. See the line after equation (27) of
-# https://doi.org/10.1016/j.icarus.2007.11.033
-function _gausstriplets1(tracklets::Vector{Tracklet{T}}, maxtriplets::Int) where {T <: Real}
-    # Number of tracklets
-    L = length(tracklets)
-    # All possible triplets of indices
-    CI = CartesianIndices((1:L-2, 2:L-1, 3:L))
-    # Allocate memory
-    triplets = zeros(Int, 3, maxtriplets)
-    τ = fill(T(Inf), maxtriplets)
-    # Check all possible triplets of indices
-    for ci in CI
-        # Indices cannot repeat
-        !allunique(ci.I) && continue
-        # Unfold indices
-        i, j, k = ci.I
-        # Timespan must be at least one day
-        Δ = (tracklets[k].date - tracklets[i].date).value
-        Δ < 86_400_000 && continue
-        # Absolute difference between τ1 and τ3
-        τ1 = (tracklets[j].date - tracklets[i].date).value
-        τ3 = (tracklets[k].date - tracklets[j].date).value
-        _τ_ = abs(τ3 - τ1)
-        # Current max τ
-        τmax, n = findmax(τ)
-        # Update triplets and τ
-        if _τ_ < τmax
-            triplets[:, n] .= i, j, k
-            τ[n] = _τ_
-        end
-    end
-    # Remove Infs and sort
-    idxs = findall(!isinf, τ)
-    perm = sortperm(view(τ, idxs))
-
-    return triplets[:, perm], τ[perm]
+function numberofdays(tracklets::AbstractVector{Tracklet{T}}) where {T <: Real}
+    dates = map(t -> extrema(date, t.radec), tracklets)
+    t0, tf = minimum(first, dates), maximum(last, dates)
+    return (tf - t0).value / 86_400_000
 end
 
 # Given a vector of three `tracklets`, update `observatories`, `dates`,
 # `α` and `δ` with the best combination of three observations for Gauss'
 # method. See the line after equation (27) of
 # https://doi.org/10.1016/j.icarus.2007.11.033
-function _gausstriplets2!(observatories::Vector{ObservatoryMPC{T}}, dates::Vector{DateTime}, α::Vector{T},
-                          δ::Vector{T}, tracklets::AbstractVector{Tracklet{T}}) where {T <: Real}
+function _gausstriplet!(observatories::Vector{ObservatoryMPC{T}},
+    dates::Vector{DateTime}, α::Vector{T}, δ::Vector{T},
+    tracklets::AbstractVector{Tracklet{T}}) where {T <: Real}
     # Unfold tracklets
     a, b, c = tracklets
     # All possible triplets of indices
@@ -428,20 +395,26 @@ function _gausstriplets2!(observatories::Vector{ObservatoryMPC{T}}, dates::Vecto
     dates .= a.radec[i].date, b.radec[j].date, c.radec[k].date
     α .= a.radec[i].α, b.radec[j].α, c.radec[k].α
     δ .= a.radec[i].δ, b.radec[j].δ, c.radec[k].δ
+    # Sort triplet
+    idxs = sortperm(dates)
+    permute!(observatories, idxs)
+    permute!(dates, idxs)
+    permute!(α, idxs)
+    permute!(δ, idxs)
 
     return nothing
 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!(q::Vector{TaylorN{T}}, jd0::T, tracklet::Tracklet,
+    params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D}
     # Exploratory propagation
     bwd, fwd, _ = propres(tracklet.radec, jd0, q(), params)
     # Admissible region
     A = AdmissibleRegion(tracklet, params)
     # Epoch [days since J2000]
-    At0 = datetime2days(A.date)
+    At0 = dtutc2days(A.date)
     # Barycentric cartesian initial condition
     if At0 <= jd0 - JD_J2000
         q0 = bwd(At0)
@@ -455,101 +428,80 @@ function _adam!(q::Vector{TaylorN{T}}, jd0::T, tracklet::Tracklet, params::NEOPa
     # ADAM
     ae, _ = adam(tracklet.radec, A, ρ, v_ρ, params; scale = :log, dynamics = dynamics)
     # Epoch [julian days] (corrected for light-time)
-    jd0 = datetime2julian(A.date) - ae[5] / c_au_per_day
+    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] + (q0[i] / 10^5) * TaylorN(i, order = params.tsaorder) for i in 1:6]
+    q .= [q0[i] + (abs(q0[i]) / 10^5) * TaylorN(i, order = params.tsaorder) for i in 1:6]
 
     return jd0
 end
 
 @doc raw"""
-    gaussinitcond(radec::Vector{RadecMPC{T}}, [tracklets::Vector{Tracklet{T}},]
-                  params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D}
+    gaussiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}},
+        params::NEOParameters{T}; kwargs...) where {T <: Real, D}
 
-Compute an orbit via Jet Transport Gauss' Method.
+Fit a preliminary orbit to `radec` via jet transport Gauss method.
 
 See also [`gauss_method`](@ref).
 
-# Arguments
+## Arguments
 
 - `radec::Vector{RadecMPC{T}}`: vector of optical astrometry.
 - `tracklets::Vector{Tracklet{T}},`: vector of tracklets.
 - `params::NEOParameters{T}`: see `Gauss' Method Parameters` of [`NEOParameters`](@ref).
-- `dynamics::D`: dynamical model.
 
-!!! warning
-    This function will change the (global) `TaylorSeries` variables.
-"""
-function gaussinitcond(radec::Vector{RadecMPC{T}}, params::NEOParameters{T};
-                       dynamics::D = newtonian!) where {T <: Real, D}
-    # Reduce tracklets by polynomial regression
-    tracklets = reduce_tracklets(radec)
-    # Set jet transport variables
-    varorder = max(params.tsaorder, params.gaussorder)
-    scaled_variables("dx", ones(T, 6); order = varorder)
-    # Gauss' Method
-    return gaussinitcond(radec, tracklets, params; dynamics)
-end
+## Keyword arguments
 
-function gaussinitcond(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}},
-                       params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D}
+- `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}
+    # Allocate memory for orbit
+    sol = zero(NEOSolution{T, T})
+    # This function requires exactly 3 tracklets
+    length(tracklets) != 3 && return sol
+    # Jet transport scaling factors
+    scalings = fill(1e-6, 6)
+    # Jet transport perturbation (ra/dec)
+    dq = [scalings[i] * TaylorN(i, order = params.gaussorder) for i in 1:6]
     # gauss_method input
     observatories = Vector{ObservatoryMPC{T}}(undef, 3)
     dates = Vector{DateTime}(undef, 3)
     α = Vector{T}(undef, 3)
     δ = Vector{T}(undef, 3)
-    # Observations triplets
-    triplets, _ = _gausstriplets1(tracklets, params.max_triplets)
-    # Jet transport scaling factors
-    scalings = fill(1e-6, 6)
-    # Jet transport perturbation (ra/dec)
-    dq = [scalings[i] * TaylorN(i, order = params.gaussorder) for i in 1:6]
-    # Best orbit
-    best_sol = zero(NEOSolution{T, T})
-    best_Q = T(Inf)
-    # Convergence flag
-    flag = false
-    # Iterate over triplets
-    for triplet in eachcol(triplets)
-        # Find best triplet of observations
-        _gausstriplets2!(observatories, dates, α, δ, view(tracklets, triplet))
-        # Julian day of middle observation
-        _jd0_ = datetime2julian(dates[2])
-        # Gauss method solution
-        solG = gauss_method(observatories, dates, α .+ dq[1:3], δ .+ dq[4:6], params)
-        # Filter non-physical (negative) rho solutions
-        filter!(x -> all(cte.(x.ρ) .> 0), solG)
-        # Filter Gauss solutions by heliocentric energy
-        filter!(x -> cte(cte(heliocentric_energy(x.statevect))) <= 0, solG)
-        # Iterate over Gauss solutions
-        for i in eachindex(solG)
-            # Epoch (corrected for light-time)
-            jd0 = _jd0_ - cte(solG[i].ρ[2]) / c_au_per_day
-            # Jet transport initial conditions
-            q = solG[i].statevect .+ params.eph_su(jd0 - JD_J2000)
-            # ADAM (if requested by the user)
-            if params.adamhelp
-                jd0 = _adam!(q, jd0, tracklets[triplet[2]], params; dynamics)
-            end
-            # Jet transport least squares
-            sol = jtls(radec, tracklets, jd0, q, triplet[2], params; dynamics)
-            # NRMS
-            Q = nrms(sol)
-            # Update best orbit
-            if Q < best_Q
-                best_sol = sol
-                best_Q = Q
-                # Break condition
-                if Q <= params.gaussQmax
-                    flag = true
-                    break
-                end
-            end
-        end
-        flag && break
+    # Find best triplet of observations
+    _gausstriplet!(observatories, dates, α, δ, tracklets)
+    # Julian day of middle observation
+    _jd0_ = dtutc2jdtdb(dates[2])
+    # Gauss method solution
+    solG = gauss_method(observatories, dates, α .+ dq[1:3], δ .+ dq[4:6], params)
+    # Filter non-physical (negative) rho solutions
+    filter!(x -> all(cte.(x.ρ) .> 0), solG)
+    # Filter Gauss solutions by heliocentric energy
+    filter!(x -> cte(cte(heliocentric_energy(x.statevect))) <= 0, solG)
+    # Iterate over Gauss solutions
+    for i in eachindex(solG)
+        # Epoch (corrected for light-time)
+        jd0 = _jd0_ - cte(solG[i].ρ[2]) / c_au_per_day
+        # 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)
+        # Update solution
+        sol = updatesol(sol, _sol_, radec)
+        # Termination condition
+        nrms(sol) <= params.gaussQmax && return sol
+        # ADAM help
+        tracklets[2].nobs < 2 && continue
+        jd0 = _adam!(q, jd0, tracklets[2], params; dynamics)
+        # Jet Transport Least Squares
+        _sol_ = jtls(radec, tracklets, jd0, q, 0, params, true; dynamics)
+        # Update solution
+        sol = updatesol(sol, _sol_, radec)
+        # Termination condition
+        nrms(sol) <= params.gaussQmax && return sol
     end
 
-    return best_sol
+    return sol
 end
\ No newline at end of file
diff --git a/src/orbitdetermination/neosolution.jl b/src/orbitdetermination/neosolution.jl
index 1a01cfc7..1c80e5ad 100644
--- a/src/orbitdetermination/neosolution.jl
+++ b/src/orbitdetermination/neosolution.jl
@@ -152,6 +152,12 @@ end
 
 iszero(x::NEOSolution{T, U}) where {T <: Real, U <: Number} = x == zero(NEOSolution{T, U})
 
+# Override Base.min
+function min(x::NEOSolution{T, T}, y::NEOSolution{T, T}) where {T <: Real}
+    nrms(x) <= nrms(y) && return x
+    return y
+end
+
 # Normalized Root Mean Square Error
 function nrms(sol::NEOSolution{T, T}) where {T <: Real}
     if iszero(sol)
@@ -190,6 +196,16 @@ Return `sol`'s initial condition signal-to-noise ratios in barycentric cartesian
 """
 snr(sol::NEOSolution{T, T}) where {T <: Real} = abs.(sol()) ./ sigmas(sol)
 
+@doc raw"""
+    minmaxdates(sol::NEOSolution{T, U}) where {T <: Real, U <: Number}
+
+Return the dates of the earliest and latest observation in `sol`.
+"""
+function minmaxdates(sol::NEOSolution{T, U}) where {T <: Real, U <: Number}
+    dates = map(t -> extrema(date, t.radec), sol.tracklets)
+    return minimum(first, dates), maximum(last, dates)
+end
+
 @doc raw"""
     jplcompare(des::String, sol::NEOSolution{T, U}) where {T <: Real, U <: Number}
 
@@ -202,7 +218,7 @@ function jplcompare(des::String, sol::NEOSolution{T, U}) where {T <: Real, U <:
     # NEOs barycentric state vector
     q1 = cte.(sol())
     # Time of first (last) observation
-    t0, tf = date(sol.tracklets[1].radec[1]), date(sol.tracklets[end].radec[end])
+    t0, tf = minmaxdates(sol)
     # Download JPL ephemerides
     bsp = smb_spk("DES = $(des);", t0, tf)
     # Object ID
diff --git a/src/orbitdetermination/orbitdetermination.jl b/src/orbitdetermination/orbitdetermination.jl
index 4da3f118..f282f70b 100644
--- a/src/orbitdetermination/orbitdetermination.jl
+++ b/src/orbitdetermination/orbitdetermination.jl
@@ -6,7 +6,7 @@ include("admissibleregion.jl")
 # Times used within propres
 function _proprestimes(radec::Vector{RadecMPC{T}}, jd0::U, params::NEOParameters{T}) where {T <: Real, U <: Number}
     # Time of first (last) observation
-    t0, tf = datetime2julian(date(radec[1])), datetime2julian(date(radec[end]))
+    t0, tf = dtutc2jdtdb(date(radec[1])), dtutc2jdtdb(date(radec[end]))
     # TDB epoch (plain)
     _jd0_ = cte(cte(jd0))
     # Years in backward (forward) integration
@@ -78,24 +78,94 @@ function propres!(res::Vector{OpticalResidual{T, U}}, radec::Vector{RadecMPC{T}}
     end
 end
 
+# Initial subset of radec for jtls
+function _initialtracklets(tracklets::Vector{Tracklet{T}}, i::Int) where {T <: Real}
+    if iszero(i)
+        tin = deepcopy(tracklets)
+        tout = Vector{Tracklet{T}}(undef, 0)
+        rin = sort!(reduce(vcat, indices.(tracklets)))
+    else
+        tin = [tracklets[i]]
+        tout = sort(tracklets, by = x -> abs( (x.date - tracklets[i].date).value ))[2:end]
+        rin = sort!(indices(tracklets[i]))
+        while length(rin) < 3 && !isempty(tout)
+            tracklet = popfirst!(tout)
+            push!(tin, tracklet)
+            sort!(tin)
+            rin = vcat(rin, tracklet.indices)
+            sort!(rin)
+        end
+    end
+    return tin, tout, rin
+end
+
+# Incrementally add observations to fit
+
+# Add as much tracklets as possible per iteration
+function addradec!(::Val{true}, rin::Vector{Int}, fit::LeastSquaresFit{T},
+    tin::Vector{Tracklet{T}}, tout::Vector{Tracklet{T}},
+    res::Vector{OpticalResidual{T, TaylorN{T}}}, x0::Vector{T},
+    params::NEOParameters{T}) where {T <: Real}
+    while !isempty(tout)
+        extra = indices(tout[1])
+        fit_new = tryls(res[rin ∪ extra], x0, params.newtoniter)
+        !fit_new.success && break
+        fit = fit_new
+        tracklet = popfirst!(tout)
+        push!(tin, tracklet)
+        sort!(tin)
+        rin = vcat(rin, extra)
+        sort!(rin)
+    end
+
+    return rin, fit
+end
+
+# Add at most one tracklet per iteration
+function addradec!(::Val{false}, rin::Vector{Int}, fit::LeastSquaresFit{T},
+    tin::Vector{Tracklet{T}}, tout::Vector{Tracklet{T}},
+    res::Vector{OpticalResidual{T, TaylorN{T}}}, x0::Vector{T},
+    params::NEOParameters{T}) where {T <: Real}
+    if nrms(res[rin], fit) < params.tsaQmax && !isempty(tout)
+        extra = indices(tout[1])
+        fit_new = tryls(res[rin ∪ extra], x0, params.newtoniter)
+        !fit_new.success && return rin, fit
+        fit = fit_new
+        tracklet = popfirst!(tout)
+        push!(tin, tracklet)
+        sort!(tin)
+        rin = vcat(rin, extra)
+        sort!(rin)
+    end
+
+    return rin, fit
+end
+
 @doc raw"""
-    jtls(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, jd0::V, q::Vector{U},
-         g0::Int, gf::Int, params::NEOParameters{T}; dynamics::D = newtonian!) where {D, T <: Real, U <: Number, V <: Number}
+    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}
 
 Compute an orbit via Jet Transport Least Squares.
 
-# Arguments
+## Arguments
 
 - `radec::Vector{RadecMPC{T}}`: vector of optical astrometry.
 - `tracklets::Vector{Tracklet{T}},`: vector of tracklets.
 - `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).
-- `dynamics::D`: dynamical model.
+- `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}; dynamics::D = newtonian!) where {D, T <: Real, V <: Number}
+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}
     # Plain initial condition
     q0 = constant_term.(q)
     # JT tail
@@ -103,15 +173,13 @@ function jtls(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, jd0::V
     # Vector of O-C residuals
     res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec))
     # Propagation buffer
-    t0, tf = datetime2days(date(radec[1])), datetime2days(date(radec[end]))
+    t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end]))
     tlim = (t0 - params.bwdoffset, tf + params.fwdoffset)
     buffer = PropagationBuffer(dynamics, jd0, tlim, q, params)
     # Origin
     x0 = zeros(T, 6)
-    # Subset of radec for orbit fit
-    tin = [tracklets[i]]
-    tout = sort(tracklets, by = x -> abs( (x.date - tracklets[i].date).value ))[2:end]
-    rin = sort!(indices(tracklets[i]))
+    # Initial subset of radec for orbit fit
+    tin, tout, rin = _initialtracklets(tracklets, i)
     # Residuals space to barycentric coordinates jacobian
     J = Matrix(TS.jacobian(dq))
     # Best orbit
@@ -130,20 +198,10 @@ function jtls(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, jd0::V
         fit = tryls(res[rin], x0, params.newtoniter)
         !fit.success && break
         # Incrementally add observations to fit
-        while !isempty(tout)
-            extra = indices(tout[1])
-            fit_new = tryls(res[rin ∪ extra], x0, params.newtoniter)
-            !fit_new.success && break
-            fit = fit_new
-            tracklet = popfirst!(tout)
-            push!(tin, tracklet)
-            sort!(tin)
-            rin = vcat(rin, extra)
-            sort!(rin)
-        end
+        rin, fit = addradec!(Val(mode), rin, fit, tin, tout, res, x0, params)
         # NRMS
         Q = nrms(res, fit)
-        if length(rin) == length(radec) && abs(best_Q - Q) < 0.1
+        if abs(best_Q - Q) < 0.1
             flag = true
         end
         # Update NRMS and initial conditions
@@ -164,6 +222,27 @@ function jtls(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}}, jd0::V
     end
 end
 
+# Default naive initial conditions for iod
+function iodinitcond(A::AdmissibleRegion{T}) where {T <: Real}
+    v_ρ = sum(A.v_ρ_domain) / 2
+    return [
+        (A.ρ_domain[1], v_ρ, :log),
+        (10^(sum(log10, A.ρ_domain) / 2), v_ρ, :log),
+        (sum(A.ρ_domain) / 2, v_ρ, :log),
+        (A.ρ_domain[2], v_ρ, :log),
+    ]
+end
+
+# Update `sol` iff `_sol_` is complete and has a lower nrms
+function updatesol(sol::NEOSolution{T, T}, _sol_::NEOSolution{T, T},
+    radec::Vector{RadecMPC{T}}) where {T <: Real}
+    if length(_sol_.res) == length(radec)
+        return min(sol, _sol_)
+    else
+        return sol
+    end
+end
+
 include("tooshortarc.jl")
 include("gaussinitcond.jl")
 
@@ -177,87 +256,76 @@ function issinglearc(radec::Vector{RadecMPC{T}}, arc::Day = Day(30)) where {T <:
     return all(diff(date.(radec)) .< arc)
 end
 
-@doc raw"""
-    isgauss(sol::NEOSolution{T, T}) where {T <: Real}
-
-Check whether `sol` was computed via [`gaussinitcond`](@ref) (`true`) or
-via [`tooshortarc`](@ref) (`false`).
-"""
-function isgauss(tracklets::Vector{Tracklet{T}}) where {T <: Real}
-    # Observing stations
-    obs = observatory.(tracklets)
-    # TSA is not well suited for satellite observatories
-    any(issatellite.(obs)) && return true
-    # Gauss cannot handle less than 3 tracklets
-    length(tracklets) < 3 && return false
-    # Time span
-    Δ = numberofdays(tracklets)
-    # Gauss approximation does not work with less than 1 day
-    return Δ > 1
-end
-
-isgauss(sol::NEOSolution{T, T}) where {T <: Real} = isgauss(sol.tracklets)
-
 @doc raw"""
     orbitdetermination(radec::Vector{RadecMPC{T}}, params::NEOParameters{T};
-                       dynamics::D = newtonian!) where {T <: Real, D}
+        kwargs...) where {T <: Real, D1, D2}
 
 Initial Orbit Determination (IOD) routine.
 
-# Arguments
+## Arguments
 
 - `radec::Vector{RadecMPC{T}}`: vector of optical astrometry.
 - `params::NEOParameters{T}`: see [`NEOParameters`](@ref).
-- `dynamics::D`: dynamical model.
+
+## 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}}`,
+    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::D = newtonian!) where {T <: Real, D}
+    dynamics::D1 = newtonian!, initcond::D2 = iodinitcond) where {T <: Real, D1, D2}
     # Allocate memory for orbit
     sol = zero(NEOSolution{T, T})
     # Eliminate observatories without coordinates
     filter!(x -> hascoord(observatory(x)), radec)
     # Cannot handle zero observations or multiple arcs
-    if iszero(length(radec)) || !issinglearc(radec)
-        return sol
-    end
+    (isempty(radec) || !issinglearc(radec)) && return sol
     # Reduce tracklets by polynomial regression
     tracklets = reduce_tracklets(radec)
     # Set jet transport variables
-    varorder = max(params.tsaorder, params.gaussorder)
+    varorder = max(params.tsaorder, params.gaussorder, params.jtlsorder)
     scaled_variables("dx", ones(T, 6); order = varorder)
-    # Case 1: Gauss' Method
-    if isgauss(tracklets)
-        sol = gaussinitcond(radec, tracklets, params; dynamics)
-    end
-    # Case 2: Too Short Arc (TSA)
-    if iszero(sol) || nrms(sol) > params.gaussQmax
-        sol = tooshortarc(radec, tracklets, params; dynamics)
-    end
+    # Gauss method
+    _sol_ = gaussiod(radec, tracklets, params; dynamics)
+    # Update solution
+    sol = updatesol(sol, _sol_, radec)
+    # Termination condition
+    nrms(sol) <= params.gaussQmax && return sol
+    # Too short arc
+    _sol_ = tsaiod(radec, tracklets, params; dynamics, initcond)
+    # Update solution
+    sol = updatesol(sol, _sol_, radec)
 
     return sol
 end
 
 @doc raw"""
-    orbitdetermination(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T}, params::NEOParameters{T};
-                       dynamics::D = newtonian!) where {T <: Real, D}
+    orbitdetermination(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T},
+        params::NEOParameters{T}; kwargs...) where {T <: Real, D}
 
-Fit `sol` to `radec` via Jet Transport Least Squares.
+Fit a least squares orbit to `radec` using `sol` as an initial condition.
 
-# Arguments
+## Arguments
 
 - `radec::Vector{RadecMPC{T}}`: vector of optical astrometry.
 - `sol::NEOSolution{T, T}:` preliminary orbit.
 - `params::NEOParameters{T}`: see [`NEOParameters`](@ref).
-- `dynamics::D`: dynamical model.
+
+## 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}
+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)
     # Reference epoch [Julian days TDB]
@@ -271,6 +339,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(datetime2days(tracklet.date) - sol.bwd.t0), tracklets)
+    _, i = findmin(tracklet -> abs(dtutc2days(tracklet.date) - sol.bwd.t0), tracklets)
     return jtls(radec, tracklets, jd0, q, i, params; dynamics)
 end
\ No newline at end of file
diff --git a/src/orbitdetermination/tooshortarc.jl b/src/orbitdetermination/tooshortarc.jl
index a3e35ab6..5dd140b8 100644
--- a/src/orbitdetermination/tooshortarc.jl
+++ b/src/orbitdetermination/tooshortarc.jl
@@ -13,7 +13,6 @@ from `(ρ, v_ρ)`.
 - `μ::T`: first moment (default: `0.75`).
 - `ν::T`: second moment (default: `0.9`).
 - `ϵ::T`: numerical stability constant (default: `1e-8`).
-- `Qtol::T`: target function relative tolerance (default: `1e-5`).
 - `adamorder::Int`: jet transport order (default: `2`).
 - `dynamics::D`: dynamical model (default: `newtonian!`).
 
@@ -22,12 +21,14 @@ from `(ρ, v_ρ)`.
 """
 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, Qtol::T = 0.001, adamorder::Int = 2,
+              μ::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]
-    jd0 = datetime2julian(A.date)
+    jd0 = dtutc2jdtdb(A.date)
     # Maximum number of iterations
     maxiter = params.adamiter
+    # Target function relative tolerance
+    Qtol = params.adamQtol
     # Allocate memory
     aes = Matrix{T}(undef, 6, maxiter+1)
     Qs = fill(T(Inf), maxiter+1)
@@ -45,7 +46,7 @@ function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T
     # Jet transport variables
     dae = [scalings[i] * TaylorN(i, order = adamorder) for i in 1:6]
     # Propagation buffer
-    t0, tf = datetime2days(date(radec[1])), datetime2days(date(radec[end]))
+    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)
     # Vector of O-C residuals
@@ -110,115 +111,67 @@ function adam(radec::Vector{RadecMPC{T}}, A::AdmissibleRegion{T}, ρ::T, v_ρ::T
     return aes[:, t], Qs[t]
 end
 
-# Order in which to check tracklets in tooshortarc
-function tsatrackletorder(x::Tracklet{T}, y::Tracklet{T}) where {T <: Real}
-    if x.nobs == y.nobs
-        return x.date > y.date
-    else
-        return x.nobs > y.nobs
-    end
-end
+@doc raw"""
+    tsaiod(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}},
+        params::NEOParameters{T}; kwargs...) where {T <: Real, D1, D2}
 
-# Point in manifold of variations -> NEOSolution{T, T}
-function _tooshortarc(A::AdmissibleRegion{T}, radec::Vector{RadecMPC{T}},
-                      tracklets::Vector{Tracklet{T}}, i::Int, params::NEOParameters{T};
-                      scale::Symbol = :log, dynamics::D = newtonian!) where {T <: Real, D}
-    # Center
-    if scale == :linear
-        ρ = sum(A.ρ_domain) / 2
-    elseif scale == :log
-        ρ = A.ρ_domain[1]
-    end
-    v_ρ = sum(A.v_ρ_domain) / 2
-    # ADAM minimization over manifold of variations
-    ae, Q = adam(tracklets[i].radec, A, ρ, v_ρ, params; scale, dynamics)
-    # ADAM failed to converge
-    if isinf(Q)
-        return zero(NEOSolution{T, T})
-    # Jet Transport Least Squares
-    else
-        # Initial time of integration [julian days]
-        # (corrected for light-time)
-        jd0 = datetime2julian(A.date) - ae[5] / c_au_per_day
-        # Convert attributable elements to barycentric cartesian coordinates
-        q0 = attr2bary(A, ae, params)
-        # Scaling factors
-        scalings = abs.(q0) ./ 10^5
-        # Jet Transport initial condition
-        q = [q0[i] + scalings[i] * TaylorN(i, order = params.tsaorder) for i in 1:6]
-        # Jet Transport Least Squares
-        return jtls(radec, tracklets, jd0, q, i, params; dynamics)
-    end
-end
 
-@doc raw"""
-    tooshortarc(radec::Vector{RadecMPC{T}}, [tracklets::Vector{Tracklet{T}},]
-                params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D}
 
 Compute an orbit by minimizing the normalized root mean square residual
 over the manifold of variations.
 
-# Arguments
+## Arguments
 
 - `radec::Vector{RadecMPC{T}}`: vector of optical astrometry.
 - `tracklets::Vector{Tracklet{T}},`: vector of tracklets.
 - `params::NEOParameters{T}`: see `Too Short Arc Parameters` of [`NEOParameters`](@ref).
-- `dynamics::D`: dynamical model.
-
-!!! warning
-    This function will change the (global) `TaylorSeries` variables.
-"""
-function tooshortarc(radec::Vector{RadecMPC{T}}, params::NEOParameters{T};
-                     dynamics::D = newtonian!) where {T <: Real, D}
-    # Reduce tracklets by polynomial regression
-    tracklets = reduce_tracklets(radec)
-    # Set jet transport variables
-    varorder = max(params.tsaorder, params.gaussorder)
-    scaled_variables("dx", ones(T, 6); order = varorder)
-    # Too Short Arc (TSA)
-    return tooshortarc(radec, tracklets, params; dynamics)
-end
 
-function tooshortarc(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}},
-                     params::NEOParameters{T}; dynamics::D = newtonian!) where {T <: Real, D}
-
-    # Allocate memory for output
-    best_sol = zero(NEOSolution{T, T})
-    # Sort tracklets by tsatrackletorder
-    idxs = sortperm(tracklets, lt = tsatrackletorder)
+## 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}}`,
+    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}
+    # Allocate memory for orbit
+    sol = zero(NEOSolution{T, T})
     # Iterate tracklets
-    for i in idxs
+    for i in eachindex(tracklets)
+        # ADAM requires a minimum of 2 observations
+        tracklets[i].nobs < 2 && continue
         # Admissible region
         A = AdmissibleRegion(tracklets[i], params)
-        iszero(A) && continue
-        # See Table 1 of https://doi.org/10.1051/0004-6361/201732104
-        if A.ρ_domain[2] < sqrt(10)
-            sol1 = _tooshortarc(A, radec, tracklets, i, params;
-                                scale = :log, dynamics)
-            # Break condition
-            nrms(sol1) < params.tsaQmax && return sol1
-            sol2 = _tooshortarc(A, radec, tracklets, i, params;
-                                scale = :linear, dynamics)
-            # Break condition
-            nrms(sol2) < params.tsaQmax && return sol2
-        else
-            sol1 = _tooshortarc(A, radec, tracklets, i, params;
-                                scale = :linear, dynamics)
-            # Break condition
-            nrms(sol1) < params.tsaQmax && return sol1
-            sol2 = _tooshortarc(A, radec, tracklets, i, params;
-                                scale = :log, dynamics)
-            # Break condition
-            nrms(sol2) < params.tsaQmax && return sol2
-        end
-        # Update best solution
-        if nrms(sol1) < nrms(best_sol)
-            best_sol = sol1
-        elseif nrms(sol2) < nrms(best_sol)
-            best_sol = sol2
+        # 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)
+            # ADAM failed to converge
+            isinf(Q) && continue
+            # Initial time of integration [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)
+            # Scaling factors
+            scalings = abs.(q0) ./ 10^5
+            # 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)
+            # Update solution
+            sol = updatesol(sol, _sol_, radec)
+            # Termination condition
+            nrms(sol) <= params.tsaQmax && return sol
         end
     end
 
-    return best_sol
+    return sol
 end
\ No newline at end of file
diff --git a/src/propagation/parameters.jl b/src/propagation/parameters.jl
index b9ef46f2..4665f4f9 100644
--- a/src/propagation/parameters.jl
+++ b/src/propagation/parameters.jl
@@ -24,6 +24,7 @@ struct NEOParameters{T <: Real}
     H_max::T
     a_max::T
     adamiter::Int
+    adamQtol::T
     tsaorder::Int
     tsaQmax::T
     # Jet Transport Least Squares Parameters
@@ -41,8 +42,8 @@ function show(io::IO, p::NEOParameters{T}) where {T <: Real}
         :maxsteps, :order, :abstol, :parse_eqs, :bwdoffset,
         :fwdoffset, :coeffstol, :max_triplets, :gaussorder,
         :adamhelp, :gaussQmax, :H_max, :a_max, :adamiter,
-        :tsaorder, :tsaQmax, :jtlsiter, :newtoniter, :jtlsorder,
-        :max_per
+        :adamQtol, :tsaorder, :tsaQmax, :jtlsiter, :newtoniter,
+        :jtlsorder, :max_per
     ]
     s = Vector{String}(undef, length(params)+1)
     s[1] = "NEOParameters{$T}:\n"
@@ -90,6 +91,7 @@ Parameters for all orbit determination functions.
 - `H_max::T`: maximum absolute magnitude (default: `34.5`).
 - `a_max::T`: maximum semimajor axis (default: `100.0`).
 - `adamiter::Int`: maximum number of iterations for `ADAM` optimizer (default: `200`).
+- `Qtol::T`: target function relative tolerance (default: `0.001`).
 - `tsaorder::Int`: order of the jet transport perturbation (default: `6`).
 - `tsaQmax::T`: nrms threshold (default: `1.5`).
 
@@ -108,9 +110,9 @@ function NEOParameters(;
     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, tsaorder::Int = 6, tsaQmax::T = 1.5,
-    jtlsiter::Int = 5, newtoniter::Int = 5, jtlsorder::Int = 5, max_per::T = 18.0
-    ) where {T <: Real}
+    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)
@@ -124,8 +126,8 @@ function NEOParameters(;
         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, tsaorder, tsaQmax, jtlsiter, newtoniter, jtlsorder,
-        max_per
+        a_max, adamiter, adamQtol, tsaorder, tsaQmax, jtlsiter, newtoniter,
+        jtlsorder, max_per
     )
 end
 
diff --git a/test/bplane.jl b/test/bplane.jl
index 3b53dbaf..4c1f7e54 100644
--- a/test/bplane.jl
+++ b/test/bplane.jl
@@ -11,13 +11,15 @@ using NEOs: propres
 @testset "2018 LA" begin
 
     # Fetch optical astrometry
-    radec = fetch_radec_mpc("designation" => "2018 LA")
+    radec = fetch_radec_mpc("2018 LA")
     # Parameters
     params = NEOParameters(coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007)
 
-    # Orbit Determination
+    # Initial Orbit Determination
     sol = orbitdetermination(radec[1:8], params)
+    # Update parameters
     params = NEOParameters(params; jtlsorder = 6)
+    # Refine orbit
     sol = orbitdetermination(radec, sol, params)
 
     # Radial velocity with respect to the Earth.
diff --git a/test/observations.jl b/test/observations.jl
index 6d47484a..f027ca64 100644
--- a/test/observations.jl
+++ b/test/observations.jl
@@ -160,7 +160,7 @@ using NEOs: src_path
 
         # Get RadecMPC
         source_file = joinpath(pkgdir(NEOs), "test", "data", "99942.txt")
-        get_radec_mpc("number" => "99942", source_file)
+        get_radec_mpc("99942", source_file)
 
         @test isfile(source_file)
 
@@ -232,7 +232,7 @@ using NEOs: src_path
         @test 0 <= object.α <= 2π
         @test -π <= object.δ <= π
         @test object.nobs > 0
-        @test object.arc > 0
+        @test object.arc >= 0
         @test object.notseen > 0
         # Check optical astrometry
         filename = get_radec_neocp(object.tmpdesig)
@@ -240,7 +240,7 @@ using NEOs: src_path
         radec1 = read_radec_mpc(filename)
         radec2 = fetch_radec_neocp(object.tmpdesig)
         rm(filename)
-        @test length(radec1) == length(radec2) == object.nobs
+        @test length(radec1) == length(radec2) >= object.nobs
         @test radec1 == radec2
         notseen = (now(UTC) - date(radec1[end])).value / 86_400_000
         @test round(notseen, RoundUp, digits = 3) >= object.notseen
@@ -250,6 +250,12 @@ using NEOs: src_path
         rm(filename)
     end
 
+    @testset "Time conversions" begin
+        t0 = now()
+        @test et2dtutc(dtutc2et(t0)) == t0
+        @test jdtdb2dtutc(dtutc2jdtdb(t0)) == t0
+    end
+
     @testset "Topocentric" begin
         using NEOs: TimeOfDay, sunriseset, obsposECEF, obsposvelECI
 
@@ -296,8 +302,8 @@ using NEOs: src_path
         # Sunrise and sunset
         radec = read_radec_mpc("99942        8C2020 12 08.15001011 20 07.510-08 02 54.20         18.50GV~4ROF094")
         sun = sunriseset(radec[1])
-        @test datetime2julian(sun[1]) ≈ datetime2julian(DateTime("2020-12-08T05:05:59.384"))
-        @test datetime2julian(sun[2]) ≈ datetime2julian(DateTime("2020-12-08T14:05:49.386"))
+        @test dtutc2jdtdb(sun[1]) ≈ dtutc2jdtdb(DateTime("2020-12-08T05:05:59.384"))
+        @test dtutc2jdtdb(sun[2]) ≈ dtutc2jdtdb(DateTime("2020-12-08T14:05:49.386"))
 
         # obsposECEF
         ecef_2 = obsposECEF.(radec_2)
@@ -321,7 +327,7 @@ using NEOs: src_path
         # Choose this example because of the discontinuity in α
 
         # Fetch optical astrometry
-        radec = fetch_radec_mpc("designation" => "2020 TJ6")
+        radec = fetch_radec_mpc("2020 TJ6")
         # Reduce tracklets
         tracklets = reduce_tracklets(radec)
 
diff --git a/test/orbitdetermination.jl b/test/orbitdetermination.jl
index 35f1cf1c..f74e95da 100644
--- a/test/orbitdetermination.jl
+++ b/test/orbitdetermination.jl
@@ -6,104 +6,127 @@ using PlanetaryEphemeris
 using LinearAlgebra
 using Test
 
-using NEOs: NEOSolution, numberofdays
+using NEOs: NEOSolution, RadecMPC, reduce_tracklets,
+    indices, numberofdays
+
+function iodsubradec(radec::Vector{RadecMPC{T}}, N::Int = 3) where {T <: Real}
+    tracklets = reduce_tracklets(radec)
+    idxs = sort!(reduce(vcat, indices.(tracklets[1:N])))
+    subradec = radec[idxs]
+    return subradec
+end
 
 @testset "Orbit Determination" begin
     @testset "Gauss Method (without ADAM)" begin
         # Load observations
         radec = read_radec_mpc(joinpath(pkgdir(NEOs), "test", "data", "RADEC_2023_DW.dat"))
+        # Subset of radec for IOD
+        subradec = iodsubradec(radec, 3)
+
+        @test length(subradec) == 9
+        @test numberofdays(subradec) < 0.18
+
         # Parameters
         params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007, parse_eqs = false)
-        params = NEOParameters(params, parse_eqs=true)
+        params = NEOParameters(params, parse_eqs = true)
 
-        # Orbit Determination
-        sol = orbitdetermination(radec, params)
+        # Initial Orbit Determination
+        sol = orbitdetermination(subradec, params)
 
-        # Values by May 31, 2024
+        # Values by Sep 25, 2024
 
-        # Vector of observations
-        @test length(radec) == 123
-        @test numberofdays(radec) < 21.0
         # Orbit solution
         @test isa(sol, NEOSolution{Float64, Float64})
         # Tracklets
-        @test length(sol.tracklets) == 44
-        @test sol.tracklets[1].radec[1] == radec[1]
-        @test sol.tracklets[end].radec[end] == radec[end]
+        @test length(sol.tracklets) == 3
+        @test sol.tracklets[1].radec[1] == subradec[1]
+        @test sol.tracklets[end].radec[end] == subradec[end]
         @test issorted(sol.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test dtutc2days(date(subradec[1])) > sol.bwd.t0 + sol.bwd.t[end]
         @test all( norm.(sol.bwd.x, Inf) .< 2 )
         @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test dtutc2days(date(subradec[end])) < sol.fwd.t0 + sol.fwd.t[end]
         @test all( norm.(sol.fwd.x, Inf) .< 2 )
         @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
         # Vector of residuals
-        @test length(sol.res) == 123
+        @test length(sol.res) == 9
         @test iszero(count(outlier.(sol.res)))
         # Least squares fit
         @test sol.fit.success
-        @test all( sigmas(sol) .< 5e-5 )
-        @test all( snr(sol) .> 4_500)
-        @test nrms(sol) < 0.36
+        @test all( sigmas(sol) .< 9e-4 )
+        @test all( snr(sol) .> 14.5)
+        @test nrms(sol) < 0.18
         # Jacobian
         @test size(sol.jacobian) == (6, 6)
         @test !isdiag(sol.jacobian)
-        @test maximum(sol.jacobian) < 8e-4
+        @test maximum(sol.jacobian) < 3e-4
         # Compatibility with JPL
-        JPL = [-1.1003236797145037, 0.20774505704014837, 0.04203643429323372,
-               -0.004736048200346307, -0.010626587751050683, -0.006016238714906758]
-        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.27)
+        JPL = [-0.9867704701732631, 0.3781890325424674, 0.14094513213009532,
+            -0.008773157203087259, -0.00947109649687576, -0.005654229864757284]
+        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.76)
+
+        # Add observations
+        subradec = iodsubradec(radec, 15)
+
+        @test length(subradec) == 43
+        @test numberofdays(subradec) < 2.76
 
         # Refine orbit
-        sol1 = orbitdetermination(radec, sol, params)
+        sol1 = orbitdetermination(subradec, sol, params)
 
         # Orbit solution
         @test isa(sol1, NEOSolution{Float64, Float64})
         # Tracklets
-        @test sol1.tracklets == sol.tracklets
+        @test length(sol1.tracklets) == 15
+        @test sol1.tracklets[1].radec[1] == subradec[1]
+        @test sol1.tracklets[end].radec[end] == subradec[end]
+        @test issorted(sol1.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[1])) > sol1.bwd.t0 + sol1.bwd.t[end]
+        @test dtutc2days(date(subradec[1])) > sol1.bwd.t0 + sol1.bwd.t[end]
         @test all( norm.(sol1.bwd.x, Inf) .< 1.2 )
         @test isempty(sol1.t_bwd) && isempty(sol1.x_bwd) && isempty(sol1.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[end])) < sol1.fwd.t0 + sol1.fwd.t[end]
+        @test dtutc2days(date(subradec[end])) < sol1.fwd.t0 + sol1.fwd.t[end]
         @test all( norm.(sol1.fwd.x, Inf) .< 1.2 )
         @test isempty(sol1.t_fwd) && isempty(sol1.x_fwd) && isempty(sol1.g_fwd)
         # Vector of residuals
-        @test length(sol1.res) == 123
+        @test length(sol1.res) == 43
         @test iszero(count(outlier.(sol1.res)))
         # Least squares fit
         @test sol1.fit.success
-        @test all( sigmas(sol1) .< 5e-5 )
-        @test all( snr(sol1) .> 4_500 )
-        @test nrms(sol1) < 0.36
+        @test all( sigmas(sol1) .< 2e-4 )
+        @test all( snr(sol1) .> 866 )
+        @test nrms(sol1) < 0.37
         # Jacobian
         @test size(sol1.jacobian) == (6, 6)
         @test isdiag(sol1.jacobian)
-        @test maximum(sol1.jacobian) < 2e-6
+        @test maximum(sol1.jacobian) < 1e-6
         # Compatibility with JPL
-        @test all(abs.(sol1() - JPL) ./ sigmas(sol1) .< 0.27)
+        @test all(abs.(sol1() - JPL) ./ sigmas(sol1) .< 0.31)
         # MPC Uncertainty Parameter
         @test uncertaintyparameter(radec, sol1, params) == 6
     end
 
     @testset "Gauss Method (with ADAM)" begin
         # Load observations
-        radec = fetch_radec_mpc("designation" => "2016 TU93")
+        radec = fetch_radec_mpc("2016 TU93")
+        # Subset of radec for IOD
+        # subradec = iodsubradec(radec, 3)
+        # In this case: radec == subradec
+
+        @test length(radec) == 9
+        @test numberofdays(radec) < 13.1
 
         # Parameters
         params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007, adamhelp = true)
 
-        # Orbit Determination
+        # Initial Orbit Determination
         sol = orbitdetermination(radec, params)
 
-        # Values by June 1, 2024
+        # Values by Sep 25, 2024
 
-        # Vector of observations
-        @test length(radec) == 9
-        @test numberofdays(radec) < 13.1
         # Orbit solution
         @test isa(sol, NEOSolution{Float64, Float64})
         # Tracklets
@@ -112,11 +135,11 @@ using NEOs: NEOSolution, numberofdays
         @test sol.tracklets[end].radec[end] == radec[end]
         @test issorted(sol.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test dtutc2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
         @test all( norm.(sol.bwd.x, Inf) .< 2 )
         @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test dtutc2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
         @test all( norm.(sol.fwd.x, Inf) .< 2 )
         @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
         # Vector of residuals
@@ -132,20 +155,20 @@ using NEOs: NEOSolution, numberofdays
         @test isdiag(sol.jacobian)
         @test maximum(sol.jacobian) < 1e-5
         # Compatibility with JPL
-        JPL = [1.0102564188486982, 0.2934743828145318, 0.10467187893161536,
-               -0.0002634434601757652, 0.01837381321202214, 0.007208485181422459]
+        JPL = [1.010256210875638, 0.29348888228332376, 0.10467756741992229,
+            -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
     end
 
     @testset "Admissible region" begin
-        using NEOs: AdmissibleRegion, reduce_tracklets, arenergydis, rangerate,
-                    rangerates, argoldensearch, arboundary, _helmaxrange, R_SI,
-                    k_gauss, μ_ES, topo2bary, bary2topo
+        using NEOs: AdmissibleRegion, arenergydis, rangerate, rangerates,
+            argoldensearch, arboundary, _helmaxrange, R_SI, k_gauss, μ_ES,
+            topo2bary, bary2topo
 
         # Fetch optical astrometry
-        radec = fetch_radec_mpc("designation" => "2024 BX1")
+        radec = fetch_radec_mpc("2024 BX1")
         # Parameters
         params = NEOParameters()
         # First tracklet
@@ -154,12 +177,13 @@ using NEOs: NEOSolution, numberofdays
         # Admissible region
         A = AdmissibleRegion(tracklet, params)
 
-        # Values by August 19, 2024
+        # Values by Sep 25, 2024
 
         # Zero AdmissibleRegion
         @test iszero(zero(AdmissibleRegion{Float64}))
         # Custom print
-        @test string(A) == "AE: [116.61547, 45.39840, -3.21667, 5.76667] t: 2024-01-20T21:50:15.360 obs: GINOP-KHK, Piszkesteto"
+        @test string(A) == "AE: [116.61547, 45.39840, -3.21667, 5.76667] \
+            t: 2024-01-20T21:50:15.360 obs: GINOP-KHK, Piszkesteto"
         # Coefficients
         @test length(A.coeffs) == 6
         @test A.coeffs[3] == A.v_α^2 * cos(A.δ)^2 + A.v_δ^2  # proper motion squared
@@ -181,7 +205,8 @@ using NEOs: NEOSolution, numberofdays
         @test isempty(rangerates(A, ρ0 + 1.0, :inner))
         @test rangerate(A, A.ρ_domain[1], :min, :outer) == A.v_ρ_domain[1]
         @test rangerate(A, A.ρ_domain[1], :max, :outer) == A.v_ρ_domain[2]
-        @test rangerate(A, A.ρ_domain[1], :min, :inner) ≈ -rangerate(A, A.ρ_domain[1], :max, :inner) atol = 1e-18
+        @test rangerate(A, A.ρ_domain[1], :min, :inner) ≈
+            -rangerate(A, A.ρ_domain[1], :max, :inner) atol = 1e-18
         # Golden section search
         ρ, v_ρ = argoldensearch(A, A.ρ_domain..., :min, :outer, 1e-20)
         @test A.ρ_domain[1] ≤ ρ ≤ A.ρ_domain[2]
@@ -246,25 +271,28 @@ using NEOs: NEOSolution, numberofdays
 
     @testset "Too Short Arc" begin
         # Fetch optical astrometry
-        radec = fetch_radec_mpc("designation" => "2008 EK68")
+        radec = fetch_radec_mpc("2008 EK68")
+        # Subset of radec for IOD
+        # subradec = iodsubradec(radec, 3)
+        # In this case: radec == subradec
+
+        @test length(radec) == 10
+        @test numberofdays(radec) < 0.05
 
         # Parameters
         params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007)
 
-        # Orbit Determination
+        # Initial Orbit Determination
         sol = orbitdetermination(radec, params)
 
-        # Values by May 31, 2024
+        # Values by Sep 25, 2024
 
-        # Vector of observations
-        @test length(radec) == 10
-        @test numberofdays(radec) < 0.05
         # Curvature
         C, Γ_C = curvature(radec)
         σ_C = sqrt.(diag(Γ_C))
         @test all( abs.(C) ./ σ_C .> 5.5)
         χ2 = C' * inv(Γ_C) * C
-        @test χ2 > 2_517
+        @test χ2 > 2_516
         # Orbit solution
         @test isa(sol, NEOSolution{Float64, Float64})
         # Tracklets
@@ -273,11 +301,11 @@ using NEOs: NEOSolution, numberofdays
         @test sol.tracklets[end].radec[end] == radec[end]
         @test issorted(sol.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test dtutc2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
         @test all( norm.(sol.bwd.x, Inf) .< 2 )
         @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test dtutc2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
         @test all( norm.(sol.fwd.x, Inf) .< 2 )
         @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
         # Vector of residuals
@@ -293,63 +321,98 @@ using NEOs: NEOSolution, numberofdays
         @test isdiag(sol.jacobian)
         @test maximum(sol.jacobian) < 1e-5
         # Compatibility with JPL
-        JPL = [-0.9698333704468041, 0.2403646120906576, 0.10288887497365079,
-               -0.009512521364762891, -0.015325432155116774, -0.008094623535119534]
-        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.015)
+        JPL = [-0.9698405495747651, 0.24035304578776012, 0.10288276585828428,
+            -0.009512301266159554, -0.01532548565855646, -0.00809464581680694]
+        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.012)
         # MPC Uncertainty Parameter
         @test uncertaintyparameter(radec, sol, params) == 11
     end
 
     @testset "Outlier Rejection" begin
         # Fetch optical astrometry
-        radec = fetch_radec_mpc("designation" => "2007 VV7")
+        radec = fetch_radec_mpc("2007 VV7")
+        # Subset of radec for IOD
+        subradec = iodsubradec(radec, 3)
+
+        @test length(subradec) == 18
+        @test numberofdays(subradec) < 2.16
 
         # Parameters
         params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007)
 
-        # Orbit Determination
-        sol = orbitdetermination(radec, params)
-        # Outlier rejection
-        sol = outlier_rejection(radec, sol, params)
+        # Initial Orbit Determination
+        sol = orbitdetermination(subradec, params)
 
-        # Values by May 31, 2024
+        # Values by Sep 25, 2024
 
-        # Vector of observations
-        @test length(radec) == 21
-        @test numberofdays(radec) < 3.03
         # Orbit solution
         @test isa(sol, NEOSolution{Float64, Float64})
         # Tracklets
-        @test length(sol.tracklets) == 4
-        @test sol.tracklets[1].radec[1] == radec[1]
-        @test sol.tracklets[end].radec[end] == radec[end]
+        @test length(sol.tracklets) == 3
+        @test sol.tracklets[1].radec[1] == subradec[1]
+        @test sol.tracklets[end].radec[end] == subradec[end]
         @test issorted(sol.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test dtutc2days(date(subradec[1])) > sol.bwd.t0 + sol.bwd.t[end]
         @test all( norm.(sol.bwd.x, Inf) .< 2 )
         @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test dtutc2days(date(subradec[end])) < sol.fwd.t0 + sol.fwd.t[end]
         @test all( norm.(sol.fwd.x, Inf) .< 2 )
         @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
         # Vector of residuals
-        @test length(sol.res) == 21
-        @test count(outlier.(sol.res)) == 2
+        @test length(sol.res) == 18
+        @test count(outlier.(sol.res)) == 0
         # Least squares fit
         @test sol.fit.success
         @test all( sigmas(sol) .< 3e-4 )
         @test all( snr(sol) .> 574)
-        @test nrms(sol) < 0.25
+        @test nrms(sol) < 0.46
         # Jacobian
         @test size(sol.jacobian) == (6, 6)
-        @test isdiag(sol.jacobian)
-        @test maximum(sol.jacobian) < 8e-7
+        @test !isdiag(sol.jacobian)
+        @test maximum(sol.jacobian) < 5e-4
+        # Compatibility with JPL
+        JPL = [0.7673366466815864, 0.6484892781853565, 0.29323267343908294,
+            -0.011023343781911974, 0.015392697071667377, 0.006528842022004942]
+        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 1.75)
+
+        # Add remaining observations
+        sol1 = orbitdetermination(radec, sol, params)
+        # Outlier rejection
+        sol1 = outlier_rejection(radec, sol1, params)
+
+        # Orbit solution
+        @test isa(sol1, NEOSolution{Float64, Float64})
+        # Tracklets
+        @test length(sol1.tracklets) == 4
+        @test sol1.tracklets[1].radec[1] == radec[1]
+        @test sol1.tracklets[end].radec[end] == radec[end]
+        @test issorted(sol1.tracklets)
+        # Backward integration
+        @test dtutc2days(date(radec[1])) > sol1.bwd.t0 + sol1.bwd.t[end]
+        @test all( norm.(sol1.bwd.x, Inf) .< 2 )
+        @test isempty(sol1.t_bwd) && isempty(sol1.x_bwd) && isempty(sol1.g_bwd)
+        # Forward integration
+        @test dtutc2days(date(radec[end])) < sol1.fwd.t0 + sol1.fwd.t[end]
+        @test all( norm.(sol1.fwd.x, Inf) .< 2 )
+        @test isempty(sol1.t_fwd) && isempty(sol1.x_fwd) && isempty(sol1.g_fwd)
+        # Vector of residuals
+        @test length(sol1.res) == 21
+        @test count(outlier.(sol1.res)) == 2
+        # Least squares fit
+        @test sol1.fit.success
+        @test all( sigmas(sol1) .< 3e-4 )
+        @test all( snr(sol1) .> 574)
+        @test nrms(sol1) < 0.25
+        # Jacobian
+        @test size(sol1.jacobian) == (6, 6)
+        @test isdiag(sol1.jacobian)
+        @test maximum(sol1.jacobian) < 8e-7
         # Compatibility with JPL
-        JPL = [0.7673449629397204, 0.6484776654615118, 0.2932277478785896,
-               -0.011023192686652665, 0.015392823966811551, 0.0065288994881745974]
-        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 7e-4)
+        @test all(abs.(sol1() - JPL) ./ sigmas(sol1) .< 1.2e-3)
         # MPC Uncertainty Parameter
-        @test uncertaintyparameter(radec, sol, params) == 8
+        @test uncertaintyparameter(radec, sol1, params) == 8
     end
 
     @testset "Interesting NEOs" begin
@@ -357,7 +420,13 @@ using NEOs: NEOSolution, numberofdays
         # 2014 AA hit the Earth around January 2, 2014, 02:49 UTC
 
         # Fetch optical astrometry
-        radec = fetch_radec_mpc("designation" => "2014 AA")
+        radec = fetch_radec_mpc("2014 AA")
+        # Subset of radec for IOD
+        # subradec = iodsubradec(radec, 3)
+        # In this case: radec == subradec
+
+        @test length(radec) == 7
+        @test numberofdays(radec) < 0.05
 
         # Parameters
         params = NEOParameters(coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007)
@@ -365,11 +434,8 @@ using NEOs: NEOSolution, numberofdays
         # Orbit Determination
         sol = orbitdetermination(radec, params)
 
-        # Values by May 31 2024
+        # Values by Sep 25, 2024
 
-        # Vector of observations
-        @test length(radec) == 7
-        @test numberofdays(radec) < 0.05
         # Curvature
         C, Γ_C = curvature(radec)
         σ_C = sqrt.(diag(Γ_C))
@@ -384,11 +450,11 @@ using NEOs: NEOSolution, numberofdays
         @test sol.tracklets[end].radec[end] == radec[end]
         @test issorted(sol.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test dtutc2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
         @test all( norm.(sol.bwd.x, Inf) .< 2 )
         @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test dtutc2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
         @test all( norm.(sol.fwd.x, Inf) .< 1e9 )
         @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
         # Vector of residuals
@@ -404,8 +470,8 @@ using NEOs: NEOSolution, numberofdays
         @test isdiag(sol.jacobian)
         @test maximum(sol.jacobian) < 9e-6
         # Compatibility with JPL
-        JPL = [-0.17932853781716676, 0.8874166708195785, 0.38414497112938667,
-               -0.017557883503263098, -0.005781328976995571, -0.0020073946372627465]
+        JPL = [-0.1793421909678032, 0.8874121750891107, 0.3841434101167349,
+            -0.017557851117612377, -0.005781634223099801, -0.0020075106081869185]
         @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.3)
         # MPC Uncertainty Parameter
         @test uncertaintyparameter(radec, sol, params) == 10
@@ -413,91 +479,176 @@ using NEOs: NEOSolution, numberofdays
         # 2008 TC3 entered the Earth's atmosphere around October 7, 2008, 02:46 UTC
 
         # Fetch optical astrometry
-        radec = fetch_radec_mpc("designation" => "2008 TC3")
+        radec = fetch_radec_mpc("2008 TC3")
+        # Subset of radec for IOD
+        subradec = iodsubradec(radec, 3)
+
+        @test length(subradec) == 18
+        @test numberofdays(subradec) < 0.34
 
         # Parameters
         params = NEOParameters(coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007)
 
-        # Observations with <1" weight
-        idxs = findall(x -> w8sveres17(x) < 1, radec)
-        # Restricted Orbit Determination
-        sol = orbitdetermination(radec[idxs], params)
+        # Initial Orbit Determination
+        sol = orbitdetermination(subradec, params)
 
-        # Values by May 31, 2024
+        # Values by Sep 25, 2024
 
-        # Vector of observations
-        @test length(radec) == 883
-        @test numberofdays(radec) < 0.80
         # Orbit solution
         @test isa(sol, NEOSolution{Float64, Float64})
         # Tracklets
-        @test length(sol.tracklets) == 2
-        @test sol.tracklets[1].radec[1] == radec[idxs[1]]
-        @test sol.tracklets[end].radec[end] == radec[idxs[end]]
+        @test length(sol.tracklets) == 3
+        @test sol.tracklets[1].radec[1] == subradec[1]
+        @test sol.tracklets[end].radec[end] == subradec[end]
         @test issorted(sol.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[idxs[1]])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test dtutc2days(date(subradec[1])) > sol.bwd.t0 + sol.bwd.t[end]
         @test all( norm.(sol.bwd.x, Inf) .< 2 )
         @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[idxs[end]])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test dtutc2days(date(subradec[end])) < sol.fwd.t0 + sol.fwd.t[end]
         @test all( norm.(sol.fwd.x, Inf) .< 1e4 )
         @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
         # Vector of residuals
-        @test length(sol.res) == 20
+        @test length(sol.res) == 18
         @test iszero(count(outlier.(sol.res)))
         # Least squares fit
         @test sol.fit.success
         @test all( sigmas(sol) .< 2e-5 )
-        @test all( snr(sol) .> 732)
-        @test nrms(sol) < 0.30
+        @test all( snr(sol) .> 645)
+        @test nrms(sol) < 0.35
         # Jacobian
         @test size(sol.jacobian) == (6, 6)
-        @test isdiag(sol.jacobian)
-        @test maximum(sol.jacobian) < 1e-5
+        @test !isdiag(sol.jacobian)
+        @test maximum(sol.jacobian) < 4e-6
         # Compatibility with JPL
-        JPL = [0.9741070120439872, 0.2151506132683409, 0.0939089782825125,
-               -0.007890445003489774, 0.016062726197895585, 0.006136042044307594]
-        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.2)
+        JPL = [0.9739760787551061, 0.21541704400792083, 0.09401075290627411,
+            -0.00789675674941779, 0.0160619782715116, 0.006135361409943397]
+        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.20)
 
-        # Update orbit with more observations
-        sol1 = orbitdetermination(radec[1:30], sol, params)
+        # Add observations
+        subradec = iodsubradec(radec, 10)
+
+        @test length(subradec) == 97
+        @test numberofdays(subradec) < 0.70
+
+        # Refine orbit
+        sol1 = orbitdetermination(subradec, sol, params)
 
         # Orbit solution
         @test isa(sol1, NEOSolution{Float64, Float64})
         # Tracklets
-        @test length(sol1.tracklets) == 5
-        @test sol1.tracklets[1].radec[1] == radec[1]
-        @test sol1.tracklets[end].radec[end] == radec[30]
+        @test length(sol1.tracklets) == 10
+        @test sol1.tracklets[1].radec[1] == subradec[1]
+        @test sol1.tracklets[end].radec[end] == subradec[93]
         @test issorted(sol1.tracklets)
         # Backward integration
-        @test datetime2days(date(radec[idxs[1]])) > sol1.bwd.t0 + sol1.bwd.t[end]
+        @test dtutc2days(date(subradec[1])) > sol1.bwd.t0 + sol1.bwd.t[end]
         @test all( norm.(sol1.bwd.x, Inf) .< 1 )
         @test isempty(sol1.t_bwd) && isempty(sol1.x_bwd) && isempty(sol1.g_bwd)
         # Forward integration
-        @test datetime2days(date(radec[idxs[end]])) < sol1.fwd.t0 + sol1.fwd.t[end]
-        @test all( norm.(sol1.fwd.x, Inf) .< 1e4 )
+        @test dtutc2days(date(subradec[end])) < sol1.fwd.t0 + sol1.fwd.t[end]
+        @test all( norm.(sol1.fwd.x, Inf) .< 1e15 )
         @test isempty(sol1.t_fwd) && isempty(sol1.x_fwd) && isempty(sol1.g_fwd)
         # Vector of residuals
-        @test length(sol1.res) == 30
+        @test length(sol1.res) == 97
         @test iszero(count(outlier.(sol1.res)))
         # Least squares fit
         @test sol1.fit.success
-        @test all( sigmas(sol1) .< 2e-6 )
-        @test all( snr(sol1) .> 4_281)
-        @test nrms(sol1) < 0.37
+        @test all( sigmas(sol1) .< 4e-7 )
+        @test all( snr(sol1) .> 21_880)
+        @test nrms(sol1) < 0.53
         # Jacobian
         @test size(sol1.jacobian) == (6, 6)
         @test isdiag(sol1.jacobian)
         @test maximum(sol1.jacobian) < 1e-6
         # Compatibility with JPL
-        @test all(abs.(sol1() - JPL) ./ sigmas(sol1) .< 1.21)
+        @test all(abs.(sol1() - JPL) ./ sigmas(sol1) .< 0.17)
         # Parameters uncertainty
         @test all(sigmas(sol1) .< sigmas(sol))
         # TODO: understand better differences wrt JPL solutions
         # @test nrms(sol1) < nrms(sol)
         # MPC Uncertainty Parameter
-        @test uncertaintyparameter(radec[1:30], sol1, params) == 6
+        @test uncertaintyparameter(subradec, sol1, params) == 5
+    end
+
+    @testset "scripts/orbitdetermination.jl" begin
+        using NEOs: iodinitcond, issatellite
+
+        function radecfilter(radec::Vector{RadecMPC{T}}) where {T <: Real}
+            # Eliminate observations before oficial discovery
+            firstobs = findfirst(r -> !isempty(r.discovery), radec)
+            isnothing(firstobs) && return false, radec
+            radec = radec[firstobs:end]
+            # Filter out incompatible observations
+            filter!(radec) do r
+                hascoord(r.observatory) && !issatellite(r.observatory) &&
+                date(r) >= Date(2000)
+            end
+            length(radec) < 3 && return false, radec
+            # Find the first set of 3 tracklets with a < 15 days timespan
+            tracklets = reduce_tracklets(radec)
+            for i in 1:length(tracklets)-2
+                numberofdays(tracklets[i:i+2]) > 15.0 && continue
+                tracklets = tracklets[i:i+2]
+                radec = reduce(vcat, getfield.(tracklets, :radec))
+                sort!(radec)
+                break
+            end
+            return numberofdays(radec) <= 15.0, radec
+        end
+
+        # Fetch and filter optical astrometry
+        radec = fetch_radec_mpc("2011 UE256")
+        _, radec = radecfilter(radec)
+
+        @test length(radec) == 14
+        @test numberofdays(radec) < 0.85
+
+        # Parameters
+        params = NEOParameters(
+            coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007,
+            gaussorder = 6, gaussQmax = 1.0,
+            adamiter = 500, adamQtol = 1e-5, tsaQmax = 2.0,
+            jtlsiter = 20, newtoniter = 10
+        )
+
+        # Initial Orbit Determination
+        sol = orbitdetermination(radec, params; initcond = iodinitcond)
+
+        # Values by Sep 25, 2024
+
+        # Orbit solution
+        @test isa(sol, NEOSolution{Float64, Float64})
+        # Tracklets
+        @test length(sol.tracklets) == 3
+        @test sol.tracklets[1].radec[1] == radec[1]
+        @test sol.tracklets[end].radec[end] == radec[end]
+        @test issorted(sol.tracklets)
+        # Backward integration
+        @test dtutc2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test all( norm.(sol.bwd.x, Inf) .< 2 )
+        @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
+        # Forward integration
+        @test dtutc2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test all( norm.(sol.fwd.x, Inf) .< 2 )
+        @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
+        # Vector of residuals
+        @test length(sol.res) == 14
+        @test iszero(count(outlier.(sol.res)))
+        # Least squares fit
+        @test sol.fit.success
+        @test all( sigmas(sol) .< 0.24 )
+        @test all( snr(sol) .> 0.40)
+        @test nrms(sol) < 0.42
+        # Jacobian
+        @test size(sol.jacobian) == (6, 6)
+        @test isdiag(sol.jacobian)
+        @test maximum(sol.jacobian) < 1e-5
+        # Compatibility with JPL
+        JPL = [1.5278518651683686, 0.9328047239405421, 0.3755795874791371,
+            -0.014356510695499359, 0.0002883092503924326, 0.002280143873309611]
+        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 1.33)
     end
 
 end
\ No newline at end of file
diff --git a/test/propagation.jl b/test/propagation.jl
index 5ce8f40c..36334b03 100644
--- a/test/propagation.jl
+++ b/test/propagation.jl
@@ -262,7 +262,7 @@ using InteractiveUtils: methodswith
         # Dynamical function
         dynamics = RNp1BP_pN_A_J23E_J2S_eph_threads!
         # Initial date of integration [Julian date TDB]
-        jd0 = datetime2julian(DateTime(2029, 4, 13, 20))
+        jd0 = dtutc2jdtdb(DateTime(2029, 4, 13, 20))
         # Time of integration [years]
         nyears = 0.02
         # Perturbation to nominal initial condition (Taylor1 jet transport)