Ultra l1b extended energy (#825)

* added functions to calculate energy

imap_processing/tests/ultra/unit/test_lookup_utils.py

@@ -57,7 +57,7 @@ def test_get_egy_norm():
 
     norm_composite_energy = get_energy_norm(np.array([2]), np.array([2]))
 
-    assert norm_composite_energy == egy_norm_df.iloc[2 * 4096 + 2]["NormEnergy"]
+    assert int(norm_composite_energy) == egy_norm_df.iloc[2 * 4096 + 2]["NormEnergy"]
 
 
 def test_get_image_params():

imap_processing/tests/ultra/unit/test_ultra_l1b_extended.py

@@ -10,6 +10,8 @@
     StopType,
     calculate_etof_xc,
     get_coincidence_positions,
+    get_energy_pulse_height,
+    get_energy_ssd,
     get_front_x_position,
     get_front_y_position,
     get_particle_velocity,
@@ -257,3 +259,33 @@ def test_get_ssd_tof(de_dataset, yf_fixture):
     np.testing.assert_allclose(
         ssd_tof, df_ssd["TOF"].astype("float"), atol=1e-05, rtol=0
     )
+
+
+def test_get_energy_ssd(de_dataset, yf_fixture):
+    """Tests get_energy_ssd function."""
+    df_filt, _, _ = yf_fixture
+    df_ssd = df_filt[df_filt["StopType"].isin(StopType.SSD.value)]
+    _, _, ssd_number = get_ssd_back_position_and_tof_offset(de_dataset)
+    energy = get_energy_ssd(de_dataset, ssd_number)
+    test_energy = df_ssd["Energy"].astype("float")
+
+    assert np.array_equal(test_energy, energy)
+
+
+def test_get_energy_pulse_height(de_dataset, yf_fixture):
+    """Tests get_energy_ssd function."""
+    df_filt, _, _ = yf_fixture
+    df_ph = df_filt[df_filt["StopType"].isin(StopType.PH.value)]
+    ph_indices = np.nonzero(
+        np.isin(de_dataset["STOP_TYPE"], [StopType.Top.value, StopType.Bottom.value])
+    )[0]
+
+    test_xb = df_filt["Xb"].astype("float").values
+    test_yb = df_filt["Yb"].astype("float").values
+
+    energy = get_energy_pulse_height(
+        de_dataset["STOP_TYPE"].data, de_dataset["ENERGY_PH"].data, test_xb, test_yb
+    )
+    test_energy = df_ph["Energy"].astype("float")
+
+    assert np.array_equal(test_energy, energy[ph_indices])

imap_processing/ultra/l1b/lookup_utils.py

@@ -140,7 +140,7 @@ def get_energy_norm(ssd: np.ndarray, composite_energy: np.ndarray) -> npt.NDArra
     """
     row_number = ssd * 4096 + composite_energy
 
-    return _ENERGY_NORM_DF["NormEnergy"].values[row_number]
+    return _ENERGY_NORM_DF["NormEnergy"].iloc[row_number]
 
 
 def get_image_params(image: str) -> np.float64:

imap_processing/ultra/l1b/ultra_l1b_extended.py

@@ -11,6 +11,7 @@
 
 from imap_processing.ultra.l1b.lookup_utils import (
     get_back_position,
+    get_energy_norm,
     get_image_params,
     get_norm,
     get_y_adjust,
@@ -41,6 +42,7 @@ class StopType(Enum):
 
     Top = 1
     Bottom = 2
+    PH: ClassVar[list[int]] = [1, 2]
     SSD: ClassVar[list[int]] = [8, 9, 10, 11, 12, 13, 14, 15]
 
 
@@ -541,3 +543,109 @@ def get_ssd_tof(de_dataset: xarray.Dataset, xf: np.ndarray) -> NDArray[np.float6
 
     # Convert TOF to tenths of a nanosecond.
     return np.asarray(tof, dtype=np.float64)
+
+
+def get_energy_pulse_height(
+    stop_type: np.ndarray, energy: np.ndarray, xb: np.ndarray, yb: np.ndarray
+) -> NDArray[np.float64]:
+    """
+    Calculate the pulse-height energy.
+
+    Calculate energy measured using the
+    pulse height from the stop anode.
+    Lookup tables (lut) are used for corrections.
+    Further description is available on pages 40-41 of
+    IMAP-Ultra Flight Software Specification document
+    (7523-9009_Rev_-.pdf).
+
+    Parameters
+    ----------
+    stop_type : np.ndarray
+        Stop type: 1=Top, 2=Bottom.
+    energy : np.ndarray
+        Energy measured using the pulse height.
+    xb : np.ndarray
+        X back position (hundredths of a millimeter).
+    yb : np.ndarray
+        Y back position (hundredths of a millimeter).
+
+    Returns
+    -------
+    energy_ph : np.ndarray
+        Energy measured using the pulse height
+        from the stop anode (DN).
+    """
+    indices_top = np.where(stop_type == 1)[0]
+    indices_bottom = np.where(stop_type == 2)[0]
+
+    xlut = np.zeros(len(stop_type), dtype=np.float64)
+    ylut = np.zeros(len(stop_type), dtype=np.float64)
+    energy_ph = np.zeros(len(stop_type), dtype=np.float64)
+
+    # Stop type 1
+    xlut[indices_top] = (xb[indices_top] / 100 - 25 / 2) * 20 / 50  # mm
+    ylut[indices_top] = (yb[indices_top] / 100 + 82 / 2) * 32 / 82  # mm
+    # Stop type 2
+    xlut[indices_bottom] = (xb[indices_bottom] / 100 + 50 + 25 / 2) * 20 / 50  # mm
+    ylut[indices_bottom] = (yb[indices_bottom] / 100 + 82 / 2) * 32 / 82  # mm
+
+    # TODO: waiting on these lookup tables: SpTpPHCorr, SpBtPHCorr
+    energy_ph[indices_top] = energy[indices_top] - get_image_params(
+        "SPTPPHOFF"
+    )  # * SpTpPHCorr[
+    # xlut[indices_top], ylut[indices_top]] / 1024
+
+    energy_ph[indices_bottom] = energy[indices_bottom] - get_image_params(
+        "SPBTPHOFF"
+    )  # * SpBtPHCorr[
+    # xlut[indices_bottom], ylut[indices_bottom]] / 1024
+
+    return energy_ph
+
+
+def get_energy_ssd(de_dataset: xarray.Dataset, ssd: np.ndarray) -> NDArray[np.float64]:
+    """
+    Get SSD energy.
+
+    For SSD events, the SSD itself provides a direct
+    measurement of the energy. To cover higher energies,
+    a so-called composite energy is calculated using the
+    SSD energy and SSD energy pulse width.
+    The result is then normalized per SSD via a lookup table.
+    Further description is available on pages 41 of
+    IMAP-Ultra Flight Software Specification document
+    (7523-9009_Rev_-.pdf).
+
+    Parameters
+    ----------
+    de_dataset : xarray.Dataset
+        Events dataset.
+    ssd : np.ndarray
+        SSD number.
+
+    Returns
+    -------
+    energy_norm : np.ndarray
+        Energy measured using the SSD.
+    """
+    # DN threshold for composite energy.
+    composite_energy_threshold = 1707
+
+    ssd_indices = np.where(de_dataset["STOP_TYPE"].data >= 8)[0]
+    energy = de_dataset["ENERGY_PH"].data[ssd_indices]
+
+    composite_energy = np.empty(len(energy), dtype=np.float64)
+
+    composite_energy[energy >= composite_energy_threshold] = (
+        composite_energy_threshold
+        + de_dataset["PULSE_WIDTH"].data[ssd_indices][
+            energy >= composite_energy_threshold
+        ]
+    )
+    composite_energy[energy < composite_energy_threshold] = energy[
+        energy < composite_energy_threshold
+    ]
+
+    energy_norm = get_energy_norm(ssd, composite_energy)
+
+    return energy_norm

imap_processing/ultra/lookup_tables/FM45_Startup1_ULTRA_IMGPARAMS_20240719.csv

@@ -1,2 +1,2 @@
 SHCOARSE,XFTSC,XFTLTOFF,XFTRTOFF,TOFSC,TOFTPOFF,TOFBTOFF,XFTTOF,XCOINTPSC,XCOINTPOFF,XCOINBTSC,XCOINBTOFF,ETOFSC,ETOFTPOFF,ETOFBTOFF,TOFDIFFTPMIN,TOFDIFFTPMAX,TOFDIFFBTMIN,TOFDIFFBTMAX,ETOFMIN,ETOFMAX,ETOFSLOPE1,ETOFOFF1,ETOFSLOPE2,ETOFOFF2,SPTPPHOFF,SPBTPHOFF,YBKSSD0,YBKSSD1,YBKSSD2,YBKSSD3,YBKSSD4,YBKSSD5,YBKSSD6,YBKSSD7,TOFSSDSC,TOFSSDLTOFF0,TOFSSDLTOFF1,TOFSSDLTOFF2,TOFSSDLTOFF3,TOFSSDLTOFF4,TOFSSDLTOFF5,TOFSSDLTOFF6,TOFSSDLTOFF7,TOFSSDRTOFF0,TOFSSDRTOFF1,TOFSSDRTOFF2,TOFSSDRTOFF3,TOFSSDRTOFF4,TOFSSDRTOFF5,TOFSSDRTOFF6,TOFSSDRTOFF7,TOFSSDTOTOFF,PATHSTEEPTHRESH,PATHMEDIUMTHRESH
-,0.172998047,49.3,48.25,0.5,-528,-525,0.01831055,0.067929688,41.75,0.067929688,-40.75,0.1,-44.5,-44.5,,,,,,,,,,,,,29.3,37.3,7.1,15.1,-15.1,-7.1,-37.3,-29.3,0.196484375,-6,-7.3,-3.8,-4.2,-3.8,-3.7,-6.3,-5,-5,-6.3,-3.7,-3.8,-4,-4.2,-7.3,-6,5.9,,
+,0.172998047,49.3,48.25,0.5,-528,-525,0.01831055,0.067929688,41.75,0.067929688,-40.75,0.1,-44.5,-44.5,,,,,,,,,,,580,580,29.3,37.3,7.1,15.1,-15.1,-7.1,-37.3,-29.3,0.196484375,-6,-7.3,-3.8,-4.2,-3.8,-3.7,-6.3,-5,-5,-6.3,-3.7,-3.8,-4,-4.2,-7.3,-6,5.9,,