From 7d422227dea48be3bea17c152b7ed898905510ce Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 22 Dec 2025 10:30:46 -0700 Subject: [PATCH 1/4] Update method used to find mean spin-axis for Pointing to just average instantaneous z-axes across Pointing --- imap_processing/spice/pointing_frame.py | 90 +++++++------------ .../tests/spice/test_pointing_frame.py | 19 ++-- 2 files changed, 42 insertions(+), 67 deletions(-) diff --git a/imap_processing/spice/pointing_frame.py b/imap_processing/spice/pointing_frame.py index 846f6962e..3346698fd 100644 --- a/imap_processing/spice/pointing_frame.py +++ b/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, @@ -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,12 @@ 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 instantaneous + spin-axes evaluated at each input ET time. Parameters ---------- @@ -328,69 +331,40 @@ def _average_quaternions(et_times: np.ndarray) -> NDArray: Returns ------- - q_avg : np.ndarray - Average quaternion. + z_avg : np.ndarray + Mean spin-axis. Shape is (n, 3) where n is the number of ET times. """ - 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) + + 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 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. y_avg = np.cross(z_avg, [0, 0, 1]) # x_avg is perpendicular to y_avg and z_avg. x_avg = np.cross(y_avg, z_avg) diff --git a/imap_processing/tests/spice/test_pointing_frame.py b/imap_processing/tests/spice/test_pointing_frame.py index 5d7d88e9f..67b8ac6c0 100644 --- a/imap_processing/tests/spice/test_pointing_frame.py +++ b/imap_processing/tests/spice/test_pointing_frame.py @@ -13,8 +13,8 @@ from imap_processing.spice.geometry import SpiceFrame 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, @@ -165,20 +165,21 @@ 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) +def test_mean_spin_axis(et_times, furnish_pointing_frame_kernels): + """Tests _mean_spin_axis function.""" + 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) + z_avg_expected = spiceypy.q2m( + np.array([-0.6611, 0.4981, -0.5019, -0.2509]) + ) @ np.array([0, 0, 1]) + 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( [ From fe40f6c7566689dd3c5da67dc487116cacfde03b Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 22 Dec 2025 10:46:08 -0700 Subject: [PATCH 2/4] Update imap_processing/spice/pointing_frame.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- imap_processing/spice/pointing_frame.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imap_processing/spice/pointing_frame.py b/imap_processing/spice/pointing_frame.py index 3346698fd..0e197f622 100644 --- a/imap_processing/spice/pointing_frame.py +++ b/imap_processing/spice/pointing_frame.py @@ -332,7 +332,7 @@ def _mean_spin_axis(et_times: np.ndarray) -> NDArray: Returns ------- z_avg : np.ndarray - Mean spin-axis. Shape is (n, 3) where n is the number of ET times. + Mean spin-axis. Shape is (3,), a single 3D vector (x, y, z). """ # we use a quick and dirty method here for sampling the instantaneous # spin-axis. Depending on how well the kernel input From d6273f3009d6f994b292480b438fe7ea1f0aadce Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 22 Dec 2025 10:53:49 -0700 Subject: [PATCH 3/4] Copilot feedback changes --- imap_processing/spice/pointing_frame.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/imap_processing/spice/pointing_frame.py b/imap_processing/spice/pointing_frame.py index 0e197f622..ee77b7d12 100644 --- a/imap_processing/spice/pointing_frame.py +++ b/imap_processing/spice/pointing_frame.py @@ -346,6 +346,10 @@ def _mean_spin_axis(et_times: np.ndarray) -> NDArray: # 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 @@ -357,7 +361,8 @@ def _create_rotation_matrix(z_avg: np.ndarray) -> NDArray: Parameters ---------- z_avg : numpy.ndarray - Average spin-axis expressed in HAE coordinates. + Average spin-axis that has been normalized to have unit length expressed + in HAE coordinates. Returns ------- From a41a28ec29a3395f2727671513ce079af0d11116 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 22 Dec 2025 12:42:26 -0700 Subject: [PATCH 4/4] Use flight data to verify calculation of mean spin axis --- imap_processing/spice/pointing_frame.py | 14 +++++-- .../tests/external_test_data_config.py | 5 +++ .../tests/spice/test_pointing_frame.py | 38 +++++++++++++++---- 3 files changed, 46 insertions(+), 11 deletions(-) diff --git a/imap_processing/spice/pointing_frame.py b/imap_processing/spice/pointing_frame.py index ee77b7d12..fbd92dd72 100644 --- a/imap_processing/spice/pointing_frame.py +++ b/imap_processing/spice/pointing_frame.py @@ -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. @@ -321,8 +321,10 @@ def _mean_spin_axis(et_times: np.ndarray) -> NDArray: """ Compute the mean spin axis for a given time range. - The mean spin-axis is computed by taking the mean of the instantaneous - spin-axes evaluated at each input ET time. + 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 ---------- @@ -370,9 +372,13 @@ def _create_rotation_matrix(z_avg: np.ndarray) -> NDArray: Rotation matrix. """ # 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]) diff --git a/imap_processing/tests/external_test_data_config.py b/imap_processing/tests/external_test_data_config.py index 8b2d434c3..6986ac4db 100644 --- a/imap_processing/tests/external_test_data_config.py +++ b/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/"), diff --git a/imap_processing/tests/spice/test_pointing_frame.py b/imap_processing/tests/spice/test_pointing_frame.py index 67b8ac6c0..8725baa4a 100644 --- a/imap_processing/tests/spice/test_pointing_frame.py +++ b/imap_processing/tests/spice/test_pointing_frame.py @@ -10,7 +10,10 @@ 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, _create_rotation_matrix, @@ -19,7 +22,7 @@ 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,14 +184,19 @@ def test_write_pointing_frame_ck( assert parent_file in lines[5] -def test_mean_spin_axis(et_times, furnish_pointing_frame_kernels): +@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 - z_avg_expected = spiceypy.q2m( - np.array([-0.6611, 0.4981, -0.5019, -0.2509]) - ) @ np.array([0, 0, 1]) + # 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)