Skip to content

Update method used to find mean spin-axis for Pointing to just averag… #2533

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
subagonsouth merged 4 commits into from
Dec 22, 2025
Merged

Update method used to find mean spin-axis for Pointing to just averag… #2533

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7d42222
Update method used to find mean spin-axis for Pointing to just averag…
subagonsouth Dec 22, 2025
fe40f6c
Update imap_processing/spice/pointing_frame.py
subagonsouth Dec 22, 2025
d6273f3
Copilot feedback changes
subagonsouth Dec 22, 2025
a41a28e
Use flight data to verify calculation of mean spin axis
subagonsouth Dec 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
105 changes: 45 additions & 60 deletions imap_processing/spice/pointing_frame.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand All @@ -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
----------
Expand All @@ -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)
Comment on lines +345 to +350
Copy link
Collaborator

@greglucas greglucas Dec 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A couple of potential thoughts/questions here.

What is et_times being requested as? Would it make sense to sample this at a higher frequency under the hood here, to try and get closer to the actual quaternions that were input.

Should we compute the standard deviation and log/warn if it is outside of some bound for awareness purposes?

Copy link
Contributor Author

@subagonsouth subagonsouth Dec 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As it is now, we are sampling ~10 times per spin. IMO, ideally, we should be pulling the quaternions straight from the attitude history files which I think are 1Hz and would give ~15/spin.

I could pull spin 70 data and add a test that checks that we get close to what GLOWS gets.

Copy link
Collaborator

@greglucas greglucas Dec 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IMO, ideally, we should be pulling the quaternions straight from the attitude history files which I think are 1Hz and would give ~15/spin

If you sample the data at 1Hz do we exactly get back the instantaneous value based on the quaternion, or is it somehow smoothed even on the exact data points so that we have a smooth curve but not match the "true" data? i.e. I'm wondering if we can just adjust our sampling period here to be what the quaternions are included at rather than this slight offset to avoid any interpolation/sampling.

# 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])
Expand Down
5 changes: 5 additions & 0 deletions imap_processing/tests/external_test_data_config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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/"),
Expand Down
49 changes: 37 additions & 12 deletions imap_processing/tests/spice/test_pointing_frame.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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."""
Expand Down Expand Up @@ -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(
[
Expand Down