Merge pull request #138 from sot/planets-performance

Performance improvements for planets
sot · Nov 15, 2023 · e25f0fc · e25f0fc
2 parents 2258972 + cf5cbc1
commit e25f0fc
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 15 deletions.
diff --git a/chandra_aca/planets.py b/chandra_aca/planets.py
@@ -24,11 +24,10 @@
 from pathlib import Path
 from typing import Optional, Union
 
-import astropy.constants as const
 import astropy.units as u
 import numpy as np
 from astropy.io import ascii
-from cxotime import CxoTime, CxoTimeLike
+from cxotime import CxoTime, CxoTimeLike, convert_time_format
 from ska_helpers.utils import LazyVal
 
 from chandra_aca.transform import eci_to_radec
@@ -134,22 +133,23 @@ def get_planet_angular_sep(
     """
     from agasc import sphere_dist
 
-    if not isinstance(time, CxoTime):
-        time = CxoTime(time)
-
     if observer_position == "earth":
         eci = get_planet_eci(body, time)
         body_ra, body_dec = eci_to_radec(eci)
     elif observer_position == "chandra":
         eci = get_planet_chandra(body, time)
         body_ra, body_dec = eci_to_radec(eci)
     elif observer_position == "chandra-horizons":
-        if time.shape == ():
-            time = CxoTime([time, time + 1000 * u.s])
+        time_secs = convert_time_format(time, "secs")
+
+        if time_secs.shape == ():
+            time_secs = [time_secs, time_secs + 1000]
             is_scalar = True
         else:
             is_scalar = False
-        pos = get_planet_chandra_horizons(body, time[0], time[1], n_times=len(time))
+        pos = get_planet_chandra_horizons(
+            body, time_secs[0], time_secs[-1], n_times=len(time_secs)
+        )
         body_ra = pos["ra"]
         body_dec = pos["dec"]
         if is_scalar:
@@ -191,8 +191,7 @@ def get_planet_barycentric(body: str, time: CxoTimeLike = None):
         )
 
     spk_pairs = BODY_NAME_TO_KERNEL_SPEC[body]
-    time = CxoTime(time)
-    time_jd = time.jd
+    time_jd = convert_time_format(time, "jd")
     pos = kernel[spk_pairs[0]].compute(time_jd)
     for spk_pair in spk_pairs[1:]:
         pos += kernel[spk_pair].compute(time_jd)
@@ -237,16 +236,17 @@ def get_planet_eci(
         Earth-Centered Inertial (ECI) position (km) as (x, y, z)
         or N x (x, y, z)
     """
-    time = CxoTime(time)
+    time_sec = convert_time_format(time, "secs")
 
-    pos_planet = get_planet_barycentric(body, time)
+    pos_planet = get_planet_barycentric(body, time_sec)
     if pos_observer is None:
         pos_observer = get_planet_barycentric("earth", time)
 
-    dist = np.sqrt(np.sum((pos_planet - pos_observer) ** 2, axis=-1)) * u.km
-    light_travel_time = (dist / const.c).to(u.s)
+    dist = np.sqrt(np.sum((pos_planet - pos_observer) ** 2, axis=-1))  # km
+    # Divide distance by the speed of light in km/s
+    light_travel_time = dist / 299792.458  # s
 
-    pos_planet = get_planet_barycentric(body, time - light_travel_time)
+    pos_planet = get_planet_barycentric(body, time_sec - light_travel_time)
 
     return pos_planet - pos_observer
 

diff --git a/chandra_aca/tests/test_maude_decom.py b/chandra_aca/tests/test_maude_decom.py
@@ -319,6 +319,8 @@ def test_vcdu_vs_level0():
             os.path.join(
                 os.path.dirname(__file__), "data", f"test_level0_{slot}.fits.gz"
             ),
+            # Undocumented way to stop UnitsWarning 'DN' did not parse warning. See
+            # https://github.com/astropy/astropy/issues/14330#issuecomment-1406538512
             unit_parse_strict="silent",
         )
         td = l0_test_data[