Initial Orbit Determination

.github/workflows/CI.yml

@@ -4,10 +4,6 @@ on:
     paths-ignore:
       - 'LICENSE.md'
       - 'README.md'
-  pull_request:
-    branches:
-      - main
-    tags: '*'
 concurrency:
   # Skip intermediate builds: always.
   # Cancel intermediate builds: only if it is a pull request build.
@@ -27,7 +23,6 @@ jobs:
           - 'nightly'
         os:
           - ubuntu-latest
-          - macos-latest
           - windows-latest
         arch:
           - x64
@@ -43,9 +38,10 @@ jobs:
         env:
           JULIA_NUM_THREADS: 2
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
-          files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false
       - uses: coverallsapp/github-action@v2
         with:
           path-to-lcov: lcov.info

.gitignore

@@ -8,3 +8,5 @@ jldeph/
 notebooks/figs
 EOP_IAU*.TXT
 Manifest.toml
+*.DS_Store
+.vscode

Project.toml

@@ -1,11 +1,13 @@
 name = "NEOs"
 uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df"
 authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"]
-version = "0.7.5"
+version = "0.8.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
+Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -14,55 +16,44 @@ HORIZONS = "5a3ac768-beb4-554a-9c98-3342fe3377f5"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 Healpix = "9f4e344d-96bc-545a-84a3-ae6b9e1b672b"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 PlanetaryEphemeris = "d83715d0-7e5f-11e9-1a59-4137b20d8363"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
-Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062"
 SatelliteToolboxTransformations = "6b019ec1-7a1e-4f04-96c7-a9db1ca5514d"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[weakdeps]
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-
-[extensions]
-DataFramesExt = ["DataFrames", "Tables"]
-
 [compat]
 Artifacts = "1"
 AutoHashEquals = "0.2"
+Clustering = "0.15"
 DataFrames = "1.5"
 DelimitedFiles = "1"
 Downloads = "1"
-HORIZONS = "0.3"
+HORIZONS = "0.4"
 HTTP = "1.9.5"
 Healpix = "4"
-IntervalArithmetic = "0.20, 0.22"
-IntervalRootFinding = "0.5.11"
 JLD2 = "0.4"
+JSON = "0.21"
 LazyArtifacts = "1"
-PlanetaryEphemeris = "0.7"
+PlanetaryEphemeris = "0.8"
 Quadmath = "0.5"
-Requires = "0.5.2, 1"
 SPICE = "0.2"
 SatelliteToolboxTransformations = "0.1"
 Scratch = "1.2"
 StatsBase = "0.33, 0.34"
 Tables = "1.10"
-TaylorIntegration = "0.14"
-TaylorSeries = "0.16"
-julia = "1.6"
-
-[extras]
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+TaylorIntegration = "0.15"
+TaylorSeries = "0.17"
+julia = "1.9"

pha/Project.toml

@@ -14,4 +14,4 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 
 [compat]
-NEOs = "0.7"
+NEOs = "0.8"

pha/apophis.jl

@@ -124,41 +124,37 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     jd0 = datetime2julian(jd0_datetime)
 
     print_header("Integrator warmup", 2)
-    _ = NEOs.propagate(dynamics, 1, jd0, nyears_fwd, q0, Val(true);
-                       order, abstol, parse_eqs)
-    _ = NEOs.propagate_root(dynamics, 1, jd0, nyears_fwd, q0, Val(true);
-                       order, abstol, parse_eqs)
+    params = NEOParameters(;maxsteps=1, order, abstol, parse_eqs)
+    _ = NEOs.propagate(dynamics, jd0, nyears_fwd, q0, params)
+    _ = NEOs.propagate_root(dynamics, jd0, nyears_fwd, q0, params)
     println()
 
     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))
 
-    sol_bwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_bwd, q0, Val(true);
-                             order, abstol, parse_eqs)
+    params = NEOParameters(;maxsteps, order, abstol, parse_eqs)
+    sol_bwd = NEOs.propagate(dynamics, jd0, nyears_bwd, q0, params)
     jldsave("Apophis_bwd.jld2"; sol_bwd)
     # sol_bwd = JLD2.load("Apophis_bwd.jld2", "sol_bwd")
 
     tmax = nyears_fwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
     println("• Final time of integration: ", julian2datetime(jd0 + tmax))
 
-    sol_fwd, tvS, xvS, gvS = NEOs.propagate_root(dynamics, maxsteps, jd0, nyears_fwd, q0, Val(true);
-                             order, abstol, parse_eqs)
+    sol_fwd, tvS, xvS, gvS = NEOs.propagate_root(dynamics, jd0, nyears_fwd, q0, params)
     jldsave("Apophis_fwd.jld2"; sol_fwd, tvS, xvS, gvS)
     # sol_fwd = JLD2.load("Apophis_fwd.jld2", "sol_bwd")
     println()
 
-    # load Solar System ephemeris
-    sseph::TaylorInterpolant{Float64,Float64,2} = loadpeeph(NEOs.sseph, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end])
-    eph_su::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, su)
-    eph_ea::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, ea)
+    # Load Sun ephemeris
+    eph_su = selecteph(NEOs.sseph, su, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end])
+    # Load Earth ephemeris
+    eph_ea = selecteph(NEOs.sseph, ea, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end])
 
-    # NEO
+    # Apophis
     # Change t, x, v units, resp., from days, au, au/day to sec, km, km/sec
-    xva_bwd(et) = auday2kmsec(sol_bwd(et/daysec)[1:6])
-    xva_fwd(et) = auday2kmsec(sol_fwd(et/daysec)[1:6])
     xva(et) = bwdfwdeph(et, sol_bwd, sol_fwd)
     # Earth
     # Change x, v units, resp., from au, au/day to km, km/sec
@@ -177,7 +173,8 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     deldop = vcat(deldop_2005_2013,deldop_2021)
 
     # Compute optical residuals
-    res_radec_all, w_radec_all = NEOs.residuals(radec; xvs, xve, xva)
+    _res_radec_all_ = NEOs.residuals(radec; xvs, xve, xva)
+    res_radec_all, w_radec_all = NEOs.unfold(_res_radec_all_)
     jldsave("Apophis_res_w_radec.jld2"; res_radec_all, w_radec_all)
     # JLD2.@load "Apophis_res_w_radec.jld2"
 
@@ -233,9 +230,9 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     res = vcat(res_radec, res_del, res_dop)
     w = vcat(w_radec, w_del, w_dop)
 
-    success, δx_OR8, Γ_OR8 = newtonls(res, w, zeros(get_numvars()), 10)
-    x_OR8 = sol_fwd(sol_fwd.t0)(δx_OR8)
-    σ_OR8 = sqrt.(diag(Γ_OR8)).*scalings
+    fit_OR8 = newtonls(res, w, zeros(get_numvars()), 10)
+    x_OR8 = sol_fwd(sol_fwd.t0)(fit_OR8.x)
+    σ_OR8 = sqrt.(diag(fit_OR8.Γ)).*scalings
 
     nradec = length(res_radec)
     res_ra = view(res_radec, 1:nradec÷2)
@@ -248,19 +245,19 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     print_header("Orbital fit (8-DOF) and post-fit statistics", 2)
 
     # orbital fit
-    println("Success flag                                                      : ", success, "\n")
+    println("Success flag                                                      : ", fit_OR8.success, "\n")
     println("Nominal solution             [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", x_OR8, "\n")
     println("1-sigma formal uncertainties [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", σ_OR8, "\n")
 
     # post-fit statistics
-    println("Normalized RMS (optical-only)               [adimensional] : ", nrms(res_radec(δx_OR8),w_radec))
-    println("Normalized RMS (radar-only)                 [adimensional] : ", nrms(vcat(res_del,res_dop)(δx_OR8),vcat(w_del,w_dop)))
-    println("Normalized RMS (combined optical and radar) [adimensional] : ", nrms(res(δx_OR8),w), "\n")
-    println("Mean weighted right-ascension residual      [arcseconds]   : ", mean(res_ra(δx_OR8), weights(w_ra)))
-    println("Mean weighted declination residual          [arcseconds]   : ", mean(res_dec(δx_OR8), weights(w_dec)))
-    println("Mean weighted time-delay residual           [micro-seconds]: ", mean(res_del(δx_OR8), weights(w_del)))
-    println("Mean weighted Doppler-shift residual        [Hz]           : ", mean(res_dop(δx_OR8), weights(w_dop)), "\n")
-    println("Chi-squared statistic (χ²):                 [adimensional] : ", chi2(res(δx_OR8),w))
+    println("Normalized RMS (optical-only)               [adimensional] : ", nrms(res_radec(fit_OR8.x),w_radec))
+    println("Normalized RMS (radar-only)                 [adimensional] : ", nrms(vcat(res_del,res_dop)(fit_OR8.x),vcat(w_del,w_dop)))
+    println("Normalized RMS (combined optical and radar) [adimensional] : ", nrms(res(fit_OR8.x),w), "\n")
+    println("Mean weighted right-ascension residual      [arcseconds]   : ", mean(res_ra(fit_OR8.x), weights(w_ra)))
+    println("Mean weighted declination residual          [arcseconds]   : ", mean(res_dec(fit_OR8.x), weights(w_dec)))
+    println("Mean weighted time-delay residual           [micro-seconds]: ", mean(res_del(fit_OR8.x), weights(w_del)))
+    println("Mean weighted Doppler-shift residual        [Hz]           : ", mean(res_dop(fit_OR8.x), weights(w_dop)), "\n")
+    println("Chi-squared statistic (χ²):                 [adimensional] : ", chi2(res(fit_OR8.x),w))
 
     return sol_bwd, sol_fwd, res_radec, res_del, res_dop, w_radec, w_del, w_dop
 end
@@ -316,3 +313,36 @@ function main()
 end
 
 main()
+
+#=
+-0.18034827489412805
+0.9406910783153754
+0.3457360118643932
+-0.016265940057745887
+4.3915296805381036e-5
+-0.0003952032399008921
+-2.883842658925719e-14
+-1.6977850502837174e-11
+
+
+6.991681792970175e-9
+3.3210758500485412e-9
+8.763895947794195e-9
+5.485729074188412e-11
+9.846074910774499e-11
+1.3782767177665416e-10
+2.525644870489942e-16
+3.2193289354129716e-12
+
+Normalized RMS (optical-only)               [adimensional] : 0.4365429773707507
+Normalized RMS (radar-only)                 [adimensional] : 0.5421268387032326
+Normalized RMS (combined optical and radar) [adimensional] : 0.43685679642513925
+
+Mean weighted right-ascension residual      [arcseconds]   : -0.002787272519210798
+Mean weighted declination residual          [arcseconds]   : -0.003376099885054764
+Mean weighted time-delay residual           [micro-seconds]: -0.0007329351962194035
+Mean weighted Doppler-shift residual        [Hz]           : -0.03025601266756564
+
+Chi-squared statistic (χ²):                 [adimensional] : 3597.4067719864497
+
+=#