diff --git a/imap_processing/glows/l1b/glows_l1b_data.py b/imap_processing/glows/l1b/glows_l1b_data.py index 28e570f32..7fcb18331 100644 --- a/imap_processing/glows/l1b/glows_l1b_data.py +++ b/imap_processing/glows/l1b/glows_l1b_data.py @@ -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) + 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 diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index 8cf3a6c4e..7cbfeb626 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -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 @@ -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: diff --git a/imap_processing/tests/glows/test_glows_l1b_data.py b/imap_processing/tests/glows/test_glows_l1b_data.py index 3d5596279..bc298c468 100644 --- a/imap_processing/tests/glows/test_glows_l1b_data.py +++ b/imap_processing/tests/glows/test_glows_l1b_data.py @@ -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 @@ -371,10 +368,8 @@ 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) @@ -382,8 +377,8 @@ def test_update_spice_parameters_spin_axis_near_wrapping_point( # 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 @@ -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 diff --git a/imap_processing/tests/glows/test_glows_l2_data.py b/imap_processing/tests/glows/test_glows_l2_data.py index bef47d06f..a34d978ba 100644 --- a/imap_processing/tests/glows/test_glows_l2_data.py +++ b/imap_processing/tests/glows/test_glows_l2_data.py @@ -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)