diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 020c11b0f2..d71ab303d3 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -28,7 +28,11 @@ get_pointing_times, interpolate_repoint_data, ) -from imap_processing.spice.spin import get_spin_data, get_spin_number +from imap_processing.spice.spin import ( + get_spin_data, + get_spin_number, + interpolate_spin_data, +) from imap_processing.spice.time import ( epoch_to_fractional_doy, et_to_utc, @@ -166,6 +170,7 @@ "spin_cycle", ] # ------------------------------------------------------------------- +DE_CLOCK_TICK_S = 4.096e-3 # seconds per DE clock tick def lo_l1b( @@ -277,13 +282,9 @@ def l1b_de( avg_spin_durations_per_cycle = get_avg_spin_durations_per_cycle(spin_data) # set the spin cycle for each direct event l1b_de = set_spin_cycle(pointing_start_met, l1a_de, l1b_de) - # get spin start times for each event - spin_start_time = get_spin_start_times(l1a_de, l1b_de, spin_data) # get the absolute met for each event - l1b_de = set_event_met( - l1a_de, l1b_de, spin_start_time, avg_spin_durations_per_cycle - ) + l1b_de = set_event_met(l1a_de, l1b_de) # set the epoch for each event l1b_de = set_each_event_epoch(l1b_de) # Set the ESA mode for each direct event @@ -794,7 +795,7 @@ def _check_sufficient_spins(spin_data: xr.Dataset) -> np.ndarray: def get_spin_start_times( - l1a_de: xr.Dataset, l1b_de: xr.Dataset, spin_data: xr.Dataset + l1a_de: xr.Dataset, ) -> xr.DataArray: """ Get the start time for the spin that each direct event is in. @@ -807,59 +808,31 @@ def get_spin_start_times( ---------- l1a_de : xr.Dataset The L1A DE dataset. - l1b_de : xr.Dataset - The L1B DE dataset. - spin_data : xr.Dataset - The L1A Spin dataset. Returns ------- - spin_start_time : xr.DataArray + spin_start_time : np.ndarray The start time for the spin that each direct event is in. """ - # align l1a_de packets with spin_data packets - de_to_spin_indices = match_science_to_spin_asc( - l1a_de["epoch"].values, spin_data["epoch"].values - ) - # Repeat this for each direct event based on the de_count - de_to_spin_indices = np.repeat(de_to_spin_indices, l1a_de["de_count"]) - - # There are 28 spins per epoch (1 aggregated science cycle) - # Set the spin_cycle_num to the spin number relative to the - # start of the ASC - spin_cycle_num = l1b_de["spin_cycle"] % 28 - asc_starts = spin_data["acq_start_sec"] + spin_data["acq_start_subsec"] * 1e-6 - avg_spin_durations = get_avg_spin_durations_per_cycle(spin_data) + # Get the actual spin start times from the spin data + # Use the individual spin start times rather than calculating from ASC averages + spin_start_times = interpolate_spin_data(l1a_de["shcoarse"].values)[ + "spin_start_met" + ].values - # Calculate the time based off of the start of the acquisition period - # then using an average spin duration to calculate the offset within the ASC - # NOTE: We don't want to use the spin start times directly from the ASC spin packet - # because we are using an average spin_cycle for the ASC and there are - # times when only half the spins were completed in an ASC. This allows us - # to still get a valid spin_cycle start time for each direct event, even - # if the average spin_cycle was after the end of the acquisition period. - # TODO: Can we do even better by knowing how many esa_steps and spins were complete? - # i.e. change the spin_cycle calculation - spin_start_time = ( - asc_starts[de_to_spin_indices] - + spin_cycle_num * avg_spin_durations[de_to_spin_indices] - ) - return spin_start_time + return spin_start_times def set_event_met( l1a_de: xr.Dataset, l1b_de: xr.Dataset, - spin_start_time: xr.DataArray, - avg_spin_durations: xr.DataArray, ) -> xr.Dataset: """ Get the event MET for each direct event. Each direct event is converted from a data number to engineering unit in seconds. - de_eu_time de_dn_time / 4096 * avg_spin_duration - where de_time is the direct event time Data Number (DN) and avg_spin_duration - is the average spin duration for the ASC that the event was measured in. + time_from_start_of_spin = de_time * DE_CLOCK_TICK_S + where de_time is the direct event time Data Number (DN). The direct event time is the time of direct event relative to the start of the spin. The event MET is the sum of the start time of the spin and the @@ -871,26 +844,18 @@ def set_event_met( The L1A DE dataset. l1b_de : xr.Dataset The L1B DE dataset. - spin_start_time : np.ndarray - The start time for the spin that each direct event is in. - avg_spin_durations : xr.DataArray - The average spin duration for each epoch. Returns ------- l1b_de : xr.Dataset The L1B DE dataset with the event MET. """ - counts = l1a_de["de_count"].values - de_time_asc_groups = np.split(l1a_de["de_time"].values, np.cumsum(counts)[:-1]) - de_times_eu = [] - for i, de_time_asc in enumerate(de_time_asc_groups): - # DE Time is 12 bit DN. The max possible value is 4095 - # divide by 4096 to get fraction of a spin duration - de_times_eu.extend(de_time_asc / 4096 * avg_spin_durations[i].values) + # get spin start times for each event + spin_start_times = get_spin_start_times(l1a_de) + # spin start + offset based on de_time ticks l1b_de["event_met"] = xr.DataArray( - spin_start_time + de_times_eu, + spin_start_times + l1a_de["de_time"].values * DE_CLOCK_TICK_S, dims=["epoch"], # attrs=attr_mgr.get_variable_attributes("epoch") ) @@ -1344,12 +1309,14 @@ def set_pointing_bin(l1b_de: xr.Dataset) -> xr.Dataset: # first column: radius (Not needed) # second column: longitude lons = direction[:, 1] + # shift to 0-360 range (spin-phase 0 should be in bin 0) + lons = (lons + 360) % 360 # third column: latitude lats = direction[:, 2] # Define bin edges # 3600 bins, 0.1° each - lon_bins = np.linspace(-180, 180, 3601) + lon_bins = np.linspace(0, 360, 3601) # 40 bins, 0.1° each lat_bins = np.linspace(-2, 2, 41) diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 9f69f3cc03..d1132196bc 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -11,6 +11,7 @@ from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.cdf.utils import load_cdf from imap_processing.lo.l1b.lo_l1b import ( + DE_CLOCK_TICK_S, calculate_de_rates, calculate_histogram_rates, calculate_star_sensor_profile_for_group, @@ -219,7 +220,9 @@ def l1a_hist(): "imap_processing.lo.l1b.lo_l1b.cartesian_to_latitudinal", return_value=np.zeros((2000, 3)), ) +@patch("imap_processing.lo.l1b.lo_l1b.interpolate_spin_data") def test_lo_l1b_de( + mock_interpolate_spin_data, mock_frame_transform, mock_lo_instrument_pointing, mocked_get_pointing_times, @@ -229,6 +232,15 @@ def test_lo_l1b_de( anc_dependencies, ): # Arrange + # Mock the spin data to provide spin start times + # Create a DataFrame covering the time range of the test data + num_events = 2000 + mock_spin_df = pd.DataFrame( + { + "spin_start_met": np.ones(num_events), + } + ) + mock_interpolate_spin_data.return_value = mock_spin_df # Add l1b_nhk dependency with pivot angle information l1b_nhk = xr.Dataset( @@ -490,62 +502,37 @@ def test_spin_cycle(mock_get_spin_number): np.testing.assert_array_equal(spin_cycle_data["spin_cycle"], spin_cycle_expected) -def test_get_spin_start_times(): +@patch("imap_processing.lo.l1b.lo_l1b.interpolate_spin_data") +def test_get_spin_start_times(mock_interpolate_spin_data): # Arrange - l1b_de = xr.Dataset( + # Mock the spin data to return specific spin start times + mock_spin_df = pd.DataFrame( { - "spin_cycle": ("epoch", [0, 1, 2, 3, 4]), - }, - coords={ - "epoch": [ - 0, - 1, - 2, - 3, - 4, - ] - }, + "spin_start_met": [10.5, 10.5, 30.1, 30.1, 30.1], + } ) + mock_interpolate_spin_data.return_value = mock_spin_df + l1a_de = xr.Dataset( { - "shcoarse": ("epoch", [0, 1]), + "shcoarse": ("epoch", [15, 35]), "de_count": ("epoch", [2, 3]), - "met": ("epoch", [0, 1]), # MET per time epoch, not per direct event - "de_time": ("direct_event", [0000, 1000, 2000, 3000, 4000]), + "de_time": ("direct_event", [0, 1000, 2000, 3000, 4000]), }, coords={"epoch": [0, 1], "direct_event": [0, 1, 2, 3, 4]}, ) - spin = xr.Dataset( - { - "shcoarse": ("epoch", [0, 1]), - "acq_start_sec": ( - "epoch", - [20, 25], - ), - "acq_start_subsec": ( - "epoch", - [0, 0], - ), - "acq_end_sec": ( - "epoch", - [25, 30], - ), - "acq_end_subsec": ( - "epoch", - [0, 0], - ), - "num_completed": ( - "epoch", - [28, 14], - ), - } - ) + # Expected: shcoarse 15 should match spin at index 0 (10 < 15 < 20) + # shcoarse 35 should match spin at index 2 (30 < 35 < 40) + # Repeated by de_count: [2, 3] -> [index0, index0, index2, index2, index2] spin_start_times_expected = np.array( - [20, 20 + 5 / 28, 25 + 5 * 2 / 14, 25 + 5 * 3 / 14, 25 + 5 * 4 / 14] + [10.5, 10.5, 30.1, 30.1, 30.1] # 10 + 0.5e6*1e-6 # 30 + 0.1e6*1e-6 ) - spin_start_times = get_spin_start_times(l1a_de, l1b_de, spin) + # Act + spin_start_times = get_spin_start_times(l1a_de) + + # Assert np.testing.assert_allclose( spin_start_times, spin_start_times_expected, @@ -553,31 +540,44 @@ def test_get_spin_start_times(): ) -def test_set_event_met(): +@patch("imap_processing.lo.l1b.lo_l1b.interpolate_spin_data") +def test_set_event_met(mock_interpolate_spin_data): # Arrange + # Mock the spin data + mock_spin_df = pd.DataFrame( + { + "spin_start_met": [10, 10, 30, 30, 30], + } + ) + mock_interpolate_spin_data.return_value = mock_spin_df + l1b_de = xr.Dataset() l1a_de = xr.Dataset( { + "shcoarse": ("epoch", [15, 35]), "de_count": ("epoch", [2, 3]), - "de_time": ("direct_event", [0000, 1000, 2000, 3000, 4000]), + "de_time": ("direct_event", [0, 1000, 2000, 3000, 4000]), }, coords={ "epoch": [0, 1], - "direct_event": [ - 0, - 1, - 2, - 3, - 4, - ], + "direct_event": [0, 1, 2, 3, 4], }, ) - avg_spin_durations = xr.DataArray([5, 10]) - spin_start_times = xr.DataArray([10, 20, 30, 40, 50]) - expected_event_met = np.array([10, 21.2207, 34.8828, 47.3242, 59.7656]) + + # shcoarse 15 -> spin_start 10, shcoarse 35 -> spin_start 30 + # event_met = spin_start + de_time * DE_CLOCK_TICK_S + expected_event_met = np.array( + [ + 10 + 0 * DE_CLOCK_TICK_S, # 10.0 + 10 + 1000 * DE_CLOCK_TICK_S, # 14.096 + 30 + 2000 * DE_CLOCK_TICK_S, # 38.192 + 30 + 3000 * DE_CLOCK_TICK_S, # 42.288 + 30 + 4000 * DE_CLOCK_TICK_S, # 46.384 + ] + ) # Act - l1b_de = set_event_met(l1a_de, l1b_de, spin_start_times, avg_spin_durations) + l1b_de = set_event_met(l1a_de, l1b_de) # Assert np.testing.assert_allclose( @@ -586,24 +586,25 @@ def test_set_event_met(): atol=1e-4, ) - def test_set_each_event_epoch(): - l1b_de = xr.Dataset( - { - "event_met": ("epoch", [10, 20, 30, 40, 50]), - }, - coords={ - "epoch": [0, 1, 2, 3, 4], - }, - ) - epoch_expected = met_to_ttj2000ns(np.array([10, 20, 30, 40, 50])) - l1b_de = set_each_event_epoch(l1b_de) +def test_set_each_event_epoch(): + l1b_de = xr.Dataset( + { + "event_met": ("epoch", [10, 20, 30, 40, 50]), + }, + coords={ + "epoch": [0, 1, 2, 3, 4], + }, + ) + epoch_expected = met_to_ttj2000ns(np.array([10, 20, 30, 40, 50])) - np.testing.assert_allclose( - l1b_de["epoch"].values, - epoch_expected, - atol=1e-4, - ) + l1b_de = set_each_event_epoch(l1b_de) + + np.testing.assert_allclose( + l1b_de["epoch"].values, + epoch_expected, + atol=1e-4, + ) def test_set_avg_spin_durations_per_event(): @@ -838,6 +839,8 @@ def test_set_direction(mock_lo_instrument_pointing, imap_ena_sim_metakernel): ) @patch( "imap_processing.lo.l1b.lo_l1b.cartesian_to_latitudinal", + # Longitudes: -180 -> 180, 0 -> 0, 90 -> 90, 180 -> 180 + # After shift to 0-360: 180, 0, 90, 180 return_value=np.array([[0, -180, -2], [0, 0, 0], [0, 90, 1], [0, 180, 2]]), ) def test_pointing_bins(mock_cartesian_to_latitudinal, mock_frame_transform): @@ -859,7 +862,8 @@ def test_pointing_bins(mock_cartesian_to_latitudinal, mock_frame_transform): ) expected_pointing_lats = np.array([0, 20, 30, 40]) - expected_pointing_lons = np.array([0, 1800, 2700, 3600]) + # Longitude bins are now in 0-360 range after the shift + expected_pointing_lons = np.array([1800, 0, 900, 1800]) # Act l1b_de = set_pointing_bin(l1b_de) diff --git a/imap_processing/tests/spice/test_spin.py b/imap_processing/tests/spice/test_spin.py index 53a1526890..808e093dc9 100644 --- a/imap_processing/tests/spice/test_spin.py +++ b/imap_processing/tests/spice/test_spin.py @@ -305,3 +305,10 @@ def test_get_instrument_spin_phase( assert np.logical_and( 0 <= inst_phase[~expected_nan_mask], inst_phase[~expected_nan_mask] < 1 ).all() + + +def test_get_spin_start_met(fake_spin_data): + # Make sure we aren't actually interpolating and rather repeating existing values + # when we have a linear ramp of data points between spins + start_mets = spin.interpolate_spin_data(np.linspace(32, 37, 3))["spin_start_met"] + np.testing.assert_array_equal(start_mets, [30, 30, 30])