Update method used to find mean spin-axis for Pointing to just averag… (IMAP-Science-Operations-Center#2533)

* Update method used to find mean spin-axis for Pointing to just average instantaneous z-axes across Pointing

* Update imap_processing/spice/pointing_frame.py

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

* Copilot feedback changes

* Use flight data to verify calculation of mean spin axis

---------

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

Author: tmplummer

imap_processing/spice/pointing_frame.py

@@ -11,7 +11,7 @@
 from imap_data_access import SPICEFilePath
 from numpy.typing import NDArray
 
-from imap_processing.spice.geometry import SpiceFrame
+from imap_processing.spice.geometry import SpiceFrame, frame_transform
 from imap_processing.spice.repoint import get_repoint_data
 from imap_processing.spice.time import (
     TICK_DURATION,
@@ -284,8 +284,8 @@ def calculate_pointing_attitude_segments(
             f"range: ({et_to_utc(pointing_start_et)}, {et_to_utc(pointing_end_et)})"
         )
 
-        # 1 spin/15 seconds; 10 quaternions / spin.
-        num_samples = (pointing_end_et - pointing_start_et) / 15 * 10
+        # Sample at 1Hz
+        num_samples = pointing_end_et - pointing_start_et
         # There were rounding errors when using spiceypy.pxform
         # so np.ceil and np.floor were used to ensure the start
         # and end times were within the ck range.
@@ -295,11 +295,11 @@ def calculate_pointing_attitude_segments(
             int(num_samples),
         )
 
-        # Get the average quaternions for the pointing
-        q_avg = _average_quaternions(et_times)
+        # Get the average spin-axis in HAE coordinates
+        z_avg = _mean_spin_axis(et_times)
 
         # Create a rotation matrix
-        rotation_matrix = _create_rotation_matrix(q_avg)
+        rotation_matrix = _create_rotation_matrix(z_avg)
 
         # Convert the rotation matrix to a quaternion.
         # https://spiceypy.readthedocs.io/en/main/documentation.html#spiceypy.spiceypy.m2q
@@ -317,9 +317,14 @@ def calculate_pointing_attitude_segments(
     return pointing_segments
 
 
-def _average_quaternions(et_times: np.ndarray) -> NDArray:
+def _mean_spin_axis(et_times: np.ndarray) -> NDArray:
     """
-    Average the quaternions.
+    Compute the mean spin axis for a given time range.
+
+    The mean spin-axis is computed by taking the mean of the spacecraft z-axis
+    expressed in HAE Cartesian coordinates at each of the input et_times. The
+    mean is computed by finding the mean of each component of the vector across
+    time.
 
     Parameters
     ----------
@@ -328,72 +333,52 @@ def _average_quaternions(et_times: np.ndarray) -> NDArray:
 
     Returns
     -------
-    q_avg : np.ndarray
-        Average quaternion.
+    z_avg : np.ndarray
+        Mean spin-axis. Shape is (3,), a single 3D vector (x, y, z).
     """
-    aggregate = np.zeros((4, 4))
-    for tdb in et_times:
-        # we use a quick and dirty method here for grabbing the quaternions
-        # from the attitude kernel.  Depending on how well the kernel input
-        # data is built and sampled, there may or may not be aliasing with this
-        # approach.  If it turns out that we need to pull the quaternions
-        # directly from the CK there are several routines that exist to do this
-        # but it's not straight forward.  We'll revisit this if needed.
-
-        # Rotation matrix from IMAP spacecraft frame to ECLIPJ2000.
-        # https://spiceypy.readthedocs.io/en/main/documentation.html#spiceypy.spiceypy.pxform
-        body_rots = spiceypy.pxform("IMAP_SPACECRAFT", "ECLIPJ2000", tdb)
-        # Convert rotation matrix to quaternion.
-        # https://spiceypy.readthedocs.io/en/main/documentation.html#spiceypy.spiceypy.m2q
-        body_quat = spiceypy.m2q(body_rots)
-
-        # Standardize the quaternion so that they may be compared.
-        body_quat = body_quat * np.sign(body_quat[0])
-        # Aggregate quaternions into a single matrix.
-        aggregate += np.outer(body_quat, body_quat)
-
-    # Reference: "On Averaging Rotations".
-    # Link: https://link.springer.com/content/pdf/10.1023/A:1011129215388.pdf
-    aggregate /= len(et_times)
-
-    # Compute eigen values and vectors of the matrix A
-    # Eigenvalues tell you how much "influence" each
-    # direction (eigenvector) has.
-    # The largest eigenvalue corresponds to the direction
-    # that has the most influence.
-    # The eigenvector corresponding to the largest
-    # eigenvalue points in the direction that has the most
-    # combined rotation influence.
-    eigvals, eigvecs = np.linalg.eig(aggregate)
-    # q0: The scalar part of the quaternion.
-    # q1, q2, q3: The vector part of the quaternion.
-    q_avg = eigvecs[:, np.argmax(eigvals)]
-
-    return q_avg
-
-
-def _create_rotation_matrix(q_avg: np.ndarray) -> NDArray:
+    # we use a quick and dirty method here for sampling the instantaneous
+    # spin-axis.  Depending on how well the kernel input
+    # data is built and sampled, there may or may not be aliasing with this
+    # approach.  If it turns out that we need to pull the quaternions
+    # directly from the CK there are several routines that exist to do this
+    # but it's not straight forward.  We'll revisit this if needed.
+    z_inertial_hae = frame_transform(
+        et_times, np.array([0, 0, 1]), SpiceFrame.IMAP_SPACECRAFT, SpiceFrame.ECLIPJ2000
+    )
+
+    # Compute the average spin axis by averaging each component across time
+    z_avg = np.mean(z_inertial_hae, axis=0)
+    # We don't need to worry about the magnitude being close to zero when
+    # normalizing because the instantaneous spin-axes will always be close
+    # to the same direction.
+    z_avg /= np.linalg.norm(z_avg)
+
+    return z_avg
+
+
+def _create_rotation_matrix(z_avg: np.ndarray) -> NDArray:
     """
-    Create a rotation matrix.
+    Create a rotation matrix from the average spin axis.
 
     Parameters
     ----------
-    q_avg : numpy.ndarray
-        Averaged quaternions for the pointing.
+    z_avg : numpy.ndarray
+        Average spin-axis that has been normalized to have unit length expressed
+        in HAE coordinates.
 
     Returns
     -------
     rotation_matrix : np.ndarray
         Rotation matrix.
     """
-    # Converts the averaged quaternion (q_avg) into a rotation matrix
-    # and get inertial z axis.
-    # https://spiceypy.readthedocs.io/en/main/documentation.html#spiceypy.spiceypy.q2m
-    z_avg = spiceypy.q2m(list(q_avg))[:, 2]
-    # y_avg is perpendicular to both z_avg and the standard Z-axis.
+    # y_avg is perpendicular to both z_avg and the HAE Z-axis.
+    # Since z_avg will never point anywhere near the HAE Z-axis, this
+    # cross-product will always work to define the Pointing Y-axis
     y_avg = np.cross(z_avg, [0, 0, 1])
+    y_avg /= np.linalg.norm(y_avg)
     # x_avg is perpendicular to y_avg and z_avg.
     x_avg = np.cross(y_avg, z_avg)
+    x_avg /= np.linalg.norm(x_avg)
 
     # Construct the rotation matrix from x_avg, y_avg, z_avg
     rotation_matrix = np.asarray([x_avg, y_avg, z_avg])

imap_processing/tests/external_test_data_config.py

@@ -138,6 +138,11 @@
     # Lo
     ("imap_lo_l1c_pset_20260101-repoint01261_v001.cdf", "lo/test_cdfs"),
 
+    # Pointing Attitude Kernel
+    ("imap_2025_338_2025_339_001.ah.bc", "spice/test_data/"),
+    ("imap_2025_339_2025_339_001.ah.bc", "spice/test_data/"),
+    ("imap_2025_339_2025_340_001.ah.bc", "spice/test_data/"),
+
     # Ultra
     ("FM90_Startup_20230711T081655.CCSDS", "ultra/data/l0/"),
     ("IMAP-Ultra45_r1_L1_V0_shortened.csv", "ultra/data/l1/"),

imap_processing/tests/spice/test_pointing_frame.py

@@ -10,16 +10,19 @@
 from imap_data_access import SPICEInput
 
 from imap_processing.spice import IMAP_SC_ID
-from imap_processing.spice.geometry import SpiceFrame
+from imap_processing.spice.geometry import (
+    SpiceFrame,
+    spherical_to_cartesian,
+)
 from imap_processing.spice.pointing_frame import (
     POINTING_SEGMENT_DTYPE,
-    _average_quaternions,
     _create_rotation_matrix,
+    _mean_spin_axis,
     calculate_pointing_attitude_segments,
     generate_pointing_attitude_kernel,
     write_pointing_frame_ck,
 )
-from imap_processing.spice.time import TICK_DURATION
+from imap_processing.spice.time import TICK_DURATION, met_to_sclkticks, sct_to_et
 
 
 @pytest.fixture
@@ -36,6 +39,22 @@ def furnish_pointing_frame_kernels(furnish_kernels, spice_test_data_path):
         yield [str(spice_test_data_path / k) for k in required_kernels]
 
 
+@pytest.fixture
+def furnish_flight_ah_kernels(furnish_kernels, spice_test_data_path):
+    """List SPICE kernels."""
+    required_kernels = [
+        "naif0012.tls",
+        "imap_sclk_0000.tsc",
+        "imap_130.tf",
+        "imap_science_100.tf",
+        "imap_2025_338_2025_339_001.ah.bc",
+        "imap_2025_339_2025_339_001.ah.bc",
+        "imap_2025_339_2025_340_001.ah.bc",
+    ]
+    with furnish_kernels(required_kernels):
+        yield [str(spice_test_data_path / k) for k in required_kernels]
+
+
 @pytest.fixture
 def et_times(furnish_pointing_frame_kernels):
     """Tests get_et_times function."""
@@ -165,20 +184,26 @@ def test_write_pointing_frame_ck(
     assert parent_file in lines[5]
 
 
-def test_average_quaternions(et_times, furnish_pointing_frame_kernels):
-    """Tests average_quaternions function."""
-    q_avg = _average_quaternions(et_times)
+@pytest.mark.external_test_data
+def test_mean_spin_axis(furnish_flight_ah_kernels):
+    """Tests _mean_spin_axis function."""
+    # Pointing 69 start/end times as defined in imap_2025_351_01.repoint
+    met_range = np.array([502624925, 502711208])
+    et_range = sct_to_et(met_to_sclkticks(met_range))
+    et_times = np.linspace(et_range[0], et_range[1], int(et_range[1] - et_range[0]))
+    z_avg = _mean_spin_axis(et_times)
 
-    # Generated from MATLAB code results
-    q_avg_expected = np.array([-0.6611, 0.4981, -0.5019, -0.2509])
-    np.testing.assert_allclose(q_avg, q_avg_expected, atol=1e-4)
+    # Generated from GLOWS average spin-axis
+    exp_z_avg_lat = 0.065
+    exp_z_avg_lon = 249.86
+    z_avg_expected = spherical_to_cartesian(np.array([1, exp_z_avg_lon, exp_z_avg_lat]))
+    np.testing.assert_allclose(z_avg, z_avg_expected, atol=1e-4)
 
 
 def test_create_rotation_matrix(et_times, furnish_pointing_frame_kernels):
     """Tests create_rotation_matrix function."""
-    q_avg = _average_quaternions(et_times)
-    rotation_matrix = _create_rotation_matrix(q_avg)
-    z_avg = spiceypy.q2m(list(q_avg))[:, 2]
+    z_avg = _mean_spin_axis(et_times)
+    rotation_matrix = _create_rotation_matrix(z_avg)
 
     rotation_matrix_expected = np.array(
         [