diff --git a/imap_processing/ancillary/ancillary_dataset_combiner.py b/imap_processing/ancillary/ancillary_dataset_combiner.py index 104255f87..8d97edb7e 100644 --- a/imap_processing/ancillary/ancillary_dataset_combiner.py +++ b/imap_processing/ancillary/ancillary_dataset_combiner.py @@ -431,13 +431,17 @@ def convert_file_to_dataset(self, filepath: str | Path) -> xr.Dataset: # noqa: lines = [line.strip() for line in f if not line.startswith("#")] identifiers = [line.split(" ", 1)[0] for line in lines] values = [float(line.split(" ", 1)[1]) for line in lines] - return xr.Dataset( + ds = xr.Dataset( { - "start_time_utc": (["time_block"], identifiers), - "cps_per_r": (["time_block"], values), - } + "cps_per_r": (["start_time_utc"], values), # floats + }, + coords={ + "start_time_utc": np.array(identifiers, dtype="datetime64[s]") + }, # (e.g. '2025-07-01T00:00:00') ) + return ds.sortby("start_time_utc") + elif filename.endswith(".json"): # Handle pipeline settings JSON file using the generic read_json method return self.convert_json_to_dataset(filepath) diff --git a/imap_processing/glows/l2/glows_l2.py b/imap_processing/glows/l2/glows_l2.py index f96ae0ef4..8e457e815 100644 --- a/imap_processing/glows/l2/glows_l2.py +++ b/imap_processing/glows/l2/glows_l2.py @@ -59,7 +59,7 @@ def glows_l2( pipeline_settings_dataset.sel(epoch=day, method="nearest") ) - l2 = HistogramL2(input_dataset, pipeline_settings) + l2 = HistogramL2(input_dataset, pipeline_settings, calibration_dataset) if l2.number_of_good_l1b_inputs == 0: logger.warning("No good data found in L1B dataset. Returning empty list.") return [] diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index 22c63950e..5949fd5ae 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -14,7 +14,12 @@ frame_transform_az_el, get_instrument_mounting_az_el, ) -from imap_processing.spice.time import met_to_sclkticks, sct_to_et +from imap_processing.spice.time import ( + et_to_datetime64, + met_to_sclkticks, + sct_to_et, + ttj2000ns_to_et, +) @dataclass @@ -46,6 +51,8 @@ class DailyLightcurve: number of bins in lightcurve l1b_data : xarray.Dataset L1B data filtered by good times, good angles, and good bins. + calibration_factor : float | None + Rayleigh calibration factor used for flux calculations. """ # All variables should have n_bin elements @@ -55,14 +62,19 @@ class DailyLightcurve: exposure_times: np.ndarray = field(init=False) flux_uncertainties: np.ndarray = field(init=False) histogram_flag_array: np.ndarray = field(init=False) - # TODO: ecliptic coordinates ecliptic_lon: np.ndarray = field(init=False) 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, position_angle: float) -> None: + calibration_factor: InitVar[float | None] + + def __post_init__( + self, + l1b_data: xr.Dataset, + position_angle: float, + calibration_factor: float | None, + ) -> None: """ Compute all the daily lightcurve variables from L1B data. @@ -74,6 +86,10 @@ def __post_init__(self, l1b_data: xr.Dataset, position_angle: float) -> None: position_angle : float The offset angle of the GLOWS instrument from the north spin point - this is used in spin angle calculations. + calibration_factor : float + Calibration factor used for flux calculations, in units of counts per second + per Rayleigh. This is used to convert from raw histograms and exposure times + to physical photon flux units. """ # number_of_bins_per_histogram is the count of valid (non-FILLVAL) bins. # Histogram arrays from L1B are always GlowsConstants.STANDARD_BIN_COUNT @@ -107,12 +123,16 @@ def __post_init__(self, l1b_data: xr.Dataset, position_angle: float) -> None: self.flux_uncertainties = np.zeros(self.number_of_bins) if ( - len(self.exposure_times) != 0 + self.number_of_bins > 0 and self.exposure_times[0] > 0 - and len(np.unique(self.exposure_times)) == 1 + and calibration_factor ): - self.photon_flux = self.raw_histograms / self.exposure_times - self.flux_uncertainties = raw_uncertainties / self.exposure_times + self.photon_flux = ( + self.raw_histograms / self.exposure_times + ) / calibration_factor + self.flux_uncertainties = ( + raw_uncertainties / self.exposure_times + ) / calibration_factor self.spin_angle = np.zeros(0) @@ -244,6 +264,8 @@ class HistogramL2: GLOWS histogram L1B dataset, as produced by glows_l1b.py. pipeline_settings : PipelineSettings Pipeline settings object read from ancillary file. + calibration_dataset : xr.Dataset + The cps-to-Rayleigh calibration dataset needed for flux calculations. Attributes ---------- @@ -327,7 +349,12 @@ class HistogramL2: spin_axis_orientation_average: np.ndarray[np.double] bad_time_flag_occurrences: np.ndarray - def __init__(self, l1b_dataset: xr.Dataset, pipeline_settings: PipelineSettings): + def __init__( + self, + l1b_dataset: xr.Dataset, + pipeline_settings: PipelineSettings, + calibration_dataset: xr.Dataset, + ) -> None: """ Given an L1B dataset, process data into an output HistogramL2 object. @@ -337,6 +364,9 @@ def __init__(self, l1b_dataset: xr.Dataset, pipeline_settings: PipelineSettings) GLOWS histogram L1B dataset, as produced by glows_l1b.py. pipeline_settings : PipelineSettings Pipeline settings object read from ancillary file. + calibration_dataset : xr.Dataset + cps-to-Rayleigh calibration dataset used for flux calculations. + coords: start_time_utc, data_vars: cps_per_r """ active_flags = np.array(pipeline_settings.active_bad_time_flags, dtype=float) @@ -438,7 +468,17 @@ def __init__(self, l1b_dataset: xr.Dataset, pipeline_settings: PipelineSettings) .data[np.newaxis, :] ) - self.daily_lightcurve = DailyLightcurve(good_data, position_angle) + # Select calibration factor corresponding to the mid-epoch in the L1B data. + if len(good_data["epoch"].data) != 0: + calibration_factor = self.get_calibration_factor( + good_data["epoch"].data, calibration_dataset + ) + else: + calibration_factor = None + + self.daily_lightcurve = DailyLightcurve( + good_data, position_angle, calibration_factor + ) def filter_bad_bins(self, histograms: NDArray, bin_exclusions: NDArray) -> NDArray: """ @@ -517,3 +557,33 @@ def compute_position_angle(self) -> float: # 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 + + @staticmethod + def get_calibration_factor( + epoch_values: np.ndarray, calibration_dataset: xr.Dataset + ) -> float: + """ + Select calibration factor for an observational day. + + The calibration factor is needed to compute flux in Rayleigh units. + There is a strong assumption that the calibration is constant for + a given observational day. + + Parameters + ---------- + epoch_values : np.ndarray + Array of epoch values from the L1B dataset, in TT J2000 nanoseconds. + calibration_dataset : xr.Dataset + Dataset containing calibration data. + + Returns + ------- + float + The calibration factor needed to compute flux in Rayleigh units. + """ + # Use the midpoint epoch for the day + mid_idx = len(epoch_values) // 2 + mid_epoch_utc = et_to_datetime64(ttj2000ns_to_et(epoch_values[mid_idx].item())) + return calibration_dataset.sel(start_time_utc=mid_epoch_utc, method="pad")[ + "cps_per_r" + ].data.item() diff --git a/imap_processing/tests/ancillary/test_ancillary_dataset_combiner.py b/imap_processing/tests/ancillary/test_ancillary_dataset_combiner.py index f7d3f0020..be784f0eb 100644 --- a/imap_processing/tests/ancillary/test_ancillary_dataset_combiner.py +++ b/imap_processing/tests/ancillary/test_ancillary_dataset_combiner.py @@ -316,15 +316,20 @@ def test_glows_exclusions_by_instr_team_combiner(glows_ancillary_filepath): def test_glows_l2_calibration_combiner(tmp_path): file_path = tmp_path / "imap_glows_l2-calibration_20251112_v001.dat" file_path.write_text( - "# header\n2025-11-13T18:12:48 1.020\n2025-11-14T09:58:04 0.849\n" + "# header\n2025-11-13T18:12:48 1.020\n" + "2025-11-14T09:58:04 0.849\n" + "2025-11-14T20:58:04 0.649\n" ) combiner = GlowsAncillaryCombiner([], "20251115") dataset = combiner.convert_file_to_dataset(file_path) - assert "start_time_utc" in dataset.data_vars + assert "start_time_utc" in dataset.coords + assert ( + np.diff(dataset.start_time_utc.values.astype("datetime64")) >= np.timedelta64(0) + ).all() assert "cps_per_r" in dataset.data_vars - assert len(dataset["cps_per_r"]) == 2 + assert len(dataset["cps_per_r"]) == 3 assert dataset["cps_per_r"].values[0] == pytest.approx(1.020) diff --git a/imap_processing/tests/glows/conftest.py b/imap_processing/tests/glows/conftest.py index 0b95cde28..8dc1118b5 100644 --- a/imap_processing/tests/glows/conftest.py +++ b/imap_processing/tests/glows/conftest.py @@ -236,6 +236,19 @@ def mock_conversion_table_dict(): return mock_dict +@pytest.fixture +def mock_calibration_dataset(): + """Create a mock CalibrationDataset object for testing.""" + return xr.Dataset( + {"cps_per_r": xr.DataArray([0.849, 1.020], dims=["start_time_utc"])}, + coords={ + "start_time_utc": np.array( + ["2011-09-19T09:58:04", "2011-09-20T18:12:48"], dtype="datetime64[s]" + ) + }, + ) + + @pytest.fixture def mock_pipeline_settings(): """Create a mock PipelineSettings dataset for testing.""" diff --git a/imap_processing/tests/glows/test_glows_l2.py b/imap_processing/tests/glows/test_glows_l2.py index faaec83ae..508f47c14 100644 --- a/imap_processing/tests/glows/test_glows_l2.py +++ b/imap_processing/tests/glows/test_glows_l2.py @@ -9,9 +9,7 @@ HistogramL1B, PipelineSettings, ) -from imap_processing.glows.l2.glows_l2 import ( - glows_l2, -) +from imap_processing.glows.l2.glows_l2 import 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 @@ -54,6 +52,7 @@ def test_glows_l2( mock_pipeline_settings, mock_conversion_table_dict, mock_ecliptic_bin_centers, + mock_calibration_dataset, caplog, ): mock_spice_function.side_effect = mock_update_spice_parameters @@ -69,7 +68,7 @@ def test_glows_l2( ) # Test case 1: L1B dataset has good times - l2 = glows_l2(l1b_hist_dataset, mock_pipeline_settings, None)[0] + l2 = glows_l2(l1b_hist_dataset, mock_pipeline_settings, mock_calibration_dataset)[0] assert l2.attrs["Logical_source"] == "imap_glows_l2_hist" assert np.allclose(l2["filter_temperature_average"].values, [57.6], rtol=0.1) @@ -79,7 +78,9 @@ def test_glows_l2( l1b_hist_dataset_no_good_times.flags.shape ) caplog.set_level("WARNING") - result = glows_l2(l1b_hist_dataset_no_good_times, mock_pipeline_settings, None) + result = glows_l2( + l1b_hist_dataset_no_good_times, mock_pipeline_settings, mock_calibration_dataset + ) assert result == [] assert any(record.levelname == "WARNING" for record in caplog.records) @@ -88,7 +89,9 @@ def test_glows_l2( l1b_hist_dataset_zero_values["spin_period_average"].data[:] = 0 l1b_hist_dataset_zero_values["number_of_spins_per_block"].data[:] = 0 caplog.set_level("WARNING") - result = glows_l2(l1b_hist_dataset_zero_values, mock_pipeline_settings, None) + result = glows_l2( + l1b_hist_dataset_zero_values, mock_pipeline_settings, mock_calibration_dataset + ) assert result == [] assert any(record.levelname == "WARNING" for record in caplog.records) @@ -109,6 +112,7 @@ def test_generate_l2( mock_pipeline_settings, mock_conversion_table_dict, mock_ecliptic_bin_centers, + mock_calibration_dataset, ): mock_spice_function.side_effect = mock_update_spice_parameters @@ -127,37 +131,38 @@ def test_generate_l2( ) # Test case 1: L1B dataset has good times - l2 = HistogramL2(l1b_hist_dataset, pipeline_settings) - - expected_values = { - "filter_temperature_average": [57.59], - "filter_temperature_std_dev": [0.21], - "hv_voltage_average": [1715.4], - "hv_voltage_std_dev": [0.0], - } - - assert np.isclose( - l2.filter_temperature_average, - expected_values["filter_temperature_average"], - 0.01, - ) - assert np.isclose( - l2.filter_temperature_std_dev, - expected_values["filter_temperature_std_dev"], - 0.01, - ) - assert np.isclose( - l2.hv_voltage_average, expected_values["hv_voltage_average"], 0.01 - ) - assert np.isclose( - l2.hv_voltage_std_dev, expected_values["hv_voltage_std_dev"], 0.01 - ) - - # Test case 2: L1B dataset has no good times (all flags 0) - l1b_hist_dataset["flags"].values = np.zeros(l1b_hist_dataset.flags.shape) - ds = HistogramL2(l1b_hist_dataset, pipeline_settings) - expected_number_of_good_l1b_inputs = 0 - assert ds.number_of_good_l1b_inputs == expected_number_of_good_l1b_inputs + with patch.object(HistogramL2, "get_calibration_factor", return_value=1): + l2 = HistogramL2(l1b_hist_dataset, pipeline_settings, mock_calibration_dataset) + + expected_values = { + "filter_temperature_average": [57.59], + "filter_temperature_std_dev": [0.21], + "hv_voltage_average": [1715.4], + "hv_voltage_std_dev": [0.0], + } + + assert np.isclose( + l2.filter_temperature_average, + expected_values["filter_temperature_average"], + 0.01, + ) + assert np.isclose( + l2.filter_temperature_std_dev, + expected_values["filter_temperature_std_dev"], + 0.01, + ) + assert np.isclose( + l2.hv_voltage_average, expected_values["hv_voltage_average"], 0.01 + ) + assert np.isclose( + l2.hv_voltage_std_dev, expected_values["hv_voltage_std_dev"], 0.01 + ) + + # Test case 2: L1B dataset has no good times (all flags 0) + l1b_hist_dataset["flags"].values = np.zeros(l1b_hist_dataset.flags.shape) + ds = HistogramL2(l1b_hist_dataset, pipeline_settings, mock_calibration_dataset) + expected_number_of_good_l1b_inputs = 0 + assert ds.number_of_good_l1b_inputs == expected_number_of_good_l1b_inputs def test_bin_exclusions(l1b_hists): diff --git a/imap_processing/tests/glows/test_glows_l2_data.py b/imap_processing/tests/glows/test_glows_l2_data.py index 5ffe6f2fa..40ae857f2 100644 --- a/imap_processing/tests/glows/test_glows_l2_data.py +++ b/imap_processing/tests/glows/test_glows_l2_data.py @@ -91,6 +91,32 @@ def l1b_dataset(): return ds +def test_get_calibration_factor(mock_calibration_dataset): + """Test selecting correct calibration factor.""" + + # The mid-epoch is after calibration timestamps, + # so the most recent (1.020) is selected. + # ['2011-09-21T00:50:15.000', '2011-09-21T00:52:15.000', '2011-09-21T00:54:15.000'] + later_epoch = np.array([369838281184000000, 369838401184000000, 369838521184000000]) + assert HistogramL2.get_calibration_factor( + later_epoch, mock_calibration_dataset + ) == pytest.approx(1.020) + + # The mid-epoch is before all calibration timestamps, + # so a KeyError is raised with the "pad" filter method. + # ['2011-09-18T19:59:08.816', '2011-09-18T20:01:08.816', '2011-09-18T20:03:08.816'] + early_epoch = np.array([369648015000000000, 369648135000000000, 369648255000000000]) + with pytest.raises(KeyError): + HistogramL2.get_calibration_factor(early_epoch, mock_calibration_dataset) + + # The mid-epoch is between the calibration times, + # so the first entry (0.849) is selected. + between_epoch = np.array([369808281184000000]) + assert HistogramL2.get_calibration_factor( + between_epoch, mock_calibration_dataset + ) == pytest.approx(0.849) + + @pytest.mark.external_kernel def test_ecliptic_coords_computation(furnish_kernels): """Test method that computes ecliptic coordinates.""" @@ -133,8 +159,16 @@ def test_ecliptic_coords_computation(furnish_kernels): 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) + """ + Flux = (sum(histograms) / sum(exposure_times)) / + Rayleigh calibration factor + + per bin (Eq. 50-53) + """ + mock_cal_factor = 2 + lc = DailyLightcurve( + l1b_dataset, position_angle=0.0, calibration_factor=mock_cal_factor + ) # l1b_exposure_time_per_bin = spin_period_average * # number_of_spins_per_block / number_of_bins_per_histogram @@ -144,7 +178,7 @@ def test_photon_flux(l1b_dataset, mock_ecliptic_bin_centers): expected_exposure = np.array( [2 * exposure_per, 2 * exposure_per, 2 * exposure_per, 2 * exposure_per] ) - expected_flux = expected_raw / expected_exposure + expected_flux = (expected_raw / expected_exposure) / mock_cal_factor assert np.allclose(lc.raw_histograms, expected_raw) assert np.allclose(lc.exposure_times, expected_exposure) @@ -152,10 +186,19 @@ def test_photon_flux(l1b_dataset, mock_ecliptic_bin_centers): 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) + """ + Uncertainty = sqrt(sum_hist) / exposure / + Rayleigh calibration factor + + per bin (Eq. 54-55).""" + mock_cal_factor = 2 + lc = DailyLightcurve( + l1b_dataset, position_angle=0.0, calibration_factor=mock_cal_factor + ) - expected_uncertainty = np.sqrt(lc.raw_histograms) / lc.exposure_times + expected_uncertainty = ( + np.sqrt(lc.raw_histograms) / lc.exposure_times + ) / mock_cal_factor assert np.allclose(lc.flux_uncertainties, expected_uncertainty) @@ -166,8 +209,11 @@ def test_zero_exposure_bins(l1b_dataset, mock_ecliptic_bin_centers): when all histogram values are masked (HISTOGRAM_FILLVAL). Flux and uncertainty are zero because the raw histogram sums are zero. """ + mock_cal_factor = 1 l1b_dataset["histogram"].values[:] = GlowsConstants.HISTOGRAM_FILLVAL - lc = DailyLightcurve(l1b_dataset, position_angle=0.0) + lc = DailyLightcurve( + l1b_dataset, position_angle=0.0, calibration_factor=mock_cal_factor + ) expected_exposure = 2 * 15.0 * 5 / 4 assert np.all(lc.photon_flux == 0) @@ -185,8 +231,12 @@ def test_zero_exposure_values(l1b_dataset, mock_ecliptic_bin_centers): l1b_dataset["spin_period_average"].data[:] = 0 l1b_dataset["number_of_spins_per_block"].data[:] = 0 + mock_cal_factor = 1 + with np.errstate(divide="raise", invalid="raise"): - lc = DailyLightcurve(l1b_dataset, position_angle=0.0) + lc = DailyLightcurve( + l1b_dataset, position_angle=0.0, calibration_factor=mock_cal_factor + ) expected = np.zeros(l1b_dataset.sizes["bins"], dtype=float) assert lc.exposure_times.shape == expected.shape @@ -199,7 +249,10 @@ def test_zero_exposure_values(l1b_dataset, mock_ecliptic_bin_centers): def test_number_of_bins(l1b_dataset, mock_ecliptic_bin_centers): - lc = DailyLightcurve(l1b_dataset, position_angle=0.0) + mock_cal_factor = 1 + lc = DailyLightcurve( + l1b_dataset, position_angle=0.0, calibration_factor=mock_cal_factor + ) assert lc.number_of_bins == 4 assert len(lc.spin_angle) == 4 assert len(lc.photon_flux) == 4 @@ -221,7 +274,10 @@ def test_histogram_flag_array_or_propagation(l1b_dataset, mock_ecliptic_bin_cent l1b_dataset["histogram_flag_array"].values[1, 1, 2] = 2 l1b_dataset["histogram_flag_array"].values[1, 0, 0] = 2 - lc = DailyLightcurve(l1b_dataset, position_angle=0.0) + mock_cal_factor = 1 + lc = DailyLightcurve( + l1b_dataset, position_angle=0.0, calibration_factor=mock_cal_factor + ) assert ( lc.histogram_flag_array[0] == 3 @@ -255,7 +311,8 @@ def test_histogram_flag_array_zero_epochs(mock_ecliptic_bin_centers): }, coords={"epoch": xr.DataArray(np.arange(0), dims=["epoch"])}, ) - lc = DailyLightcurve(ds, position_angle=0.0) + mock_cal_factor = 1 + lc = DailyLightcurve(ds, position_angle=0.0, calibration_factor=mock_cal_factor) # if the dataset is empty, there is no way to infer the number_of_bins assert len(lc.histogram_flag_array) == 0 @@ -287,7 +344,10 @@ def test_spin_angle_offset_formula(l1b_dataset, mock_ecliptic_bin_centers): 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) + mock_cal_factor = 1 + lc = DailyLightcurve( + l1b_dataset, position_angle=90.0, calibration_factor=mock_cal_factor + ) expected = np.array([0.0, 90.0, 180.0, 270.0]) assert np.allclose(lc.spin_angle, expected) @@ -299,7 +359,10 @@ def test_spin_angle_starts_at_minimum(l1b_dataset, mock_ecliptic_bin_centers): 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) + mock_cal_factor = 1 + lc = DailyLightcurve( + l1b_dataset, position_angle=45.0, calibration_factor=mock_cal_factor + ) 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])) @@ -370,15 +433,24 @@ def l1b_dataset_full(): def test_position_angle_offset_average( - l1b_dataset_full, pipeline_settings, mock_ecliptic_bin_centers + l1b_dataset_full, + pipeline_settings, + mock_ecliptic_bin_centers, + mock_calibration_dataset, ): """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) + mock_cal_factor = 1 + + with ( + patch.object(HistogramL2, "compute_position_angle", return_value=mock_pa), + patch.object( + HistogramL2, "get_calibration_factor", return_value=mock_cal_factor + ), + ): + l2 = HistogramL2(l1b_dataset_full, pipeline_settings, mock_calibration_dataset) - assert l2.position_angle_offset_average == pytest.approx(mock_pa) + assert l2.position_angle_offset_average == pytest.approx(mock_pa)