Upgrade synthetic test [ADAM-57]

[moeyensj]

This PR adds two new scripts:

• get_orbits.py: query JPL Small-Body Database for orbits of sample objects and saves to a CSV.
• make_observations.py: read orbits and generate an observations file saved as HDF5 for the sample objects.

Now test_precovery() tests the following:

• ra, ra_sigma, dec, dec_sigma, mag, mag_sigma, mjd_utc, obs_id, exposure_id, observatory_code, filter are equal between input observations and indexed observations
• input ra, dec is equal to predicted ra and dec reported by precovery (this fails at the moment, debugging to understand)
• distances calculated by precovery are zero (this fails at the moment because of the above)

[review-notebook-app[bot]]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[moeyensj]

While testing the above, I stumbled into a time scale bug. The Orbit class was using the time scale with which it was initialized as the time scale for both propagation and ephemeris generation. So when you initialized an orbit with MJDs in TT then all propagations and ephemeris generation times were assumed to be TT. However, our observations are in UTC. This was solved by adding a time_scale kwarg to both orbit.propagate and orbit.compute_ephemeris.

The new unit test is still failing but I've narrowed it down to floating point error accumulation as a result of a different number of computations used to generate the test data compared to how precovery finds observations.

[moeyensj]

I've added a README to explain what the scripts do and some of the dataset assumptions. Let me know what you think, @akoumjian.

I've been able to narrow down the discrepancy between the predicted RA and Dec returned by precovery and the input RA and Dec generated by make_observations.py to the accumulation of floating point errors. Precovery will propagate, store, then propagate the resulting orbit, and so on. The number of propagations precovery will conduct will depend at least on the number of frames and the window_size parameter. On the other hand, generating the synthetic observations uses a single propagation from the reference epoch to the observation time. The differences between the two positions here are of order ~1e-12, so that's the tolerance to which I have set the comparison between the inputs and predicted outputs in the test. I tested this by adding a few additional propagations inside the make_observations.py and modifying the window size and at that point we start seeing changes in the difference between the two arrays, which tells me it comes down to the number of propagations.

Here is a summary of some of the "interesting" features I've encountered while debugging the above:

• PYOORB cannot invert a Keplerian orbit with a negative semi-major axis (hyperbolic orbits, e > 1.0). These orbits need to be expressed in Cometary elements for that inversion to work. I've kept our ISO sample object "commented out" to bypass this for the time being.
• PYOORB has inconsistent behavior between Cometary and Keplerian element inversion. When we generate the observations using the Cometary elements and not the default Keplerian elements, the unit tests fail.

The second case can be tested by running:

python make_observations.py --orbit_type=cometary
pytest test_precovery_index.py

This will cause the unit test to fail:

AssertionError: 
>           np.testing.assert_allclose(
                results[["pred_ra_deg", "pred_dec_deg"]].values,
                object_observations[["ra", "dec"]].values,
                atol=1e-12,
                rtol=1e-12,
            )
E           AssertionError: 
E           Not equal to tolerance rtol=1e-12, atol=1e-12
E           
E           Mismatched elements: 120 / 120 (100%)
E           Max absolute difference: 3.18053139e-09
E           Max relative difference: 3.97825811e-10
E            x: array([[143.643415,  17.771414],
E                  [143.675804,  17.755943],
E                  [143.708188,  17.740464],...
E            y: array([[143.643415,  17.771414],
E                  [143.675804,  17.755943],
E                  [143.708188,  17.740464],...

test_precovery_index.py:117: AssertionError

In other words, with Cometary elements, the best we can do is get agreement to within ~360 nanoarseconds. While with Keplerian we can get agreement to within ~3.6 nanoarcseconds. This is likely not a significant issue since observational errors in practice will be ~10s of milliarcseconds. The difference here may come down to differences in the algorithm used to invert the two representations to Cartesian.

Note: fundamentally Keplerian elements cannot be used to express parabolic orbits since the semi-major axis is undefined when e = 1.0. These orbits are always definable with cometary elements, as is the case with ellipses (e < 1.0), circles (e = 0.0), and hyperbolas (e > 1.0).