Skip to content

Commit

Permalink
Parameterize the generate_ephemeris_2body test and update tolerance f…
Browse files Browse the repository at this point in the history 
…or new data
  • Loading branch information
@moeyensj
moeyensj committed Oct 4, 2023
1 parent 7792b6d commit 66087d1
Show file tree
Hide file tree
Showing 2 changed files with 121 additions and 107 deletions.
11 changes: 2 additions & 9 deletions adam_core/dynamics/tests/conftest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,16 +19,9 @@ def orbital_elements():

@pytest.fixture
def propagated_orbits():
propagated_orbits_file = files("adam_core.utils.helpers.data").joinpath(
"propagated_orbits.csv"
return Orbits.from_parquet(
files("adam_core.utils.helpers.data").joinpath("propagated_orbits.parquet")
)
df = pd.read_csv(
propagated_orbits_file,
index_col=False,
float_precision="round_trip",
dtype={"orbit_id": str},
)
return Orbits.from_dataframe(df, frame="ecliptic")


@pytest.fixture
Expand Down
217 changes: 119 additions & 98 deletions adam_core/dynamics/tests/test_ephemeris.py
View file
Original file line number Diff line number Diff line change
@@ -1,108 +1,129 @@
import numpy as np
import pytest
import quivr as qv
from astropy import units as u
from astropy.time import Time

from ...observers import Observers
from ..ephemeris import generate_ephemeris_2body


def test_generate_ephemeris_2body(propagated_orbits, ephemeris):
OBJECT_IDS = [
"594913 'Aylo'chaxnim (2020 AV2)",
"163693 Atira (2003 CP20)",
"(2010 TK7)",
"3753 Cruithne (1986 TO)",
"54509 YORP (2000 PH5)",
"2063 Bacchus (1977 HB)",
"1221 Amor (1932 EA1)",
"433 Eros (A898 PA)",
"3908 Nyx (1980 PA)",
"434 Hungaria (A898 RB)",
"1876 Napolitania (1970 BA)",
"2001 Einstein (1973 EB)",
"2 Pallas (A802 FA)",
"6 Hebe (A847 NA)",
"6522 Aci (1991 NQ)",
"10297 Lynnejones (1988 RJ13)",
"17032 Edlu (1999 FM9)",
"202930 Ivezic (1998 SG172)",
"911 Agamemnon (A919 FB)",
"1143 Odysseus (1930 BH)",
"1172 Aneas (1930 UA)",
"3317 Paris (1984 KF)",
"5145 Pholus (1992 AD)",
"5335 Damocles (1991 DA)",
"15760 Albion (1992 QB1)",
"15788 (1993 SB)",
"15789 (1993 SC)",
"1I/'Oumuamua (A/2017 U1)",
]

TOLERANCES = {
# range (m), angular difference (mas), light_time difference (microseconds)
# Default range is about consistent with our ability to get
# observatory state vector positions compared to Horizons
"default": (20, 0.1, 10),
# 594913 'Aylo'chaxnim (2020 AV2) : Vatira
# Range difference for first 30 observations between
# 32.8 m and 80.4 m, for last 60 observations range is
# between 0.13 m and 26.3 m
"594913 'Aylo'chaxnim (2020 AV2)": (100, 1, 10),
# 3753 Cruithne (1986 TO) : NEO (ATEN)
# Range anywhere between 0.01 m and 40.5 m throughout
# the 90 observations
"3753 Cruithne (1986 TO)": (50, 1, 10),
# 6522 Aci (1991 NQ) : Inner Main Belt asteroid
# 89/90 RA, Decs have a mean offset of 0.007 mas
# 1/90 RA, Decs has an offset of 1.24 mas
# 88/90 ranges have a mean offset of 10.6 m
# 2/90 ranges are 89.5, 174.3 m
"6522 Aci (1991 NQ)": (200, 2, 10),
# 1143 Odysseus (1930 BH) : Jupiter Trojan
# 89/90 light times have a mean offst of 0.14 microseconds
# 1/90 light times has an offset of 38.4 microseconds
# 89/90 RA, Decs have a mean offset of 0.0056 mas
# 1/90 RA, Decs has an offset of 8.9 mas
# 89/90 ranges have a mean offset of 6.9 m
# 1/90 ranges has na offset of 11575 m? ...
"1143 Odysseus (1930 BH)": (11600, 10, 50),
}


@pytest.mark.parametrize("object_id", OBJECT_IDS)
def test_generate_ephemeris_2body(object_id, propagated_orbits, ephemeris):
# For our catalog of test orbits, generate ephemeris using Horizons generated state vectors
# and compare the results to the Horizons generated ephemeris
orbit_ids = propagated_orbits.orbit_id.unique().to_numpy(zero_copy_only=False)
for orbit_id in orbit_ids:

propagated_orbit = propagated_orbits.select("orbit_id", orbit_id)
ephemeris_orbit = ephemeris[ephemeris["orbit_id"].isin([orbit_id])]

observers_list = []
for observatory_code in ephemeris_orbit["observatory_code"].unique():
observatory_mask = ephemeris_orbit["observatory_code"] == observatory_code
observer_i = Observers.from_code(
observatory_code,
Time(
ephemeris_orbit[observatory_mask]["mjd_utc"].values,
format="mjd",
scale="utc",
),
)
observers_list.append(observer_i)

observers = qv.concatenate(observers_list)

ephemeris_orbit_2body = generate_ephemeris_2body(propagated_orbit, observers)

# Extract only the ephemeris table
ephemeris_orbit_2body = ephemeris_orbit_2body.left_table

# Calculate the offset in light travel time
light_time_diff = np.abs(
ephemeris_orbit_2body.light_time.to_numpy()
- (ephemeris_orbit["lighttime"].values * 60 / 86400)
) * (u.d).to(u.microsecond)

# Assert that the light time is less than 10 microseconds
np.testing.assert_array_less(light_time_diff, 10)

# Calculate the difference in RA and Dec
angular_diff = np.abs(
ephemeris_orbit_2body.coordinates.values[:, 1:3]
- ephemeris_orbit[["RA", "DEC"]].values
) * (u.deg).to(u.milliarcsecond)

# Topocentric difference
range_diff = np.abs(
ephemeris_orbit_2body.coordinates.values[:, 0]
- ephemeris_orbit["delta"].values
) * (u.au).to(u.m)

# Orbit 00014 is 6522 Aci (1991 NQ) an inner main-belt asteroid
# with a diameter of about 5 km. What's interesting here is that of the 90
# ephemeris positions only 2 show a range difference of more than 15 m
# (many are less than 10 m). Their values are 92.17 and 177.24 m.
# These are the same observations that have offsets in RA and Dec greater
# than 1 milliarcsecond. The maximum offset for those observations not
# corresponding to the ones with a discrepant range offset is
# only 0.013 milliarcseconds!
if orbit_id == "00014":
# Assert that the position is less than 2 milliarcsecond
np.testing.assert_array_less(angular_diff, 2)

# Assert that the range is less than 200 m
np.testing.assert_array_less(range_diff, 200)

# Orbit 00019 is 1143 Odysseus (1930 BH), a Jupiter Trojan with a
# diameter of about 115 km. What's interesting here is that of the 90
# ephemeris positions only 2 show a range difference of more than 10 m.
# Their values are 231.79 and 468.26 m. These are the same observations
# that have offsets in RA and Dec greater than 0.5 milliarcsecond.
# However, the maximum offset for those observations not corresponding to the
# ones with a discrepant range offset is only 0.004 milliarcseconds!
elif orbit_id == "00019":
# Assert that the position is less than 1 milliarcsecond
np.testing.assert_array_less(angular_diff, 1)

# Assert that the range is less than 500 m
np.testing.assert_array_less(range_diff, 500)

# Orbit 00000 is 594913 'Aylo'chaxnim (2020 AV2) (an orbit completely
# interior to Venus' orbit). It would not surprise me that there are
# some dynamics that we are not modeling fully for this orbit in terms
# of light travel time and/or the ephemeris generation.
elif orbit_id == "00000":
# Assert that the position is less than 1 milliarcsecond
np.testing.assert_array_less(angular_diff, 1)

# Assert that the range is less than 200 m
np.testing.assert_array_less(range_diff, 200)

# These tolerance apply to the rest of the orbits (25/28) and
# show excellent results. RA, DEC to within 0.1 milliarcsecond
# and range to within 20 m.
else:
# Assert that the position is less than 0.1 milliarcsecond
np.testing.assert_array_less(angular_diff, 0.1)

# Assert that the range is less than 20 m
np.testing.assert_array_less(range_diff, 20)
propagated_orbit = propagated_orbits.select("object_id", object_id)
ephemeris_orbit = ephemeris[ephemeris["targetname"].isin([object_id])]

observers_list = []
for observatory_code in ephemeris_orbit["observatory_code"].unique():
observatory_mask = ephemeris_orbit["observatory_code"] == observatory_code
observer_i = Observers.from_code(
observatory_code,
Time(
ephemeris_orbit[observatory_mask]["mjd_utc"].values,
format="mjd",
scale="utc",
),
)
observers_list.append(observer_i)

observers = qv.concatenate(observers_list)

ephemeris_orbit_2body = generate_ephemeris_2body(propagated_orbit, observers)

# Extract only the ephemeris table
ephemeris_orbit_2body = ephemeris_orbit_2body.left_table

# Get the tolerances for this orbit
if object_id in TOLERANCES:
range_tolerance, angular_tolerance, light_time_tolerance = TOLERANCES[object_id]
else:
range_tolerance, angular_tolerance, light_time_tolerance = TOLERANCES["default"]

# Calculate the offset in light travel time
light_time_diff = np.abs(
ephemeris_orbit_2body.light_time.to_numpy()
- (ephemeris_orbit["lighttime"].values * 60 / 86400)
) * (u.d).to(u.microsecond)

# Assert that the light time is less than the tolerance
np.testing.assert_array_less(light_time_diff, light_time_tolerance)

# Calculate the difference in RA and Dec
angular_diff = np.abs(
ephemeris_orbit_2body.coordinates.values[:, 1:3]
- ephemeris_orbit[["RA", "DEC"]].values
) * (u.deg).to(u.milliarcsecond)

# Topocentric difference
range_diff = np.abs(
ephemeris_orbit_2body.coordinates.values[:, 0] - ephemeris_orbit["delta"].values
) * (u.au).to(u.m)

# Assert that the position is less than the tolerance
np.testing.assert_array_less(angular_diff, angular_tolerance)

# Assert that the range is less than the tolerance
np.testing.assert_array_less(range_diff, range_tolerance)

0 comments on commit 66087d1

Please sign in to comment.