Add phase angle (alpha) to ephemeris table

Remove variable shadowing in PYOORB._generate_ephemeris while we are at it Co-authored-by: Alec Koumjian <[email protected]>
B612-Asteroid-Institute · Oct 17, 2023 · 9e91de2 · 9e91de2
1 parent 2c63f2d
commit 9e91de2
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 24 deletions.
diff --git 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
@@ -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)