diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index d5860c112..c9fb7cf6d 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -9,7 +9,12 @@ 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 +from imap_processing.spice.geometry import ( + SpiceFrame, + frame_transform_az_el, + get_instrument_mounting_az_el, +) +from imap_processing.spice.time import met_to_sclkticks, sct_to_et @dataclass @@ -122,8 +127,6 @@ def __post_init__(self, l1b_data: xr.Dataset, position_angle: float) -> None: ) else: self.histogram_flag_array = np.zeros(self.number_of_bins, dtype=np.uint8) - 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°, @@ -144,8 +147,18 @@ def __post_init__(self, l1b_data: xr.Dataset, position_angle: float) -> None: 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) + + # Get the midpoint start time covered by repointing kernels + # needed to compute ecliptic coordinates + mid_idx = len(l1b_data["imap_start_time"]) // 2 + pointing_midpoint_time_et = sct_to_et( + met_to_sclkticks(l1b_data["imap_start_time"][mid_idx].data) + ) + self.ecliptic_lon, self.ecliptic_lat = ( + self.compute_ecliptic_coords_of_bin_centers( + pointing_midpoint_time_et, self.spin_angle + ) + ) @staticmethod def calculate_histogram_sums(histograms: NDArray) -> NDArray: @@ -167,6 +180,53 @@ def calculate_histogram_sums(histograms: NDArray) -> NDArray: histograms[histograms == GlowsConstants.HISTOGRAM_FILLVAL] = 0 return np.sum(histograms, axis=0, dtype=np.int64) + @staticmethod + def compute_ecliptic_coords_of_bin_centers( + data_time_et: float, spin_angle_bin_centers: NDArray + ) -> tuple[np.ndarray, np.ndarray]: + """ + Compute the ecliptic coordinates of the histogram bin centers. + + This method transforms the instrument pointing direction for each bin + center from the IMAP Pointing frame (IMAP_DPS) to the ECLIPJ2000 frame. + + Parameters + ---------- + data_time_et : float + Ephemeris time corresponding to the midpoint of the histogram accumulation. + + spin_angle_bin_centers : numpy.ndarray + Spin angle bin centers for the histogram bins, measured in the IMAP frame, + with shape (n_bins,), and already corrected for the northernmost point of + the scanning circle. + + Returns + ------- + tuple[numpy.ndarray, numpy.ndarray] + Longitude and latitudes in the ECLIPJ2000 frame representing the pointing + direction of each histogram bin center, with shape (n_bins,). + """ + # In the IMAP frame, the azimuth corresponds to the spin angle bin centers + azimuth = spin_angle_bin_centers + + # Get elevation from instrument pointing direction in the DPS frame. + az_el_instrument_mounting = get_instrument_mounting_az_el(SpiceFrame.IMAP_GLOWS) + elevation = az_el_instrument_mounting[1] + + # Create array of azimuth, elevation coordinates in the DPS frame (n_bins, 2) + az_el = np.stack((azimuth, np.full_like(azimuth, elevation)), axis=-1) + + # Transform coordinates to ECLIPJ2000 frame using SPICE transformations. + ecliptic_coords = frame_transform_az_el( + data_time_et, + az_el, + SpiceFrame.IMAP_DPS, + SpiceFrame.ECLIPJ2000, + ) + + # Return ecliptic longitudes and latitudes as separate arrays. + return ecliptic_coords[:, 0], ecliptic_coords[:, 1] + @dataclass class HistogramL2: diff --git a/imap_processing/tests/glows/conftest.py b/imap_processing/tests/glows/conftest.py index c2d573f73..0b95cde28 100644 --- a/imap_processing/tests/glows/conftest.py +++ b/imap_processing/tests/glows/conftest.py @@ -13,6 +13,7 @@ AncillaryParameters, ) from imap_processing.glows.l2.glows_l2 import glows_l2 +from imap_processing.glows.l2.glows_l2_data import DailyLightcurve @pytest.fixture @@ -278,6 +279,23 @@ def mock_pipeline_settings(): return mock_pipeline_dataset +@pytest.fixture +def mock_ecliptic_bin_centers(monkeypatch): + """Mock ecliptic coordinates for bin centers.""" + + def _mock_compute_coords( + _data_start_time_et: float, spin_angle: np.ndarray + ) -> tuple[np.ndarray, np.ndarray]: + n_bins = len(spin_angle) + return np.zeros(n_bins, dtype=float), np.zeros(n_bins, dtype=float) + + monkeypatch.setattr( + DailyLightcurve, + "compute_ecliptic_coords_of_bin_centers", + staticmethod(_mock_compute_coords), + ) + + def mock_update_spice_parameters(self, *args, **kwargs): self.spin_period_ground_average = np.float64(0.0) self.spin_period_ground_std_dev = np.float64(0.0) diff --git a/imap_processing/tests/glows/test_glows_l2.py b/imap_processing/tests/glows/test_glows_l2.py index ab7ab336b..d7632f79d 100644 --- a/imap_processing/tests/glows/test_glows_l2.py +++ b/imap_processing/tests/glows/test_glows_l2.py @@ -15,7 +15,9 @@ 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 +from imap_processing.tests.glows.conftest import ( + mock_update_spice_parameters, +) @pytest.fixture @@ -51,6 +53,7 @@ def test_glows_l2( mock_ancillary_exclusions, mock_pipeline_settings, mock_conversion_table_dict, + mock_ecliptic_bin_centers, caplog, ): mock_spice_function.side_effect = mock_update_spice_parameters @@ -93,6 +96,7 @@ def test_generate_l2( mock_ancillary_exclusions, mock_pipeline_settings, mock_conversion_table_dict, + mock_ecliptic_bin_centers, ): mock_spice_function.side_effect = mock_update_spice_parameters diff --git a/imap_processing/tests/glows/test_glows_l2_data.py b/imap_processing/tests/glows/test_glows_l2_data.py index 8885320a5..d343fd630 100644 --- a/imap_processing/tests/glows/test_glows_l2_data.py +++ b/imap_processing/tests/glows/test_glows_l2_data.py @@ -7,6 +7,7 @@ from imap_processing.glows.l1b.glows_l1b_data import PipelineSettings 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 met_to_sclkticks, sct_to_et @pytest.fixture @@ -78,6 +79,7 @@ def l1b_dataset(): "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), + "imap_start_time": (["epoch"], [0.0, 1.0]), "histogram_flag_array": ( ["epoch", "bad_angle_flags", "bins"], histogram_flag_array, @@ -89,7 +91,48 @@ def l1b_dataset(): return ds -def test_photon_flux(l1b_dataset): +@pytest.mark.external_kernel +def test_ecliptic_coords_computation(furnish_kernels): + """Test method that computes ecliptic coordinates.""" + + # Use a met value within the SPICE kernel coverage (2026-01-01). + data_start_time_et = sct_to_et(met_to_sclkticks(504975603.125)) + n_bins = 4 + spin_angle = np.linspace(0, 270, n_bins) + + kernels = [ + "naif0012.tls", + "imap_sclk_0000.tsc", + "imap_130.tf", + "imap_science_120.tf", + "sim_1yr_imap_pointing_frame.bc", + ] + + with furnish_kernels(kernels): + ecliptic_lon, ecliptic_lat = ( + DailyLightcurve.compute_ecliptic_coords_of_bin_centers( + data_start_time_et, spin_angle + ) + ) + + # ecliptic_lon and ecliptic_lat must have one entry per bin + assert len(ecliptic_lon) == n_bins + assert len(ecliptic_lat) == n_bins + + # ecliptic longitude must be in [0, 360) + assert np.all(ecliptic_lon >= 0.0) + assert np.all(ecliptic_lon < 360.0) + + # ecliptic latitude must be in [-90, 90] + assert np.all(ecliptic_lat >= -90.0) + assert np.all(ecliptic_lat <= 90.0) + + # values must be finite (no NaN / Inf from SPICE) + assert np.all(np.isfinite(ecliptic_lon)) + assert np.all(np.isfinite(ecliptic_lat)) + + +def test_photon_flux(l1b_dataset, mock_ecliptic_bin_centers): """Flux = sum(histograms) / sum(exposure_times) per bin (Eq. 50).""" lc = DailyLightcurve(l1b_dataset, position_angle=0.0) @@ -108,7 +151,7 @@ def test_photon_flux(l1b_dataset): assert np.allclose(lc.photon_flux, expected_flux) -def test_flux_uncertainty(l1b_dataset): +def test_flux_uncertainty(l1b_dataset, mock_ecliptic_bin_centers): """Uncertainty = sqrt(sum_hist) / exposure per bin (Eq. 54).""" lc = DailyLightcurve(l1b_dataset, position_angle=0.0) @@ -116,7 +159,7 @@ def test_flux_uncertainty(l1b_dataset): assert np.allclose(lc.flux_uncertainties, expected_uncertainty) -def test_zero_exposure_bins(l1b_dataset): +def test_zero_exposure_bins(l1b_dataset, mock_ecliptic_bin_centers): """Bins with all-masked histograms get zero flux and uncertainty. Exposure time still accumulates uniformly from each good-time file even @@ -132,16 +175,18 @@ def test_zero_exposure_bins(l1b_dataset): assert np.allclose(lc.exposure_times, expected_exposure) -def test_number_of_bins(l1b_dataset): +def test_number_of_bins(l1b_dataset, mock_ecliptic_bin_centers): 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 assert len(lc.flux_uncertainties) == 4 assert len(lc.exposure_times) == 4 + assert len(lc.ecliptic_lon) == 4 + assert len(lc.ecliptic_lat) == 4 -def test_histogram_flag_array_or_propagation(l1b_dataset): +def test_histogram_flag_array_or_propagation(l1b_dataset, mock_ecliptic_bin_centers): """histogram_flag_array is OR'd across all L1B epochs and flag rows per bin. Per Section 12.3.4: a flag is True in L2 if it is True in any L1B block. @@ -163,7 +208,7 @@ def test_histogram_flag_array_or_propagation(l1b_dataset): assert lc.histogram_flag_array[3] == 0 # no flags on bin 3 -def test_histogram_flag_array_zero_epochs(): +def test_histogram_flag_array_zero_epochs(mock_ecliptic_bin_centers): """histogram_flag_array is all zeros when the input dataset is empty. Note: this is NEVER expected to happen in production @@ -212,7 +257,7 @@ def test_filter_good_times(): # ── spin_angle tests ────────────────────────────────────────────────────────── -def test_spin_angle_offset_formula(l1b_dataset): +def test_spin_angle_offset_formula(l1b_dataset, mock_ecliptic_bin_centers): """spin_angle = (imap_spin_angle_bin_cntr - position_angle + 360) % 360. Fixture spin_angle_bin_cntr = [0, 90, 180, 270], position_angle = 90. @@ -224,7 +269,7 @@ def test_spin_angle_offset_formula(l1b_dataset): assert np.allclose(lc.spin_angle, expected) -def test_spin_angle_starts_at_minimum(l1b_dataset): +def test_spin_angle_starts_at_minimum(l1b_dataset, mock_ecliptic_bin_centers): """After rolling, lc.spin_angle[0] is the minimum value. Fixture spin_angle_bin_cntr = [0, 90, 180, 270], position_angle = 45. @@ -301,7 +346,9 @@ def l1b_dataset_full(): ) -def test_position_angle_offset_average(l1b_dataset_full, pipeline_settings): +def test_position_angle_offset_average( + l1b_dataset_full, pipeline_settings, mock_ecliptic_bin_centers +): """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.