(20505) Makadia orbit determination test

Let’s start by importing the necessary libraries

from grss import fit
import numpy as np
np.set_printoptions(precision=40, linewidth=np.inf)

We’ll then retrieve the cometary state of the asteroid (from JPL SBDB) plus any nongravitational accelerations acting on it.

body_id = '20505'
init_sol, init_cov, nongrav_info = fit.get_sbdb_info(body_id)
de_kernel = 440

Next, we’ll retrieve the observations from different sources (MPC, JPL, Gaia Data Releases) and prepare them for the orbit determination process.

add_gaia_obs = True
optical_obs_file = None
t_min_tdb = None
t_max_tdb = None
debias_lowres = True
deweight = True
eliminate = False
num_obs_per_night = 4
verbose = True
obs_df = fit.get_optical_obs(body_id, optical_obs_file, t_min_tdb, t_max_tdb, debias_lowres, deweight, eliminate, num_obs_per_night, verbose)
obs_df = fit.add_radar_obs(obs_df, t_min_tdb, t_max_tdb, verbose)
if add_gaia_obs:
    gaia_dr = 'gaiafpr'
    obs_df = fit.add_gaia_obs(obs_df, t_min_tdb, t_max_tdb, gaia_dr, verbose)

Read in 4291 observations from the MPC.
        Filtered to 4291 observations that satisfy the time range and accepted observatory constraints.
Applying Eggl et al. (2020) debiasing scheme to the observations.
        Unknown star catalog: UNK
        No debiasing needed for 2762 observations.
        Debiased 1443 observations.
        No bias information for 86 observations.
Applying Vereš et al. (2017) weighting scheme to the observations.
        Using 3723 CCD observations with station-specific weight rules.
Applying sqrt(N/4) deweighting scheme.
        Deweighted 1037 observations.
Read in 250 Gaia observations from gaiafpr
        Filtered to 250 observations that satisfy the time range constraints.

All we need to do now is initialize the OD simulation and run the filter.

n_iter_max = 10
fit_sim = fit.FitSimulation(init_sol, obs_df, init_cov, n_iter_max=n_iter_max, de_kernel=de_kernel, nongrav_info=nongrav_info)

fit_sim.filter_lsq()

Iteration               Unweighted RMS          Weighted RMS            Chi-squared             Reduced Chi-squared
1                        0.456                   0.451                   1845.416                        0.203
2                        0.456                   0.451                   1845.388                        0.203
Converged without rejecting outliers. Starting outlier rejection now...
3                        0.367                   0.438                   1739.031                        0.192
4                        0.349                   0.413                   1546.658                        0.171
5                        0.349                   0.413                   1546.636                        0.171
Converged after rejecting outliers. Rejected 16 out of 4541 optical observations.

Let’s print some summary statistics and plot some results.

fit_sim.print_summary()

Summary of the orbit fit calculations after postfit pass:
==============================================================
RMS unweighted: 0.3490486724130132
RMS weighted: 0.4126703642911768
chi-squared: 1546.6358061021785
reduced chi-squared: 0.1710123624615412
square root of reduced chi-squared: 0.4135364100796219
--------------------------------------------------------------
Solution Time: MJD 58205.000 TDB = 2018-03-28 00:00:00.000 TDB
Solution Observation Arc: 27605.39 days (75.58 years)
--------------------------------------------------------------
Fitted Variable         Initial Value                   Uncertainty                     Fitted Value                    Uncertainty                     Change                          Change (sigma)
e                       7.56295317434e-02               6.49222933173e-10               7.56295317395e-02               6.23546253029e-10               -3.86995990809e-12              -0.006
q                       2.92366099974e+00               3.22439017915e-09               2.92366099967e+00               3.10337264134e-09               -7.32263139014e-11              -0.024
tp                      5.73464369975e+04               4.48124464705e-06               5.73464369974e+04               4.34662273989e-06               -1.23691279441e-07              -0.028
om                      3.25422007829e+02               6.01043159569e-07               3.25422007820e+02               5.79839525373e-07               -8.68203642312e-09              -0.015
w                       3.29371665601e+02               9.89492916636e-07               3.29371665584e+02               9.60511113175e-07               -1.72169052348e-08              -0.018
i                       5.29547237699e+00               6.25060698811e-08               5.29547237823e+00               6.06056503833e-08               +1.24181553929e-09              +0.020

fit_sim.plot_summary(auto_close=True)

fit_sim.iters[-1].plot_iteration_summary(title='Postfit Residuals', auto_close=True)

mean_0 = np.array(list(init_sol.values())[1:])
cov_0 = init_cov
mean_f = np.array(list(fit_sim.x_nom.values()))
cov_f = fit_sim.covariance

maha_dist_f, maha_dist_0, bhattacharya, bhatt_coeff = fit.get_similarity_stats(mean_0, cov_0, mean_f, cov_f)
print(f'Mahalonobis distance between JPL and GRSS solution: {maha_dist_f:0.2f}')
print(f'Mahalonobis distance between GRSS and JPL solution: {maha_dist_0:0.2f}')
print(f'Bhattacharya distance between JPL and GRSS solution: {bhattacharya:0.4f}')
print(f'Bhattacharya coefficient between JPL and GRSS solution: {bhatt_coeff:0.4f}')

Mahalonobis distance between JPL and GRSS solution: 0.05
Mahalonobis distance between GRSS and JPL solution: 0.06
Bhattacharya distance between JPL and GRSS solution: 0.0021
Bhattacharya coefficient between JPL and GRSS solution: 0.9980

Finally, we’ll make sure the GRSS solution is statistically consistent with the JPL SBDB solution

assert maha_dist_f < 5.0
assert maha_dist_0 < 5.0
assert bhattacharya < 0.10
assert bhatt_coeff > 0.90