Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
9 changes: 5 additions & 4 deletions imap_processing/glows/l1b/glows_l1b_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -953,18 +953,19 @@ def update_spice_parameters(self, spin_offset_correction: float = 0.0) -> None:

# Calculate spin axis orientation

spin_axis_all_times = geometry.cartesian_to_latitudinal(
spin_axis_all_times = geometry.cartesian_to_spherical(
geometry.frame_transform(
time_range,
np.array([0, 0, 1]),
SpiceFrame.IMAP_SPACECRAFT,
SpiceFrame.ECLIPJ2000,
),
degrees=False,
)
# Calculate circular statistics for longitude (wraps around)
lon_mean = circmean(spin_axis_all_times[..., 1], low=-np.pi, high=np.pi)
lon_std = circstd(spin_axis_all_times[..., 1], low=-np.pi, high=np.pi)
# Convert to radians for use with circular statistics methods
spin_axis_all_times = np.radians(spin_axis_all_times)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

is converting to radian ok for other two variables that haven't changed in this work? Eg. lat_mean and lat_std

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Yes. Just below on lines 972-973, both latitude and longitude get converted back to degrees. The scipy.circmean and scipy.circstd functions require that the inputs be in radians.

lon_mean = circmean(spin_axis_all_times[..., 1], low=0, high=2 * np.pi)
lon_std = circstd(spin_axis_all_times[..., 1], low=0, high=2 * np.pi)
lat_mean = circmean(spin_axis_all_times[..., 2], low=-np.pi, high=np.pi)
lat_std = circstd(spin_axis_all_times[..., 2], low=-np.pi, high=np.pi)
# Convert circular statistics to degrees and store
Expand Down
26 changes: 16 additions & 10 deletions imap_processing/glows/l2/glows_l2_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import numpy as np
import xarray as xr
from numpy.typing import NDArray
from scipy.stats import circmean, circstd

from imap_processing.glows import FLAG_LENGTH
from imap_processing.glows.l1b.glows_l1b_data import PipelineSettings
Expand Down Expand Up @@ -470,16 +471,21 @@ def __init__(
.std(dim="epoch")
.data[np.newaxis, :]
)
self.spin_axis_orientation_average = (
good_data["spin_axis_orientation_average"]
.mean(dim="epoch")
.data[np.newaxis, :]
)
self.spin_axis_orientation_std_dev = (
good_data["spin_axis_orientation_average"]
.std(dim="epoch")
.data[np.newaxis, :]
)
spin_axis_data = good_data[
"spin_axis_orientation_average"
].data # (n_epochs, 2)
if spin_axis_data.shape[0] > 0:
# Use circular statistics for longitude since it will be near the 0->360
# boundary. Latitude will never be near the pole, so standard mean and
# std functions are appropriate
lon_avg = circmean(np.radians(spin_axis_data[:, 0]), low=0, high=2 * np.pi)
lon_std = circstd(np.radians(spin_axis_data[:, 0]), low=0, high=360)
lat_avg = float(np.mean(spin_axis_data[:, 1]))
lat_std = float(np.std(spin_axis_data[:, 1]))
else:
lon_avg = lon_std = lat_avg = lat_std = np.nan
self.spin_axis_orientation_average = np.array([[np.degrees(lon_avg), lat_avg]])
self.spin_axis_orientation_std_dev = np.array([[np.degrees(lon_std), lat_std]])

# Select calibration factor corresponding to the mid-epoch in the L1B data.
if len(good_data["epoch"].data) != 0:
Expand Down
45 changes: 18 additions & 27 deletions imap_processing/tests/glows/test_glows_l1b_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -342,20 +342,17 @@ def test_update_spice_parameters_spin_axis_near_wrapping_point(
mock_get_instrument_spin_phase,
mock_imap_state,
):
"""Test spin axis orientation calculation near the longitude wrapping point.
"""Test spin axis orientation longitude is in [0, 360] near the 0/360 boundary.

This test verifies that cartesian_to_latitudinal is called with degrees=False
so that circmean/circstd receive values in radians. The bug was that without
degrees=False, values in degrees would be passed to circmean/circstd which
expect radians with low=-pi, high=pi.
Verifies that circmean correctly averages longitudes that straddle the 0/360
degree boundary and that the result is wrapped to [0, 360) rather than [-180, 180].

Test conditions:
- Longitude values straddling +/-pi (wrapping point at 180 degrees)
- Longitude values straddling 0/360 degrees
- Latitude near equator (-4 degrees)

If the bug existed (degrees=True or default), circmean would receive values
like 179 or -179 degrees when it expects radians in [-pi, pi]. This would
produce nonsensical results because 179 >> pi.
Arithmetic mean of [359, 1, 358, 2, 0] would give ~144 degrees (wrong).
Correct circular mean should give ~0 degrees.
"""
# Mock time conversions - creates a time range of 5 seconds
mock_met_to_sclkticks.return_value = 1000
Expand All @@ -371,19 +368,17 @@ def test_update_spice_parameters_spin_axis_near_wrapping_point(
# Mock instrument spin phase
mock_get_instrument_spin_phase.return_value = 0.5

# Create cartesian vectors that straddle the longitude wrapping point.
# Some points at +179 degrees and some at -179 degrees (which should
# average to ~180 degrees when using proper circular mean).
# Latitude near equator at -4 degrees.
# Create cartesian vectors that straddle the 0/360 degree boundary.
# Points at [359, 1, 358, 2, 0] degrees longitude, latitude near equator.
#
# For spherical to cartesian (r=1):
# x = cos(lat) * cos(lon)
# y = cos(lat) * sin(lon)
# z = sin(lat)
lat_rad = np.deg2rad(-4.0) # Near equator

# Create 5 time steps: [+179, -179, +178, -178, +180] degrees longitude
longitudes_deg = np.array([179.0, -179.0, 178.0, -178.0, 180.0])
# Equivalent in [-180, 180] range: [-1, 1, -2, 2, 0] degrees
longitudes_deg = np.array([359.0, 1.0, 358.0, 2.0, 0.0])
longitudes_rad = np.deg2rad(longitudes_deg)

# Build cartesian vectors for each time step
Expand Down Expand Up @@ -424,18 +419,14 @@ def __init__(self):
lon_std = mock_hist.spin_axis_orientation_std_dev[0]
lat_std = mock_hist.spin_axis_orientation_std_dev[1]

# The circular mean of [179, -179, 178, -178, 180] should be ~180 degrees
# (or equivalently -180 degrees). The key test is that the result is NOT
# near 0 degrees, which would happen if the values weren't properly handled
# as circular data near the wrapping point.
#
# With the bug (degrees=True), circmean would receive [179, -179, ...] and
# interpret these as radians, giving completely wrong results.
#
# Check that longitude is near 180 degrees (could be reported as -180)
assert abs(abs(lon_result) - 180.0) < 5.0, (
f"Longitude {lon_result} should be near +/-180 degrees. "
"If near 0, the circular mean failed at the wrapping point."
# Output longitude must be in [0, 360).
assert 0.0 <= lon_result < 360.0, f"Longitude {lon_result} is outside [0, 360)."

# Circular mean of values near 0/360 boundary should be near 0 degrees,
# not near 144 degrees (arithmetic mean) or 180 degrees.
assert lon_result < 5.0 or lon_result > 355.0, (
f"Longitude {lon_result} should be near 0/360 degrees. "
"If near 144 or 180, the circular mean failed at the wrapping point."
)

# Latitude should be near -4 degrees
Expand Down
42 changes: 42 additions & 0 deletions imap_processing/tests/glows/test_glows_l2_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -521,3 +521,45 @@ def test_position_angle_offset_average(
l2 = HistogramL2(l1b_dataset_full, pipeline_settings, mock_calibration_dataset)

assert l2.position_angle_offset_average == pytest.approx(mock_pa)


# ── spin_axis_orientation_average tests ──────────────────────────────────────


def test_spin_axis_orientation_average_wrapping(
l1b_dataset_full,
pipeline_settings,
mock_ecliptic_bin_centers,
mock_calibration_dataset,
):
"""Longitude averaging uses circular mean to handle the 0/360 boundary correctly.

Arithmetic mean of [358, 2] gives 180 (wrong).
Circular mean correctly gives ~0.
"""
wrapping_dataset = l1b_dataset_full.copy()
wrapping_dataset["spin_axis_orientation_average"] = xr.DataArray(
np.array([[358.0, 5.0], [2.0, 5.0]]),
dims=["epoch", "lonlat"],
)

with (
patch.object(HistogramL2, "compute_position_angle", return_value=42.5),
patch.object(HistogramL2, "get_calibration_factor", return_value=1.0),
):
l2 = HistogramL2(wrapping_dataset, pipeline_settings, mock_calibration_dataset)

lon_avg = l2.spin_axis_orientation_average[0, 0]
lat_avg = l2.spin_axis_orientation_average[0, 1]

# Result must be in [0, 360).
assert 0.0 <= lon_avg < 360.0, f"Longitude {lon_avg} is outside [0, 360)."

# Circular mean of 358 and 2 should be near 0 (both within 2° of boundary),
# not 180 (arithmetic mean).
assert lon_avg < 5.0 or lon_avg > 355.0, (
f"Longitude {lon_avg} should be near 0/360, not near 180."
)

# Latitude is not circular — arithmetic mean of [5, 5] = 5.
assert lat_avg == pytest.approx(5.0)
Loading