From 9e91de2e16d80a8fac232df6c6b9ae4e89e4ea26 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Tue, 17 Oct 2023 08:26:53 -0700 Subject: [PATCH] Add phase angle (alpha) to ephemeris table Remove variable shadowing in PYOORB._generate_ephemeris while we are at it Co-authored-by: Alec Koumjian --- adam_core/orbits/ephemeris.py | 1 + adam_core/propagator/pyoorb.py | 53 +++++++++++++++++++--------------- 2 files changed, 30 insertions(+), 24 deletions(-) diff --git a/adam_core/orbits/ephemeris.py b/adam_core/orbits/ephemeris.py index 8c3dc2d7..e4e8250e 100644 --- a/adam_core/orbits/ephemeris.py +++ b/adam_core/orbits/ephemeris.py @@ -20,6 +20,7 @@ class Ephemeris(qv.Table): # will be different from its actual geometric coordinates at the time of observation. # Aberrated coordinates are coordinates that account for the light travel time # from the time of emission/reflection to the time of observation + alpha = qv.Float64Column(nullable=True) light_time = qv.Float64Column(nullable=True) aberrated_coordinates = CartesianCoordinates.as_column(nullable=True) diff --git a/adam_core/propagator/pyoorb.py b/adam_core/propagator/pyoorb.py index 19b7cf73..19db4b91 100644 --- a/adam_core/propagator/pyoorb.py +++ b/adam_core/propagator/pyoorb.py @@ -377,7 +377,7 @@ def _generate_ephemeris( epochs_pyoorb = self._configure_epochs(mjd_utc, OpenOrbTimescale.UTC) # Generate ephemeris - ephemeris, err = oo.pyoorb.oorb_ephemeris_full( + ephemeris_array_3D, err = oo.pyoorb.oorb_ephemeris_full( in_orbits=orbits_pyoorb, in_obscode=code_i, in_date_ephems=epochs_pyoorb, @@ -390,7 +390,7 @@ def _generate_ephemeris( ) # Stack 3D ephemeris into single 2D array - ephemeris = np.vstack(ephemeris) + ephemeris_array = np.vstack(ephemeris_array_3D) # PYOORB returns ephemerides for each orbit, so lets reconstruct orbit IDs ids = np.arange(0, len(orbits)) @@ -433,14 +433,14 @@ def _generate_ephemeris( # "obs_z", # "TrueAnom", # ] - codes = np.empty(len(ephemeris), dtype=object) + codes = np.empty(len(ephemeris_array), dtype=object) codes[:] = code_i # Check to make sure the desired times are within 1 microsecond # of the times returned by PYOORB _assert_times_almost_equal( np.tile(mjd_utc, len(orbits)), - ephemeris[:, 0], + ephemeris_array[:, 0], tolerance=0.001, # FIXME: This is 0.001 days, which is 86.4 seconds, not 1 microsecond? ) @@ -459,36 +459,40 @@ def _generate_ephemeris( object_id=object_ids, coordinates=SphericalCoordinates.from_kwargs( time=Timestamp.from_mjd( - ephemeris[:, 0], + ephemeris_array[:, 0], scale="utc", ), rho=None, # PYOORB rho (delta_au) is geocentric not topocentric - lon=ephemeris[:, 1], - lat=ephemeris[:, 2], - vlon=ephemeris[:, 3] / np.cos(np.radians(ephemeris[:, 2])), - vlat=ephemeris[:, 4], + lon=ephemeris_array[:, 1], + lat=ephemeris_array[:, 2], + vlon=ephemeris_array[:, 3] + / np.cos(np.radians(ephemeris_array[:, 2])), + vlat=ephemeris_array[:, 4], vrho=None, # PYOORB doesn't calculate observer velocity so it can't calulate vrho origin=Origin.from_kwargs(code=codes), frame="equatorial", ), + alpha=ephemeris_array[:, 5], aberrated_coordinates=CartesianCoordinates.from_kwargs( - x=ephemeris[:, 24], - y=ephemeris[:, 25], - z=ephemeris[:, 26], - vx=ephemeris[:, 27], - vy=ephemeris[:, 28], - vz=ephemeris[:, 29], + x=ephemeris_array[:, 24], + y=ephemeris_array[:, 25], + z=ephemeris_array[:, 26], + vx=ephemeris_array[:, 27], + vy=ephemeris_array[:, 28], + vz=ephemeris_array[:, 29], time=Timestamp.from_mjd( - ephemeris[:, 0], + ephemeris_array[:, 0], scale="utc", ), origin=Origin.from_kwargs( - code=["SUN" for i in range(len(codes))] + code=["SUN" for i in range(len(ephemeris_array))] ), frame="ecliptic", ), ) + ephemeris_list.append(ephemeris) + elif isinstance(orbits, VariantOrbits): # Map the object and orbit IDs back to the input arrays @@ -501,19 +505,20 @@ def _generate_ephemeris( weights = orbits.weights.to_numpy()[orbit_ids_idx] weights_cov = orbits.weights_cov.to_numpy()[orbit_ids_idx] - ephemeris = VariantEphemeris.from_kwargs( + variant_ephemeris = VariantEphemeris.from_kwargs( orbit_id=orbit_ids, object_id=object_ids, coordinates=SphericalCoordinates.from_kwargs( time=Timestamp.from_mjd( - ephemeris[:, 0], + ephemeris_array[:, 0], scale="utc", ), rho=None, # PYOORB rho (delta_au) is geocentric not topocentric - lon=ephemeris[:, 1], - lat=ephemeris[:, 2], - vlon=ephemeris[:, 3] / np.cos(np.radians(ephemeris[:, 2])), - vlat=ephemeris[:, 4], + lon=ephemeris_array[:, 1], + lat=ephemeris_array[:, 2], + vlon=ephemeris_array[:, 3] + / np.cos(np.radians(ephemeris_array[:, 2])), + vlat=ephemeris_array[:, 4], vrho=None, # PYOORB doesn't calculate observer velocity so it can't calulate vrho origin=Origin.from_kwargs(code=codes), frame="equatorial", @@ -522,6 +527,6 @@ def _generate_ephemeris( weights_cov=weights_cov, ) - ephemeris_list.append(ephemeris) + ephemeris_list.append(variant_ephemeris) return qv.concatenate(ephemeris_list)