From 25d7488f9ec60a98b9351e11bfa4d8639931e840 Mon Sep 17 00:00:00 2001 From: Maxine Hartnett Date: Tue, 17 Mar 2026 14:42:40 -0600 Subject: [PATCH 1/4] checkpoint commit --- imap_processing/glows/l2/glows_l2.py | 17 +++++++++++---- imap_processing/glows/l2/glows_l2_data.py | 26 ++++++++++++++++------- 2 files changed, 31 insertions(+), 12 deletions(-) diff --git a/imap_processing/glows/l2/glows_l2.py b/imap_processing/glows/l2/glows_l2.py index 29aa589d0..9eebfc12b 100644 --- a/imap_processing/glows/l2/glows_l2.py +++ b/imap_processing/glows/l2/glows_l2.py @@ -11,6 +11,7 @@ PipelineSettings, ) from imap_processing.glows.l2.glows_l2_data import HistogramL2 +from imap_processing.glows.utils.constants import GlowsConstants from imap_processing.spice.time import et_to_datetime64, ttj2000ns_to_et @@ -85,7 +86,7 @@ def create_l2_dataset( ) bins = xr.DataArray( - np.arange(histogram_l2.daily_lightcurve.number_of_bins, dtype=np.uint32), + np.arange(GlowsConstants.STANDARD_BIN_COUNT, dtype=np.uint32), name="bins", dims=["bins"], attrs=attrs.get_variable_attributes("bins_dim", check_schema=False), @@ -169,19 +170,27 @@ def create_l2_dataset( attrs=attrs.get_variable_attributes(key), ) + n_bins = histogram_l2.daily_lightcurve.number_of_bins for key, value in dataclasses.asdict(histogram_l2.daily_lightcurve).items(): if key == "number_of_bins": - # number_of_bins does not have n_bins dimensions. + # number_of_bins does not have a bins dimension. output[key] = xr.DataArray( np.array([value]), dims=["epoch"], attrs=attrs.get_variable_attributes(key), ) else: + # Bin arrays are chopped to number_of_bins in DailyLightcurve to + # avoid operating on FILLVAL data. Re-expand to STANDARD_BIN_COUNT + # here, filling unused bins with the variable's CDF FILLVAL. + var_attrs = attrs.get_variable_attributes(key) + fillval = var_attrs["FILLVAL"] + padded = np.full(GlowsConstants.STANDARD_BIN_COUNT, fillval) + padded[:n_bins] = value output[key] = xr.DataArray( - np.array([value]), + np.array([padded]), dims=["epoch", "bins"], - attrs=attrs.get_variable_attributes(key), + attrs=var_attrs, ) return output diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index 52a7510b3..4f87ed7b0 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -8,6 +8,7 @@ from imap_processing.glows import FLAG_LENGTH from imap_processing.glows.l1b.glows_l1b_data import PipelineSettings +from imap_processing.glows.utils.constants import GlowsConstants @dataclass @@ -65,9 +66,18 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None: L1B data filtered by good times, good angles, and good bins for one observation day. """ - self.raw_histograms = self.calculate_histogram_sums(l1b_data["histogram"].data) - - self.number_of_bins = l1b_data["histogram"].shape[1] + # number_of_bins_per_histogram is the count of valid (non-FILLVAL) bins. + # Histogram arrays from L1B are always GlowsConstants.STANDARD_BIN_COUNT + # (3600) long, with unused bins filled with GlowsConstants.HISTOGRAM_FILLVAL. + # All bin-dimensioned arrays here are chopped to number_of_bins so that + # computations only operate on valid data. glows_l2.py is responsible for + # re-expanding these arrays back to STANDARD_BIN_COUNT, filling unused bins + # with the appropriate CDF FILLVAL before writing to output. + self.number_of_bins = l1b_data["number_of_bins_per_histogram"].data[0] + + self.raw_histograms = self.calculate_histogram_sums( + l1b_data["histogram"].data + )[: self.number_of_bins] exposure_per_epoch = ( l1b_data["spin_period_average"].data @@ -79,16 +89,15 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None: self.exposure_times = np.full(self.number_of_bins, np.sum(exposure_per_epoch)) raw_uncertainties = np.sqrt(self.raw_histograms) - self.photon_flux = np.zeros(len(self.raw_histograms)) - self.flux_uncertainties = np.zeros(len(self.raw_histograms)) + self.photon_flux = np.zeros(self.number_of_bins) + self.flux_uncertainties = np.zeros(self.number_of_bins) # TODO: Only where exposure counts != 0 if len(self.exposure_times) != 0: self.photon_flux = self.raw_histograms / self.exposure_times self.flux_uncertainties = raw_uncertainties / self.exposure_times - # TODO: Average this, or should they all be the same? - self.spin_angle = np.average(l1b_data["imap_spin_angle_bin_cntr"].data, axis=0) + self.spin_angle = l1b_data["imap_spin_angle_bin_cntr"].data[0][: self.number_of_bins] self.histogram_flag_array = np.zeros(self.number_of_bins) self.ecliptic_lon = np.zeros(self.number_of_bins) @@ -110,7 +119,8 @@ def calculate_histogram_sums(histograms: NDArray) -> NDArray: Sum of valid histograms across all timestamps. """ histograms = histograms.copy() - histograms[histograms == -1] = 0 + # Zero out areas where HISTOGRAM_FILLVAL (i.e. unused bins) + histograms[histograms == GlowsConstants.HISTOGRAM_FILLVAL] = 0 return np.sum(histograms, axis=0, dtype=np.int64) From 5a7692bb6bb0e60b664cee9b064331346feada74 Mon Sep 17 00:00:00 2001 From: Maxine Hartnett Date: Wed, 18 Mar 2026 10:56:23 -0600 Subject: [PATCH 2/4] Correct GLOWS position angle average and std dev --- imap_processing/glows/l2/glows_l2_data.py | 53 +++++++++++++++++------ 1 file changed, 40 insertions(+), 13 deletions(-) diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index 6ffb5dcbf..19c242f01 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -9,6 +9,7 @@ from imap_processing.glows import FLAG_LENGTH from imap_processing.glows.l1b.glows_l1b_data import PipelineSettings from imap_processing.glows.utils.constants import GlowsConstants +from imap_processing.spice.geometry import SpiceFrame, get_instrument_mounting_az_el @dataclass @@ -48,7 +49,6 @@ class DailyLightcurve: raw_histograms: np.ndarray = field(init=False) exposure_times: np.ndarray = field(init=False) flux_uncertainties: np.ndarray = field(init=False) - # TODO: flag array histogram_flag_array: np.ndarray = field(init=False) # TODO: ecliptic coordinates ecliptic_lon: np.ndarray = field(init=False) @@ -73,7 +73,12 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None: # computations only operate on valid data. glows_l2.py is responsible for # re-expanding these arrays back to STANDARD_BIN_COUNT, filling unused bins # with the appropriate CDF FILLVAL before writing to output. - self.number_of_bins = l1b_data["number_of_bins_per_histogram"].data[0] + + self.number_of_bins = ( + l1b_data["number_of_bins_per_histogram"].data[0] + if len(l1b_data["number_of_bins_per_histogram"].data) != 0 + else 0 + ) self.raw_histograms = self.calculate_histogram_sums(l1b_data["histogram"].data)[ : self.number_of_bins @@ -97,10 +102,13 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None: self.photon_flux = self.raw_histograms / self.exposure_times self.flux_uncertainties = raw_uncertainties / self.exposure_times - self.spin_angle = l1b_data["imap_spin_angle_bin_cntr"].data[0][ - : self.number_of_bins - ] + if self.number_of_bins: + self.spin_angle = l1b_data["imap_spin_angle_bin_cntr"].data[0][ + : self.number_of_bins + ] + else: + self.spin_angle = np.zeros(0) # Apply 'OR' operation to histogram_flag_array across all # good-time L1B blocks per bin. # Per Section 12.3.4: a flag is True in L2 if it is True in any L1B block. @@ -311,15 +319,14 @@ def __init__(self, l1b_dataset: xr.Dataset, pipeline_settings: PipelineSettings) self.spin_period_ground_std_dev = ( good_data["spin_period_ground_average"].std(dim="epoch", keepdims=True).data ) - self.position_angle_offset_average = ( - good_data["position_angle_offset_average"] - .mean(dim="epoch", keepdims=True) - .data + + self.position_angle_offset_average = np.full( + (self.number_of_good_l1b_inputs), self.compute_position_angle() ) - self.position_angle_offset_std_dev = ( - good_data["position_angle_offset_average"] - .std(dim="epoch", keepdims=True) - .data + + # Always zero - per algorithm doc 10.6 + self.position_angle_offset_std_dev = np.full( + (self.number_of_good_l1b_inputs), 0.0 ) self.spacecraft_location_average = ( good_data["spacecraft_location_average"] @@ -404,3 +411,23 @@ def return_good_times(flags: xr.DataArray, active_flags: NDArray) -> NDArray: # where all the active indices == 1. good_times = np.where(np.all(flags[:, active_flags == 1] == 1, axis=1))[0] return good_times + + def compute_position_angle(self) -> float: + """ + Compute the position angle based on the instrument mounting. + + This number is not expected to change significantly. It is the same for all L1B + blocks (epoch values). + + Returns + ------- + float + The GLOWS mounting position angle. + """ + # Calculation described in algorithm doc 10.6 (Eq. 30): + # psi_G_eff = 360 - psi_GLOWS + # where psi_GLOWS is the azimuth of the GLOWS boresight in the + # IMAP spacecraft frame, measured from the spacecraft x-axis. + # delta_psi_G_eff is assumed to be 0 per instrument team decision. + glows_mounting_azimuth, _ = get_instrument_mounting_az_el(SpiceFrame.IMAP_GLOWS) + return (360.0 - glows_mounting_azimuth) % 360.0 From 5cccdb5019b17cfec0159b166ee41dc7a199c812 Mon Sep 17 00:00:00 2001 From: Maxine Hartnett Date: Wed, 18 Mar 2026 13:50:21 -0600 Subject: [PATCH 3/4] spin angle/positioning calculations --- imap_processing/glows/l2/glows_l2_data.py | 54 ++++--- imap_processing/tests/glows/test_glows_l2.py | 13 +- .../tests/glows/test_glows_l2_data.py | 137 ++++++++++++++++-- 3 files changed, 169 insertions(+), 35 deletions(-) diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index 19c242f01..2f4725c47 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -55,8 +55,9 @@ class DailyLightcurve: ecliptic_lat: np.ndarray = field(init=False) number_of_bins: int = field(init=False) l1b_data: InitVar[xr.Dataset] + position_angle: InitVar[float] - def __post_init__(self, l1b_data: xr.Dataset) -> None: + def __post_init__(self, l1b_data: xr.Dataset, position_angle: float) -> None: """ Compute all the daily lightcurve variables from L1B data. @@ -65,6 +66,11 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None: l1b_data : xarray.Dataset L1B data filtered by good times, good angles, and good bins for one observation day. + position_angle : float + The offset angle of the GLOWS instrument from the spin plane - this is used + in spin angle calculations. This number may change with different SPICE + kernels, but does not vary by time inside the code, so we can just use a + fixed value. """ # number_of_bins_per_histogram is the count of valid (non-FILLVAL) bins. # Histogram arrays from L1B are always GlowsConstants.STANDARD_BIN_COUNT @@ -102,13 +108,8 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None: self.photon_flux = self.raw_histograms / self.exposure_times self.flux_uncertainties = raw_uncertainties / self.exposure_times - if self.number_of_bins: - self.spin_angle = l1b_data["imap_spin_angle_bin_cntr"].data[0][ - : self.number_of_bins - ] + self.spin_angle = np.zeros(0) - else: - self.spin_angle = np.zeros(0) # Apply 'OR' operation to histogram_flag_array across all # good-time L1B blocks per bin. # Per Section 12.3.4: a flag is True in L2 if it is True in any L1B block. @@ -126,6 +127,28 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None: self.ecliptic_lon = np.zeros(self.number_of_bins) self.ecliptic_lat = np.zeros(self.number_of_bins) + if self.number_of_bins: + # imap_spin_angle_bin_cntr is the raw IMAP spin angle ψ (0 - 360°, + # bin midpoints). + spin_angle_bin_cntr = l1b_data["imap_spin_angle_bin_cntr"].data[0][ + : self.number_of_bins + ] + # Convert ψ → ψPA (Eq. 29): position angle measured from the + # northernmost point of the scanning circle. + self.spin_angle = (spin_angle_bin_cntr - position_angle + 360.0) % 360.0 + + # Roll all bin arrays so bin 0 corresponds to the northernmost + # point (minimum ψPA). + roll = -np.argsort(self.spin_angle)[0] + self.spin_angle = np.roll(self.spin_angle, roll) + self.raw_histograms = np.roll(self.raw_histograms, roll) + self.photon_flux = np.roll(self.photon_flux, roll) + self.exposure_times = np.roll(self.exposure_times, roll) + self.flux_uncertainties = np.roll(self.flux_uncertainties, roll) + self.histogram_flag_array = np.roll(self.histogram_flag_array, roll) + self.ecliptic_lon = np.roll(self.ecliptic_lon, roll) + self.ecliptic_lat = np.roll(self.ecliptic_lat, roll) + @staticmethod def calculate_histogram_sums(histograms: NDArray) -> NDArray: """ @@ -233,8 +256,8 @@ class HistogramL2: pulse_length_std_dev: np.ndarray[np.double] spin_period_ground_average: np.ndarray[np.double] spin_period_ground_std_dev: np.ndarray[np.double] - position_angle_offset_average: np.ndarray[np.double] - position_angle_offset_std_dev: np.ndarray[np.double] + position_angle_offset_average: np.double + position_angle_offset_std_dev: np.double spin_axis_orientation_std_dev: np.ndarray[np.double] spacecraft_location_average: np.ndarray[np.double] spacecraft_location_std_dev: np.ndarray[np.double] @@ -265,8 +288,6 @@ def __init__(self, l1b_dataset: xr.Dataset, pipeline_settings: PipelineSettings) # TODO: filter bad bins out. Needs to happen here while everything is still # per-timestamp. - self.daily_lightcurve = DailyLightcurve(good_data) - self.total_l1b_inputs = len(l1b_dataset["epoch"]) self.number_of_good_l1b_inputs = len(good_data["epoch"]) self.identifier = -1 # TODO: retrieve from spin table @@ -320,14 +341,11 @@ def __init__(self, l1b_dataset: xr.Dataset, pipeline_settings: PipelineSettings) good_data["spin_period_ground_average"].std(dim="epoch", keepdims=True).data ) - self.position_angle_offset_average = np.full( - (self.number_of_good_l1b_inputs), self.compute_position_angle() - ) + position_angle = self.compute_position_angle() + self.position_angle_offset_average: np.double = np.double(position_angle) # Always zero - per algorithm doc 10.6 - self.position_angle_offset_std_dev = np.full( - (self.number_of_good_l1b_inputs), 0.0 - ) + self.position_angle_offset_std_dev = np.double(0.0) self.spacecraft_location_average = ( good_data["spacecraft_location_average"] .mean(dim="epoch") @@ -359,6 +377,8 @@ def __init__(self, l1b_dataset: xr.Dataset, pipeline_settings: PipelineSettings) .data[np.newaxis, :] ) + self.daily_lightcurve = DailyLightcurve(good_data, position_angle) + def filter_bad_bins(self, histograms: NDArray, bin_exclusions: NDArray) -> NDArray: """ Filter out bad bins from the histogram. diff --git a/imap_processing/tests/glows/test_glows_l2.py b/imap_processing/tests/glows/test_glows_l2.py index a647397e4..ab7ab336b 100644 --- a/imap_processing/tests/glows/test_glows_l2.py +++ b/imap_processing/tests/glows/test_glows_l2.py @@ -13,6 +13,7 @@ glows_l2, ) from imap_processing.glows.l2.glows_l2_data import DailyLightcurve, HistogramL2 +from imap_processing.glows.utils.constants import GlowsConstants from imap_processing.spice.time import et_to_datetime64, ttj2000ns_to_et from imap_processing.tests.glows.conftest import mock_update_spice_parameters @@ -24,10 +25,10 @@ def l1b_hists(): hist = xr.DataArray( np.ones((4, 5)), dims=["epoch", "bins"], coords={"epoch": epoch, "bins": bins} ) - hist[1, 0] = -1 - hist[2, 0] = -1 - hist[1, 1] = -1 - hist[2, 3] = -1 + hist[1, 0] = GlowsConstants.HISTOGRAM_FILLVAL + hist[2, 0] = GlowsConstants.HISTOGRAM_FILLVAL + hist[1, 1] = GlowsConstants.HISTOGRAM_FILLVAL + hist[2, 3] = GlowsConstants.HISTOGRAM_FILLVAL input = xr.Dataset(coords={"epoch": epoch, "bins": bins}) input["histogram"] = hist @@ -35,6 +36,7 @@ def l1b_hists(): return input +@patch.object(HistogramL2, "compute_position_angle", return_value=42.0) @patch.object( HistogramL1B, "flag_uv_and_excluded", @@ -44,6 +46,7 @@ def l1b_hists(): def test_glows_l2( mock_spice_function, mock_flag_uv_and_excluded, + mock_compute_position_angle, l1a_dataset, mock_ancillary_exclusions, mock_pipeline_settings, @@ -75,6 +78,7 @@ def test_glows_l2( assert any(record.levelname == "WARNING" for record in caplog.records) +@patch.object(HistogramL2, "compute_position_angle", return_value=42.0) @patch.object( HistogramL1B, "flag_uv_and_excluded", @@ -84,6 +88,7 @@ def test_glows_l2( def test_generate_l2( mock_spice_function, mock_flag_uv_and_excluded, + mock_compute_position_angle, l1a_dataset, mock_ancillary_exclusions, mock_pipeline_settings, diff --git a/imap_processing/tests/glows/test_glows_l2_data.py b/imap_processing/tests/glows/test_glows_l2_data.py index 5b1cdfaf3..8885320a5 100644 --- a/imap_processing/tests/glows/test_glows_l2_data.py +++ b/imap_processing/tests/glows/test_glows_l2_data.py @@ -1,3 +1,5 @@ +from unittest.mock import patch + import numpy as np import pytest import xarray as xr @@ -42,8 +44,12 @@ def pipeline_settings(): ] pipeline_dataset = xr.Dataset( { - "active_bad_time_flags": xr.DataArray(active_bad_time_flags), - "active_bad_angle_flags": xr.DataArray(active_bad_angle_flags), + "active_bad_time_flags": xr.DataArray( + active_bad_time_flags, dims=["time_flag_index"] + ), + "active_bad_angle_flags": xr.DataArray( + active_bad_angle_flags, dims=["angle_flag_index"] + ), } ) return PipelineSettings(pipeline_dataset) @@ -62,9 +68,7 @@ def l1b_dataset(): epoch = xr.DataArray(np.arange(n_epochs), dims=["epoch"]) bins = xr.DataArray(np.arange(n_bins), dims=["bins"]) - histogram = np.array( - [[10, 20, 30, fillval], [10, 20, 30, 40]], dtype=float - ) + histogram = np.array([[10, 20, 30, fillval], [10, 20, 30, 40]], dtype=float) spin_angle = np.tile(np.linspace(0, 270, n_bins), (n_epochs, 1)) histogram_flag_array = np.zeros((n_epochs, n_flags, n_bins), dtype=np.uint8) @@ -78,7 +82,7 @@ def l1b_dataset(): ["epoch", "bad_angle_flags", "bins"], histogram_flag_array, ), - "number_of_bins_per_histogram": (["epoch"], [n_bins, n_bins]) + "number_of_bins_per_histogram": (["epoch"], [n_bins, n_bins]), }, coords={"epoch": epoch, "bins": bins}, ) @@ -87,7 +91,7 @@ def l1b_dataset(): def test_photon_flux(l1b_dataset): """Flux = sum(histograms) / sum(exposure_times) per bin (Eq. 50).""" - lc = DailyLightcurve(l1b_dataset) + lc = DailyLightcurve(l1b_dataset, position_angle=0.0) # l1b_exposure_time_per_bin = spin_period_average * # number_of_spins_per_block / number_of_bins_per_histogram @@ -106,7 +110,7 @@ def test_photon_flux(l1b_dataset): def test_flux_uncertainty(l1b_dataset): """Uncertainty = sqrt(sum_hist) / exposure per bin (Eq. 54).""" - lc = DailyLightcurve(l1b_dataset) + lc = DailyLightcurve(l1b_dataset, position_angle=0.0) expected_uncertainty = np.sqrt(lc.raw_histograms) / lc.exposure_times assert np.allclose(lc.flux_uncertainties, expected_uncertainty) @@ -120,7 +124,7 @@ def test_zero_exposure_bins(l1b_dataset): uncertainty are zero because the raw histogram sums are zero. """ l1b_dataset["histogram"].values[:] = GlowsConstants.HISTOGRAM_FILLVAL - lc = DailyLightcurve(l1b_dataset) + lc = DailyLightcurve(l1b_dataset, position_angle=0.0) expected_exposure = 2 * 15.0 * 5 / 4 assert np.all(lc.photon_flux == 0) @@ -129,7 +133,7 @@ def test_zero_exposure_bins(l1b_dataset): def test_number_of_bins(l1b_dataset): - lc = DailyLightcurve(l1b_dataset) + lc = DailyLightcurve(l1b_dataset, position_angle=0.0) assert lc.number_of_bins == 4 assert len(lc.spin_angle) == 4 assert len(lc.photon_flux) == 4 @@ -149,7 +153,7 @@ def test_histogram_flag_array_or_propagation(l1b_dataset): l1b_dataset["histogram_flag_array"].values[1, 1, 2] = 2 l1b_dataset["histogram_flag_array"].values[1, 0, 0] = 2 - lc = DailyLightcurve(l1b_dataset) + lc = DailyLightcurve(l1b_dataset, position_angle=0.0) assert ( lc.histogram_flag_array[0] == 3 @@ -163,7 +167,7 @@ def test_histogram_flag_array_zero_epochs(): """histogram_flag_array is all zeros when the input dataset is empty. Note: this is NEVER expected to happen in production - TODO: remove tests and allow failures if empty imput?""" + """ n_bins, n_flags = 4, 4 histogram = np.empty((0, n_bins), dtype=float) spin_angle = np.empty((0, n_bins), dtype=float) @@ -179,11 +183,11 @@ def test_histogram_flag_array_zero_epochs(): ["epoch", "bad_angle_flags", "bins"], histogram_flag_array, ), - "number_of_bins_per_histogram": (["epoch"], []) + "number_of_bins_per_histogram": (["epoch"], []), }, coords={"epoch": xr.DataArray(np.arange(0), dims=["epoch"])}, ) - lc = DailyLightcurve(ds) + lc = DailyLightcurve(ds, position_angle=0.0) # if the dataset is empty, there is no way to infer the number_of_bins assert len(lc.histogram_flag_array) == 0 @@ -203,3 +207,108 @@ def test_filter_good_times(): expected_good_times = [0, 2, 3] assert np.array_equal(good_times, expected_good_times) + + +# ── spin_angle tests ────────────────────────────────────────────────────────── + + +def test_spin_angle_offset_formula(l1b_dataset): + """spin_angle = (imap_spin_angle_bin_cntr - position_angle + 360) % 360. + + Fixture spin_angle_bin_cntr = [0, 90, 180, 270], position_angle = 90. + Expected before roll: [270, 0, 90, 180]. + Minimum is at index 1, so roll = -1 -> [0, 90, 180, 270]. + """ + lc = DailyLightcurve(l1b_dataset, position_angle=90.0) + expected = np.array([0.0, 90.0, 180.0, 270.0]) + assert np.allclose(lc.spin_angle, expected) + + +def test_spin_angle_starts_at_minimum(l1b_dataset): + """After rolling, lc.spin_angle[0] is the minimum value. + + Fixture spin_angle_bin_cntr = [0, 90, 180, 270], position_angle = 45. + Before roll: [315, 45, 135, 225]; minimum 45 is at index 1 -> roll = -1 + -> [45, 135, 225, 315]. + """ + lc = DailyLightcurve(l1b_dataset, position_angle=45.0) + assert lc.spin_angle[0] == np.min(lc.spin_angle) + assert np.allclose(lc.spin_angle, np.array([45.0, 135.0, 225.0, 315.0])) + + +# ── position_angle_offset_average tests ────────────────────────────────────── + + +def test_compute_position_angle(): + """compute_position_angle returns (360 - azimuth) % 360 (Eq. 30).""" + target_module = ( + "imap_processing.glows.l2.glows_l2_data.get_instrument_mounting_az_el" + ) + with patch(target_module, return_value=(270.0, 0.0)): + result = HistogramL2.compute_position_angle(None) + assert result == pytest.approx(90.0) + + +@pytest.fixture +def l1b_dataset_full(): + """Minimal L1B dataset with all variables required by HistogramL2. + + Two epochs, four bins, 17 flags (all good). + """ + n_epochs, n_bins, n_angle_flags, n_time_flags = 2, 4, 4, 17 + fillval = GlowsConstants.HISTOGRAM_FILLVAL + epoch = np.arange(n_epochs, dtype=np.float64) + + histogram = np.array([[10, 20, 30, fillval], [10, 20, 30, 40]], dtype=float) + spin_angle = np.tile(np.linspace(0, 270, n_bins), (n_epochs, 1)) + histogram_flag_array = np.zeros((n_epochs, n_angle_flags, n_bins), dtype=np.uint8) + + return xr.Dataset( + { + "histogram": (["epoch", "bins"], histogram), + "spin_period_average": (["epoch"], [15.0, 15.0]), + "number_of_spins_per_block": (["epoch"], [5, 5]), + "imap_spin_angle_bin_cntr": (["epoch", "bins"], spin_angle), + "histogram_flag_array": ( + ["epoch", "bad_angle_flags", "bins"], + histogram_flag_array, + ), + "number_of_bins_per_histogram": (["epoch"], [n_bins, n_bins]), + "flags": ( + ["epoch", "flag_index"], + np.ones((n_epochs, n_time_flags)), + ), + "filter_temperature_average": (["epoch"], [20.0, 21.0]), + "hv_voltage_average": (["epoch"], [1000.0, 1000.0]), + "pulse_length_average": (["epoch"], [5.0, 5.0]), + "spin_period_ground_average": (["epoch"], [15.0, 15.0]), + "spacecraft_location_average": ( + ["epoch", "xyz"], + np.ones((n_epochs, 3)), + ), + "spacecraft_velocity_average": ( + ["epoch", "xyz"], + np.ones((n_epochs, 3)), + ), + "spin_axis_orientation_average": ( + ["epoch", "lonlat"], + np.ones((n_epochs, 2)), + ), + "imap_start_time": (["epoch"], [0.0, 1.0]), + "imap_time_offset": (["epoch"], [60.0, 60.0]), + }, + coords={"epoch": xr.DataArray(epoch, dims=["epoch"])}, + ) + + +def test_position_angle_offset_average(l1b_dataset_full, pipeline_settings): + """position_angle_offset_average is a scalar equal to the result of + compute_position_angle (Eq. 30, Section 10.6). It is constant across the + observational day since it depends only on instrument mounting geometry. + """ + mock_pa = 42.5 + target = "imap_processing.glows.l2.glows_l2_data.HistogramL2.compute_position_angle" + with patch(target, return_value=mock_pa): + l2 = HistogramL2(l1b_dataset_full, pipeline_settings) + + assert l2.position_angle_offset_average == pytest.approx(mock_pa) From e8fb2cc4792bc1e645ac143339967143fbdfac8d Mon Sep 17 00:00:00 2001 From: Maxine Hartnett Date: Tue, 24 Mar 2026 13:31:17 -0600 Subject: [PATCH 4/4] PR updates --- imap_processing/glows/l2/glows_l2_data.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index 2f4725c47..d5860c112 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -67,10 +67,8 @@ def __post_init__(self, l1b_data: xr.Dataset, position_angle: float) -> None: L1B data filtered by good times, good angles, and good bins for one observation day. position_angle : float - The offset angle of the GLOWS instrument from the spin plane - this is used - in spin angle calculations. This number may change with different SPICE - kernels, but does not vary by time inside the code, so we can just use a - fixed value. + The offset angle of the GLOWS instrument from the north spin point - this + is used in spin angle calculations. """ # number_of_bins_per_histogram is the count of valid (non-FILLVAL) bins. # Histogram arrays from L1B are always GlowsConstants.STANDARD_BIN_COUNT @@ -139,7 +137,7 @@ def __post_init__(self, l1b_data: xr.Dataset, position_angle: float) -> None: # Roll all bin arrays so bin 0 corresponds to the northernmost # point (minimum ψPA). - roll = -np.argsort(self.spin_angle)[0] + roll = -np.argmin(self.spin_angle) self.spin_angle = np.roll(self.spin_angle, roll) self.raw_histograms = np.roll(self.raw_histograms, roll) self.photon_flux = np.roll(self.photon_flux, roll) @@ -448,6 +446,11 @@ def compute_position_angle(self) -> float: # psi_G_eff = 360 - psi_GLOWS # where psi_GLOWS is the azimuth of the GLOWS boresight in the # IMAP spacecraft frame, measured from the spacecraft x-axis. - # delta_psi_G_eff is assumed to be 0 per instrument team decision. + # This angle does not change with time, as it is in the spinning IMAP frame. + # It basically defines the angle between x=0 in the IMAP frame and x=0 in the + # GLOWS instrument frame, and is defined by the physical mounting location of + # the instrument. + # delta_psi_G_eff is assumed to be 0 per instrument team decision (aka this + # doesn't move from the SPICE determined mounting angle. glows_mounting_azimuth, _ = get_instrument_mounting_az_el(SpiceFrame.IMAP_GLOWS) return (360.0 - glows_mounting_azimuth) % 360.0