From 6e8c9c58909a2007299aa17d4bd314d169c39008 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Fri, 26 Jul 2024 12:29:46 -0400 Subject: [PATCH 01/11] Add generic Ephemeris implementation to Propagator with light time correction --- src/adam_core/_version.py | 1 + src/adam_core/orbits/query/horizons.py | 5 +- src/adam_core/propagator/propagator.py | 162 +++++++++++++++++++++++-- 3 files changed, 156 insertions(+), 12 deletions(-) create mode 100644 src/adam_core/_version.py diff --git a/src/adam_core/_version.py b/src/adam_core/_version.py new file mode 100644 index 00000000..23f4fc81 --- /dev/null +++ b/src/adam_core/_version.py @@ -0,0 +1 @@ +__version__ = '0.2.1.dev1+g7f7fcd5.d20240725' \ No newline at end of file diff --git a/src/adam_core/orbits/query/horizons.py b/src/adam_core/orbits/query/horizons.py index ec9d02ce..876303a1 100644 --- a/src/adam_core/orbits/query/horizons.py +++ b/src/adam_core/orbits/query/horizons.py @@ -157,10 +157,11 @@ def _get_horizons_ephemeris( as seen from the observer location at the given times. """ dfs = [] + mjd_utc = times.rescale("utc").mjd().to_numpy(zero_copy_only=False) for i, obj_id in enumerate(object_ids): obj = Horizons( id=obj_id, - epochs=times.rescale("utc").mjd().to_numpy(zero_copy_only=False), + epochs=mjd_utc, location=location, id_type=id_type, ) @@ -171,7 +172,7 @@ def _get_horizons_ephemeris( cache=False, ).to_pandas() ephemeris.insert(0, "orbit_id", f"{i:05d}") - ephemeris.insert(2, "mjd_utc", times.utc.mjd) + ephemeris.insert(2, "mjd_utc", mjd_utc) ephemeris.insert(3, "observatory_code", location) dfs.append(ephemeris) diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 0f23da35..9e9cda26 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -6,12 +6,16 @@ import numpy.typing as npt import quivr as qv -from adam_core.ray_cluster import initialize_use_ray - +from ..coordinates.cartesian import CartesianCoordinates +from ..coordinates.origin import Origin, OriginCodes +from ..coordinates.spherical import SphericalCoordinates +from ..coordinates.transform import transform_coordinates +from ..dynamics.aberrations import C from ..observers.observers import Observers from ..orbits.ephemeris import Ephemeris from ..orbits.orbits import Orbits from ..orbits.variants import VariantEphemeris, VariantOrbits +from ..ray_cluster import initialize_use_ray from ..time import Timestamp from .utils import _iterate_chunks @@ -89,17 +93,155 @@ class EphemerisMixin: Subclasses should implement the _generate_ephemeris method. """ - @abstractmethod + def _add_light_time( + self, + orbits, + observers, + lt_tol: float = 1e-10, + ): + """ + Accounts for light time correction in the ephemeris generation. + + This differs from the _add_light_time method in the aberrations module + in that is uses the underlying propagator to propagate the orbits as opposed + to the _propagate_2body function. + """ + orbits_aberrated = orbits.empty() + lts = np.zeros(len(orbits)) + for i, (orbit, observer) in enumerate(zip(orbits, observers)): + lt0 = 0 + dlt = 1 + iterations = 0 + while dlt > lt_tol: + iterations += 1 + if iterations > 3: + break + observer_position = observer.coordinates.values[0,:3] + orbit_i = orbit + t0 = orbit_i.coordinates.time.mjd()[0].as_py() + + rho = np.linalg.norm(orbit_i.coordinates.values[0,:3] - observer_position) + + lt = rho / C + + dlt = np.abs(lt - lt0) + + t1 = t0 - lt + t1 = Timestamp.from_mjd([t1], scale="tdb") + orbit_propagated = self.propagate_orbits(orbit, t1) + + orbit_i = orbit_propagated + t0 = t1 + lt0 = lt + dlt = dlt + orbits_aberrated = qv.concatenate([orbits_aberrated, orbit]) + lts[i] = lt + + return orbits_aberrated, lts + def _generate_ephemeris( - self, orbits: EphemerisType, observers: ObserverType + self, orbits: OrbitType, observers: ObserverType, lt_tol: float = 1e-10 ) -> EphemerisType: + """ - Generate ephemerides for the given orbits as observed by - the observers. - - THIS FUNCTION SHOULD BE DEFINED BY THE USER. + A generic ephemeris implementation, which can be used or overridden by subclasses. """ - pass + + if isinstance(orbits, Orbits): + ephemeris_total = Ephemeris.empty() + elif isinstance(orbits, VariantOrbits): + ephemeris_total = VariantEphemeris.empty() + + for orbit in orbits: + propagated_orbits = self.propagate_orbits( + orbit, observers.coordinates.time + ) + + # Transform both the orbits and observers to the barycenter if they are not already. + propagated_orbits_barycentric = propagated_orbits.set_column( + "coordinates", + transform_coordinates( + propagated_orbits.coordinates, + CartesianCoordinates, + frame_out="ecliptic", + origin_out=OriginCodes.SOLAR_SYSTEM_BARYCENTER, + ), + ) + observers_barycentric = observers.set_column( + "coordinates", + transform_coordinates( + observers.coordinates, + CartesianCoordinates, + frame_out="ecliptic", + origin_out=OriginCodes.SOLAR_SYSTEM_BARYCENTER, + ), + ) + + num_orbits = len(propagated_orbits_barycentric.orbit_id.unique()) + + observer_codes = np.tile( + observers.code.to_numpy(zero_copy_only=False), num_orbits + ) + + propagated_orbits_aberrated, light_time = self._add_light_time( + propagated_orbits_barycentric, + observers_barycentric, + lt_tol=lt_tol, + ) + + topocentric_coordinates = CartesianCoordinates.from_kwargs( + x=propagated_orbits_aberrated.coordinates.values[:, 0] + - observers_barycentric.coordinates.values[:, 0], + y=propagated_orbits_aberrated.coordinates.values[:, 1] + - observers_barycentric.coordinates.values[:, 1], + z=propagated_orbits_aberrated.coordinates.values[:, 2] + - observers_barycentric.coordinates.values[:, 2], + vx=propagated_orbits_aberrated.coordinates.values[:, 3] + - observers_barycentric.coordinates.values[:, 3], + vy=propagated_orbits_aberrated.coordinates.values[:, 4] + - observers_barycentric.coordinates.values[:, 4], + vz=propagated_orbits_aberrated.coordinates.values[:, 5] + - observers_barycentric.coordinates.values[:, 5], + covariance=None, + time=propagated_orbits_aberrated.coordinates.time, + origin=Origin.from_kwargs(code=observer_codes), + frame="ecliptic", + ) + + spherical_coordinates = SphericalCoordinates.from_cartesian( + topocentric_coordinates + ) + light_time = np.array(light_time) + + spherical_coordinates = transform_coordinates( + spherical_coordinates, SphericalCoordinates, frame_out="equatorial" + ) + + if isinstance(orbits, Orbits): + + ephemeris = Ephemeris.from_kwargs( + orbit_id=propagated_orbits_barycentric.orbit_id, + object_id=propagated_orbits_barycentric.object_id, + coordinates=spherical_coordinates, + light_time=light_time, + aberrated_coordinates=propagated_orbits_aberrated.coordinates, + ) + + elif isinstance(orbits, VariantOrbits): + weights = orbits.weights + weights_cov = orbits.weights_cov + + ephemeris = VariantEphemeris.from_kwargs( + orbit_id=propagated_orbits_barycentric.orbit_id, + object_id=propagated_orbits_barycentric.object_id, + coordinates=spherical_coordinates, + weights=weights, + weights_cov=weights_cov, + ) + + ephemeris_total = qv.concatenate([ephemeris_total, ephemeris]) + + return ephemeris_total def generate_ephemeris( self, @@ -261,7 +403,7 @@ def generate_ephemeris( ) -class Propagator(ABC): +class Propagator(ABC, EphemerisMixin): """ Abstract class for propagating orbits and related functions. From 0e9e8c2a54fe923670497e6deea60d7fbe93b6cf Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Fri, 26 Jul 2024 13:30:16 -0400 Subject: [PATCH 02/11] typing --- src/adam_core/coordinates/residuals.py | 2 +- src/adam_core/coordinates/transform.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/adam_core/coordinates/residuals.py b/src/adam_core/coordinates/residuals.py index 71490b83..e8c0cbec 100644 --- a/src/adam_core/coordinates/residuals.py +++ b/src/adam_core/coordinates/residuals.py @@ -194,7 +194,7 @@ def calculate( raise TypeError( f"Predicted coordinates must be one of {SUPPORTED_COORDINATES}, not {type(predicted)}." ) - if type(observed) != type(predicted): + if type(observed) is not type(predicted): raise TypeError( "Observed and predicted coordinates must be the same type, " f"not {type(observed)} and {type(predicted)}." diff --git a/src/adam_core/coordinates/transform.py b/src/adam_core/coordinates/transform.py index d4265f47..f531d5c9 100644 --- a/src/adam_core/coordinates/transform.py +++ b/src/adam_core/coordinates/transform.py @@ -1368,7 +1368,7 @@ def transform_coordinates( # `~adam_core.coordinates.origin.OriginCodes` so we can compare them directly. # If its not an OriginCodes enum then origin_out will be an array of strings which # also can be checked for equality. - if type(coords) == representation_out_: + if type(coords) is representation_out_: if coord_frame == frame_out and np.all(coord_origin == origin_out): return coords From a8e93ec959c36cb7aad05ab33fbe8b8db463a293 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Fri, 26 Jul 2024 13:30:58 -0400 Subject: [PATCH 03/11] try now --- src/adam_core/propagator/propagator.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 9e9cda26..bb402a87 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -94,7 +94,7 @@ class EphemerisMixin: """ def _add_light_time( - self, + self, orbits, observers, lt_tol: float = 1e-10, @@ -116,11 +116,13 @@ def _add_light_time( iterations += 1 if iterations > 3: break - observer_position = observer.coordinates.values[0,:3] + observer_position = observer.coordinates.values[0, :3] orbit_i = orbit t0 = orbit_i.coordinates.time.mjd()[0].as_py() - rho = np.linalg.norm(orbit_i.coordinates.values[0,:3] - observer_position) + rho = np.linalg.norm( + orbit_i.coordinates.values[0, :3] - observer_position + ) lt = rho / C @@ -142,7 +144,6 @@ def _add_light_time( def _generate_ephemeris( self, orbits: OrbitType, observers: ObserverType, lt_tol: float = 1e-10 ) -> EphemerisType: - """ A generic ephemeris implementation, which can be used or overridden by subclasses. """ @@ -153,9 +154,7 @@ def _generate_ephemeris( ephemeris_total = VariantEphemeris.empty() for orbit in orbits: - propagated_orbits = self.propagate_orbits( - orbit, observers.coordinates.time - ) + propagated_orbits = self.propagate_orbits(orbit, observers.coordinates.time) # Transform both the orbits and observers to the barycenter if they are not already. propagated_orbits_barycentric = propagated_orbits.set_column( From 4d3da0c1be9b28ff2e03f8dcad7f97de0971534d Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Fri, 26 Jul 2024 14:00:00 -0400 Subject: [PATCH 04/11] _version.py formatting --- src/adam_core/_version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adam_core/_version.py b/src/adam_core/_version.py index 23f4fc81..73efbb21 100644 --- a/src/adam_core/_version.py +++ b/src/adam_core/_version.py @@ -1 +1 @@ -__version__ = '0.2.1.dev1+g7f7fcd5.d20240725' \ No newline at end of file +__version__ = "0.2.1.dev1+g7f7fcd5.d20240725" From 5d8d3bdcd7eba1f47dd6730d5652b0dd22979832 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 1 Aug 2024 11:45:09 -0400 Subject: [PATCH 05/11] Update horizons to use adam-core clases and update light time calculations --- src/adam_core/_version.py | 2 +- src/adam_core/coordinates/origin.py | 12 ++++ src/adam_core/orbits/query/horizons.py | 40 ++++++++++--- src/adam_core/propagator/propagator.py | 77 ++++++++++++++------------ src/adam_core/time/time.py | 3 +- 5 files changed, 86 insertions(+), 48 deletions(-) diff --git a/src/adam_core/_version.py b/src/adam_core/_version.py index 73efbb21..3261ed66 100644 --- a/src/adam_core/_version.py +++ b/src/adam_core/_version.py @@ -1 +1 @@ -__version__ = "0.2.1.dev1+g7f7fcd5.d20240725" +__version__ = '0.2.1.dev5+g4d3da0c.d20240731' \ No newline at end of file diff --git a/src/adam_core/coordinates/origin.py b/src/adam_core/coordinates/origin.py index 38775339..b75f0702 100644 --- a/src/adam_core/coordinates/origin.py +++ b/src/adam_core/coordinates/origin.py @@ -86,6 +86,18 @@ def SOLAR_SYSTEM_BARYCENTER(cls) -> float: class Origin(qv.Table): code = qv.LargeStringColumn() + def as_OriginCodes(self) -> OriginCodes: + """ + Convert the origin codes to an `~adam_core.coordinates.origin.OriginCodes` object. + + Returns + ------- + OriginCodes + Origin codes as an `~adam_core.coordinates.origin.OriginCodes` object. + """ + assert len(self.code.unique()) == 1, "Only one origin code can be converted at a time." + return OriginCodes[self.code.unique()[0].as_py()] + def __eq__(self, other: object) -> np.ndarray: if isinstance(other, (str, np.ndarray)): codes = self.code.to_numpy(zero_copy_only=False) diff --git a/src/adam_core/orbits/query/horizons.py b/src/adam_core/orbits/query/horizons.py index 876303a1..9bbe4866 100644 --- a/src/adam_core/orbits/query/horizons.py +++ b/src/adam_core/orbits/query/horizons.py @@ -2,14 +2,17 @@ import numpy.typing as npt import pandas as pd +import pyarrow as pa from astroquery.jplhorizons import Horizons from ...coordinates.cartesian import CartesianCoordinates from ...coordinates.cometary import CometaryCoordinates from ...coordinates.keplerian import KeplerianCoordinates from ...coordinates.origin import Origin +from ...coordinates.spherical import SphericalCoordinates from ...observers import Observers from ...time import Timestamp +from ..ephemeris import Ephemeris from ..orbits import Orbits @@ -53,7 +56,7 @@ def _get_horizons_vectors( for i, obj_id in enumerate(object_ids): obj = Horizons( id=obj_id, - epochs=times.rescale("tdb").mjd().to_numpy(zero_copy_only=False), + epochs=times.rescale("tdb").jd().to_numpy(zero_copy_only=False), location=location, id_type=id_type, ) @@ -157,11 +160,11 @@ def _get_horizons_ephemeris( as seen from the observer location at the given times. """ dfs = [] - mjd_utc = times.rescale("utc").mjd().to_numpy(zero_copy_only=False) + jd_utc = times.rescale("utc").jd().to_numpy(zero_copy_only=False) for i, obj_id in enumerate(object_ids): obj = Horizons( id=obj_id, - epochs=mjd_utc, + epochs=jd_utc, location=location, id_type=id_type, ) @@ -172,7 +175,7 @@ def _get_horizons_ephemeris( cache=False, ).to_pandas() ephemeris.insert(0, "orbit_id", f"{i:05d}") - ephemeris.insert(2, "mjd_utc", mjd_utc) + ephemeris.insert(2, "jd_utc", jd_utc) ephemeris.insert(3, "observatory_code", location) dfs.append(ephemeris) @@ -188,7 +191,7 @@ def _get_horizons_ephemeris( def query_horizons_ephemeris( object_ids: Union[List, npt.ArrayLike], observers: Observers -) -> pd.DataFrame: +) -> Ephemeris: """ Query JPL Horizons (through astroquery) for an object's predicted ephemeris as seen from a given location at the given times. @@ -209,19 +212,38 @@ def query_horizons_ephemeris( """ dfs = [] for observatory_code, observers_i in observers.iterate_codes(): - ephemeris = _get_horizons_ephemeris( + _ephemeris = _get_horizons_ephemeris( object_ids, observers_i.coordinates.time, observatory_code, ) - dfs.append(ephemeris) + dfs.append(_ephemeris) - ephemeris = pd.concat(dfs, ignore_index=True) - ephemeris.sort_values( + dfs = pd.concat(dfs, ignore_index=True) + dfs.sort_values( by=["orbit_id", "datetime_jd", "observatory_code"], inplace=True, ignore_index=True, ) + + # Horizons produces UTC but we use tdb everywhere + epochs = Timestamp.from_jd(pa.array(dfs["datetime_jd"]), scale="utc").rescale("tdb") + + ephemeris = Ephemeris.from_kwargs( + orbit_id=dfs["orbit_id"], + object_id=dfs["targetname"], + # Convert from minutes to days + light_time=dfs["lighttime"] / 1440, + alpha=dfs["alpha"], + coordinates=SphericalCoordinates.from_kwargs( + time=epochs, + lon=dfs["RA"], + lat=dfs["DEC"], + origin=Origin.from_kwargs(code=dfs["observatory_code"]), + frame="ecliptic", + ), + ) + return ephemeris diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index bb402a87..4aa63261 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -94,53 +94,56 @@ class EphemerisMixin: """ def _add_light_time( - self, - orbits, - observers, - lt_tol: float = 1e-10, - ): - """ - Accounts for light time correction in the ephemeris generation. - - This differs from the _add_light_time method in the aberrations module - in that is uses the underlying propagator to propagate the orbits as opposed - to the _propagate_2body function. - """ + self, + orbits, + observers, + lt_tol: float = 1e-12, + max_iter: int = 10, + ): orbits_aberrated = orbits.empty() lts = np.zeros(len(orbits)) for i, (orbit, observer) in enumerate(zip(orbits, observers)): - lt0 = 0 - dlt = 1 + # Set the running variables + lt_prev = 0 + dlt = float('inf') + orbit_i = orbit + lt = 0 + + # Extract the observer's position which remains + # constant for all iterations + observer_position = observer.coordinates.r + + # Calculate the orbit's current epoch (the epoch from which + # the light travel time will be calculated) + t0 = orbit_i.coordinates.time.rescale("tdb").mjd()[0].as_py() + + iterations = 0 - while dlt > lt_tol: + while dlt > lt_tol and iterations < max_iter: iterations += 1 - if iterations > 3: - break - observer_position = observer.coordinates.values[0, :3] - orbit_i = orbit - t0 = orbit_i.coordinates.time.mjd()[0].as_py() - - rho = np.linalg.norm( - orbit_i.coordinates.values[0, :3] - observer_position - ) + # Calculate the topocentric distance + rho = np.linalg.norm(orbit_i.coordinates.r - observer_position) + + # Calculate the light travel time lt = rho / C - dlt = np.abs(lt - lt0) + # Calculate the change in light travel time since the previous iteration + dlt = np.abs(lt - lt_prev) + + # Calculate the new epoch and propagate the initial orbit to that epoch + orbit_i = self.propagate_orbits(orbit, Timestamp.from_mjd([t0 - lt], scale="tdb")) - t1 = t0 - lt - t1 = Timestamp.from_mjd([t1], scale="tdb") - orbit_propagated = self.propagate_orbits(orbit, t1) + # Update the previous light travel time to this iteration's light travel time + lt_prev = lt - orbit_i = orbit_propagated - t0 = t1 - lt0 = lt - dlt = dlt - orbits_aberrated = qv.concatenate([orbits_aberrated, orbit]) + orbits_aberrated = qv.concatenate([orbits_aberrated, orbit_i]) lts[i] = lt return orbits_aberrated, lts + + def _generate_ephemeris( self, orbits: OrbitType, observers: ObserverType, lt_tol: float = 1e-10 ) -> EphemerisType: @@ -175,7 +178,6 @@ def _generate_ephemeris( origin_out=OriginCodes.SOLAR_SYSTEM_BARYCENTER, ), ) - num_orbits = len(propagated_orbits_barycentric.orbit_id.unique()) observer_codes = np.tile( @@ -202,7 +204,9 @@ def _generate_ephemeris( vz=propagated_orbits_aberrated.coordinates.values[:, 5] - observers_barycentric.coordinates.values[:, 5], covariance=None, - time=propagated_orbits_aberrated.coordinates.time, + # The ephemeris times are at the point of the observer, + # not the aberated orbit + time=observers.coordinates.time, origin=Origin.from_kwargs(code=observer_codes), frame="ecliptic", ) @@ -210,6 +214,7 @@ def _generate_ephemeris( spherical_coordinates = SphericalCoordinates.from_cartesian( topocentric_coordinates ) + light_time = np.array(light_time) spherical_coordinates = transform_coordinates( @@ -379,7 +384,7 @@ def generate_ephemeris( ephemeris_variants = None else: - ephemeris = self._generate_ephemeris(orbits, observers) + ephemeris = self._generate_ephemeris(orbits, observers, lt_tol=1e-20) if covariance is True and not orbits.coordinates.covariance.is_all_nan(): variants = VariantOrbits.create( diff --git a/src/adam_core/time/time.py b/src/adam_core/time/time.py index c5d57f03..4f18f425 100644 --- a/src/adam_core/time/time.py +++ b/src/adam_core/time/time.py @@ -392,8 +392,7 @@ def add_fractional_days( nano_part = pc.subtract(fractional_days, day_part) days = pc.cast(day_part, pa.int64()) - nanos = pc.cast(pc.multiply(nano_part, 86400 * 1e9), pa.int64()) - + nanos = pc.cast(pc.multiply(nano_part, 86400 * 1e9), options=pc.CastOptions(target_type=pa.int64(), allow_float_truncate=True)) return self.add_days(days).add_nanos(nanos) def difference_scalar( From 0a07b7b43a6dcde90d658c8304be4136228e3ecd Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 1 Aug 2024 11:56:34 -0400 Subject: [PATCH 06/11] linting --- src/adam_core/_version.py | 2 +- src/adam_core/coordinates/origin.py | 4 +++- src/adam_core/propagator/propagator.py | 23 +++++++++++------------ src/adam_core/time/time.py | 5 ++++- 4 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/adam_core/_version.py b/src/adam_core/_version.py index 3261ed66..d7bee229 100644 --- a/src/adam_core/_version.py +++ b/src/adam_core/_version.py @@ -1 +1 @@ -__version__ = '0.2.1.dev5+g4d3da0c.d20240731' \ No newline at end of file +__version__ = "0.2.1.dev5+g4d3da0c.d20240731" diff --git a/src/adam_core/coordinates/origin.py b/src/adam_core/coordinates/origin.py index b75f0702..fa7231c1 100644 --- a/src/adam_core/coordinates/origin.py +++ b/src/adam_core/coordinates/origin.py @@ -95,7 +95,9 @@ def as_OriginCodes(self) -> OriginCodes: OriginCodes Origin codes as an `~adam_core.coordinates.origin.OriginCodes` object. """ - assert len(self.code.unique()) == 1, "Only one origin code can be converted at a time." + assert ( + len(self.code.unique()) == 1 + ), "Only one origin code can be converted at a time." return OriginCodes[self.code.unique()[0].as_py()] def __eq__(self, other: object) -> np.ndarray: diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 4aa63261..92385282 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -94,18 +94,18 @@ class EphemerisMixin: """ def _add_light_time( - self, - orbits, - observers, - lt_tol: float = 1e-12, - max_iter: int = 10, - ): + self, + orbits, + observers, + lt_tol: float = 1e-12, + max_iter: int = 10, + ): orbits_aberrated = orbits.empty() lts = np.zeros(len(orbits)) for i, (orbit, observer) in enumerate(zip(orbits, observers)): # Set the running variables lt_prev = 0 - dlt = float('inf') + dlt = float("inf") orbit_i = orbit lt = 0 @@ -117,7 +117,6 @@ def _add_light_time( # the light travel time will be calculated) t0 = orbit_i.coordinates.time.rescale("tdb").mjd()[0].as_py() - iterations = 0 while dlt > lt_tol and iterations < max_iter: iterations += 1 @@ -132,7 +131,9 @@ def _add_light_time( dlt = np.abs(lt - lt_prev) # Calculate the new epoch and propagate the initial orbit to that epoch - orbit_i = self.propagate_orbits(orbit, Timestamp.from_mjd([t0 - lt], scale="tdb")) + orbit_i = self.propagate_orbits( + orbit, Timestamp.from_mjd([t0 - lt], scale="tdb") + ) # Update the previous light travel time to this iteration's light travel time lt_prev = lt @@ -142,8 +143,6 @@ def _add_light_time( return orbits_aberrated, lts - - def _generate_ephemeris( self, orbits: OrbitType, observers: ObserverType, lt_tol: float = 1e-10 ) -> EphemerisType: @@ -214,7 +213,7 @@ def _generate_ephemeris( spherical_coordinates = SphericalCoordinates.from_cartesian( topocentric_coordinates ) - + light_time = np.array(light_time) spherical_coordinates = transform_coordinates( diff --git a/src/adam_core/time/time.py b/src/adam_core/time/time.py index 4f18f425..beabc283 100644 --- a/src/adam_core/time/time.py +++ b/src/adam_core/time/time.py @@ -392,7 +392,10 @@ def add_fractional_days( nano_part = pc.subtract(fractional_days, day_part) days = pc.cast(day_part, pa.int64()) - nanos = pc.cast(pc.multiply(nano_part, 86400 * 1e9), options=pc.CastOptions(target_type=pa.int64(), allow_float_truncate=True)) + nanos = pc.cast( + pc.multiply(nano_part, 86400 * 1e9), + options=pc.CastOptions(target_type=pa.int64(), allow_float_truncate=True), + ) return self.add_days(days).add_nanos(nanos) def difference_scalar( From e7e766d8c2a8b5057061e4974751cb6ca1c42f7b Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 1 Aug 2024 12:59:26 -0400 Subject: [PATCH 07/11] Update src/adam_core/propagator/propagator.py Co-authored-by: Joachim Moeyens --- src/adam_core/propagator/propagator.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 92385282..446db155 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -10,7 +10,9 @@ from ..coordinates.origin import Origin, OriginCodes from ..coordinates.spherical import SphericalCoordinates from ..coordinates.transform import transform_coordinates -from ..dynamics.aberrations import C +from ..constants import Constants as c + +C = c.C from ..observers.observers import Observers from ..orbits.ephemeris import Ephemeris from ..orbits.orbits import Orbits From 4289ad29e7334612d30df73f595d407af1230165 Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 1 Aug 2024 13:23:52 -0400 Subject: [PATCH 08/11] Update src/adam_core/propagator/propagator.py Co-authored-by: Joachim Moeyens --- src/adam_core/propagator/propagator.py | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 446db155..124d2386 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -191,22 +191,17 @@ def _generate_ephemeris( lt_tol=lt_tol, ) + topocentric_state = propagated_orbits_aberrated.coordinates.values - observers_barycentric.coordinates.values topocentric_coordinates = CartesianCoordinates.from_kwargs( - x=propagated_orbits_aberrated.coordinates.values[:, 0] - - observers_barycentric.coordinates.values[:, 0], - y=propagated_orbits_aberrated.coordinates.values[:, 1] - - observers_barycentric.coordinates.values[:, 1], - z=propagated_orbits_aberrated.coordinates.values[:, 2] - - observers_barycentric.coordinates.values[:, 2], - vx=propagated_orbits_aberrated.coordinates.values[:, 3] - - observers_barycentric.coordinates.values[:, 3], - vy=propagated_orbits_aberrated.coordinates.values[:, 4] - - observers_barycentric.coordinates.values[:, 4], - vz=propagated_orbits_aberrated.coordinates.values[:, 5] - - observers_barycentric.coordinates.values[:, 5], + x=topocentric_state[:, 0], + y=topocentric_state[:, 1], + z=topocentric_state[:, 2], + vx=topocentric_state[:, 3], + vy=topocentric_state[:, 4], + vz=topocentric_state[:, 5], covariance=None, # The ephemeris times are at the point of the observer, - # not the aberated orbit + # not the aberrated orbit time=observers.coordinates.time, origin=Origin.from_kwargs(code=observer_codes), frame="ecliptic", From 6708ee24df5d8a7aa14f0dfde69c3d15900e568f Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 1 Aug 2024 13:46:03 -0400 Subject: [PATCH 09/11] Push lt_tol down to light time method --- src/adam_core/propagator/propagator.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 124d2386..2eb2c6de 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -6,11 +6,11 @@ import numpy.typing as npt import quivr as qv +from ..constants import Constants as c from ..coordinates.cartesian import CartesianCoordinates from ..coordinates.origin import Origin, OriginCodes from ..coordinates.spherical import SphericalCoordinates from ..coordinates.transform import transform_coordinates -from ..constants import Constants as c C = c.C from ..observers.observers import Observers @@ -146,7 +146,7 @@ def _add_light_time( return orbits_aberrated, lts def _generate_ephemeris( - self, orbits: OrbitType, observers: ObserverType, lt_tol: float = 1e-10 + self, orbits: OrbitType, observers: ObserverType ) -> EphemerisType: """ A generic ephemeris implementation, which can be used or overridden by subclasses. @@ -188,10 +188,13 @@ def _generate_ephemeris( propagated_orbits_aberrated, light_time = self._add_light_time( propagated_orbits_barycentric, observers_barycentric, - lt_tol=lt_tol, + lt_tol=1e-12, ) - topocentric_state = propagated_orbits_aberrated.coordinates.values - observers_barycentric.coordinates.values + topocentric_state = ( + propagated_orbits_aberrated.coordinates.values + - observers_barycentric.coordinates.values + ) topocentric_coordinates = CartesianCoordinates.from_kwargs( x=topocentric_state[:, 0], y=topocentric_state[:, 1], @@ -380,7 +383,7 @@ def generate_ephemeris( ephemeris_variants = None else: - ephemeris = self._generate_ephemeris(orbits, observers, lt_tol=1e-20) + ephemeris = self._generate_ephemeris(orbits, observers) if covariance is True and not orbits.coordinates.covariance.is_all_nan(): variants = VariantOrbits.create( From 1988db5f5dd47b39acf1f04ae2db623bf4fdb1aa Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 1 Aug 2024 13:51:24 -0400 Subject: [PATCH 10/11] linting seemed to care about this --- src/adam_core/propagator/propagator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 2eb2c6de..69443045 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -11,8 +11,6 @@ from ..coordinates.origin import Origin, OriginCodes from ..coordinates.spherical import SphericalCoordinates from ..coordinates.transform import transform_coordinates - -C = c.C from ..observers.observers import Observers from ..orbits.ephemeris import Ephemeris from ..orbits.orbits import Orbits @@ -23,6 +21,8 @@ logger = logging.getLogger(__name__) +C = c.C + RAY_INSTALLED = False try: import ray From ac69b6a45a6713dc613e38b2b9bf63bea491c3dc Mon Sep 17 00:00:00 2001 From: Alec Koumjian Date: Thu, 1 Aug 2024 14:07:24 -0400 Subject: [PATCH 11/11] When dealing with ephemeris, always provide in UTC --- src/adam_core/orbits/query/horizons.py | 5 +---- src/adam_core/propagator/propagator.py | 6 ++++++ 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/adam_core/orbits/query/horizons.py b/src/adam_core/orbits/query/horizons.py index 9bbe4866..c4091d0e 100644 --- a/src/adam_core/orbits/query/horizons.py +++ b/src/adam_core/orbits/query/horizons.py @@ -226,9 +226,6 @@ def query_horizons_ephemeris( ignore_index=True, ) - # Horizons produces UTC but we use tdb everywhere - epochs = Timestamp.from_jd(pa.array(dfs["datetime_jd"]), scale="utc").rescale("tdb") - ephemeris = Ephemeris.from_kwargs( orbit_id=dfs["orbit_id"], object_id=dfs["targetname"], @@ -236,7 +233,7 @@ def query_horizons_ephemeris( light_time=dfs["lighttime"] / 1440, alpha=dfs["alpha"], coordinates=SphericalCoordinates.from_kwargs( - time=epochs, + time=Timestamp.from_jd(pa.array(dfs["datetime_jd"]), scale="utc"), lon=dfs["RA"], lat=dfs["DEC"], origin=Origin.from_kwargs(code=dfs["observatory_code"]), diff --git a/src/adam_core/propagator/propagator.py b/src/adam_core/propagator/propagator.py index 69443045..713438b4 100644 --- a/src/adam_core/propagator/propagator.py +++ b/src/adam_core/propagator/propagator.py @@ -220,6 +220,12 @@ def _generate_ephemeris( spherical_coordinates, SphericalCoordinates, frame_out="equatorial" ) + # Ephemeris are generally compared in UTC, so rescale the time + spherical_coordinates = spherical_coordinates.set_column( + "time", + spherical_coordinates.time.rescale("utc"), + ) + if isinstance(orbits, Orbits): ephemeris = Ephemeris.from_kwargs(