FIX: Lo update spin time calculation

- use universal spin table to get true onboard timing
  The Lo spin packet is not accurate enough or has different
  time information in it, so we need to reference DE events to
  the onboard timekeeping (spin table), which the Lo instrument
  references to reset its counter.

- Change formula for DE time tag counter. It is a fixed length
  pulse of 4.096ms, not varying with spin period. This means we
  need to multiply each DE tick by that length of time.

Author: greglucas

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)
 

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,94 +502,82 @@ 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,
         atol=1e-4,
     )
 
 
-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)