Working test for get_hpc_info

stixpy/coordinates/tests/test_transforms.py

@@ -5,7 +5,6 @@
 import pytest
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.time import Time
-from numpy.testing import assert_allclose
 from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
 
 from stixpy.coordinates.frames import STIXImaging
@@ -63,34 +62,35 @@ def test_get_aux_data():
 @pytest.mark.remote_data
 def test_get_hpc_info():
     with pytest.warns(match="SAS solution not available.*"):
-        roll, solo_heeq, stix_pointing = get_hpc_info(Time(["2022-08-28T16:00:16", "2022-08-28T16:00:17.000"]))
+        roll, solo_heeq, stix_pointing = get_hpc_info(Time(["2022-08-28T16:00:17", "2022-08-28T16:00:18.000"]))
 
-    assert_quantity_allclose(-3.3733306 * u.deg, roll)
-    assert_quantity_allclose([25.61525646, 57.914266] * u.arcsec, stix_pointing)
-    assert_quantity_allclose([-97868816.0, 62984852.0, -5531986.5] * u.km, solo_heeq)
+    assert_quantity_allclose(roll, -3.3733292 * u.deg)
+    assert_quantity_allclose(stix_pointing[0, :], [25.622768, 57.920166] * u.arcsec)
+    assert_quantity_allclose(solo_heeq[0, :], [-97868803.82, 62984839.12, -5531987.445] * u.km)
 
     roll, solo_heeq, stix_pointing = get_hpc_info(Time("2022-08-28T16:00:00"), end_time=Time("2022-08-28T16:03:59.575"))
-    assert_quantity_allclose(-3.3732104 * u.deg, roll)
-    assert_quantity_allclose([25.923784, 58.179676] * u.arcsec, stix_pointing)
-    assert_quantity_allclose([-9.7867768e07, 6.2983744e07, -5532067.5] * u.km, solo_heeq)
+    assert_quantity_allclose(
+        roll,
+        -3.3732104 * u.deg,
+    )
+    assert_quantity_allclose(
+        stix_pointing[0, :],
+        [25.923784, 58.179676] * u.arcsec,
+    )
+    assert_quantity_allclose(solo_heeq, [-9.7867768e07, 6.2983744e07, -5532067.5] * u.km)
 
     roll, solo_heeq, stix_pointing = get_hpc_info(Time("2022-08-28T21:00:00"), end_time=Time("2022-08-28T21:03:59.575"))
-    assert_quantity_allclose(-3.3619127 * u.deg, roll)
-    assert_quantity_allclose([-26.070198, 173.48871] * u.arcsec, stix_pointing)
+    assert_quantity_allclose(roll, -3.3619127 * u.deg)
+    assert_quantity_allclose(stix_pointing[0, :], [-26.070198, 173.48871] * u.arcsec)
     assert_quantity_allclose([-9.7671984e07, 6.2774768e07, -5547166.0] * u.km, solo_heeq)
 
 
 @pytest.mark.remote_data
 def test_get_hpc_info_shapes():
     t = Time("2022-08-28T16:00:00") + np.arange(10) * u.min
     roll1, solo_heeq1, stix_pointing1 = get_hpc_info(t)
+    roll2, solo_heeq2, stix_pointing2 = get_hpc_info(t[5])
 
-    roll2, solo_heeq2, stix_pointing2 = get_hpc_info(t[[0, -1]])
-
-    roll3, solo_heeq3, stix_pointing3 = get_hpc_info(t[5])
-
-    assert_allclose(np.mean(roll1, axis=0), roll2, rtol=1e-6)
-    assert_allclose(np.mean(roll1, axis=0), roll3, rtol=5e-5)
-
-    assert_allclose(np.mean(solo_heeq1, axis=0), solo_heeq2.value, rtol=1e-6)
-    assert_allclose(np.mean(solo_heeq1, axis=0), solo_heeq3[0, :], rtol=1e-5)
+    assert_quantity_allclose(roll1[5], roll2)
+    assert_quantity_allclose(solo_heeq1[5, :], solo_heeq2[0, :])
+    assert_quantity_allclose(stix_pointing1[5, :], stix_pointing2[0, :])

stixpy/coordinates/transforms.py

@@ -55,7 +55,7 @@ def _get_rotation_matrix_and_position(obstime):
     return rmatrix, solo_position_heeq
 
 
-def get_hpc_info(times):
+def get_hpc_info(times, end_time=None):
     r"""
     Get STIX pointing and SO location from L2 aspect files.
 
@@ -71,9 +71,11 @@ def get_hpc_info(times):
     aux = _get_aux_data(times.min(), times.max())
 
     indices = np.argwhere((aux["time"] >= times.min()) & (aux["time"] <= times.max()))
+    if end_time is not None:
+        indices = np.argwhere((aux["time"] >= times.min()) & (aux["time"] <= end_time))
     indices = indices.flatten()
 
-    if times.size == 2 and indices.size >= 2:
+    if end_time is not None and times.size == 1 and indices.size >= 2:
         # mean
         aux = aux[indices]
 
@@ -84,7 +86,7 @@ def get_hpc_info(times):
         good_sas = aux[good_solution]
 
         if len(good_sas) == 0:
-            warnings.warn(f"No good SAS solution found for time range: {times[0]} to {times[1]}.")
+            warnings.warn(f"No good SAS solution found for time range: {times} to {end_time}.")
             sas_x = 0
             sas_y = 0
         else:
@@ -97,7 +99,7 @@ def get_hpc_info(times):
                 warnings.warn(f"Pointing unstable: StD(X) = {sigma_x}, StD(Y) = {sigma_y}.")
     else:
         if indices.size < 2:
-            logger.info("Only one data point found interpolating between two closest times.")
+            logger.info("Only one data contained in time interval found interpolating between two closest times.")
             # Times contained in one time or only contains one time
             if times.size == 2:
                 times = times[0] + (times[1] - times[0]) * 0.5
@@ -113,13 +115,13 @@ def get_hpc_info(times):
             aux = aux[start_ind:end_ind]
             indices = [0, 1]
 
-        # Interpolate
+        # Interpolate all times
         x = (times - aux["time"][0]).to_value(u.s)
         xp = (aux["time"] - aux["time"][0]).to_value(u.s)
 
         roll = np.interp(x, xp, aux["roll_angle_rpy"][:, 0].value) << aux["roll_angle_rpy"].unit
-        yaw = np.interp(x, xp, aux["roll_angle_rpy"][:, 1].value) << aux["roll_angle_rpy"].unit
-        pitch = np.interp(x, xp, aux["roll_angle_rpy"][:, 2].value) << aux["roll_angle_rpy"].unit
+        pitch = np.interp(x, xp, aux["roll_angle_rpy"][:, 1].value) << aux["roll_angle_rpy"].unit
+        yaw = np.interp(x, xp, aux["roll_angle_rpy"][:, 2].value) << aux["roll_angle_rpy"].unit
 
         solo_heeq = (
             np.vstack(
@@ -177,27 +179,27 @@ def get_hpc_info(times):
     # Convert the spacecraft pointing to STIX frame
     rotated_yaw = -yaw * np.cos(roll) + pitch * np.sin(roll)
     rotated_pitch = yaw * np.sin(roll) + pitch * np.cos(roll)
-    spacecraft_pointing = np.hstack([STIX_X_SHIFT + rotated_yaw, STIX_Y_SHIFT + rotated_pitch])
+    spacecraft_pointing = np.vstack([STIX_X_SHIFT + rotated_yaw, STIX_Y_SHIFT + rotated_pitch])
     stix_pointing = spacecraft_pointing
-    sas_pointing = np.hstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET])
+    sas_pointing = np.vstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET])
 
-    pointing_diff = np.linalg.norm(spacecraft_pointing - sas_pointing)
+    pointing_diff = np.sqrt(np.sum((spacecraft_pointing - sas_pointing) ** 2, axis=0))
     if np.all(np.isfinite(sas_pointing)) and len(good_sas) > 0:
-        if pointing_diff < 200 * u.arcsec:
+        if np.max(pointing_diff) < 200 * u.arcsec:
             logger.info(f"Using SAS pointing: {sas_pointing}")
-            stix_pointing = sas_pointing
+            stix_pointing = np.where(pointing_diff < 200 * u.arcsec, sas_pointing, spacecraft_pointing)
         else:
             warnings.warn(
                 f"Using spacecraft pointing: {spacecraft_pointing}" f" large difference between SAS and spacecraft."
             )
     else:
         warnings.warn(f"SAS solution not available using spacecraft pointing: {stix_pointing}.")
 
-    return roll, solo_heeq, stix_pointing
+    return roll, solo_heeq, stix_pointing.T
 
 
 @lru_cache
-def _get_aux_data(start_time, end_time):
+def _get_aux_data(start_time, end_time=None):
     r"""
     Search, download and read L2 pointing data.
 
@@ -210,6 +212,8 @@ def _get_aux_data(start_time, end_time):
     -------
 
     """
+    if end_time is None:
+        end_time = start_time
     # Find, download, read aux file with pointing, sas and position information
     logger.debug(f"Searching for AUX data: {start_time} - {end_time}")
     query = Fido.search(