I-ALiRT - MAG l1b equivalent (IMAP-Science-Operations-Center#1661)

Author: laspsandoval

imap_processing/ialirt/l0/parse_mag.py

@@ -1,6 +1,7 @@
 """Functions to support I-ALiRT MAG packet parsing."""
 
 import logging
+from typing import Union
 
 import numpy as np
 import xarray as xr
@@ -13,6 +14,13 @@
 )
 from imap_processing.ialirt.utils.grouping import find_groups
 from imap_processing.ialirt.utils.time import calculate_time
+from imap_processing.mag.l1a.mag_l1a_data import TimeTuple
+from imap_processing.mag.l1b.mag_l1b import (
+    calibrate_vector,
+    retrieve_matrix_from_l1b_calibration,
+    shift_time,
+)
+from imap_processing.spice.time import met_to_ttj2000ns
 
 logger = logging.getLogger(__name__)
 
@@ -132,7 +140,13 @@ def extract_magnetic_vectors(science_values: xr.DataArray) -> dict:
     return vectors
 
 
-def get_time(grouped_data: xr.Dataset, group: int, pkt_counter: xr.DataArray) -> dict:
+def get_time(
+    grouped_data: xr.Dataset,
+    group: int,
+    pkt_counter: xr.DataArray,
+    time_shift_mago: xr.DataArray,
+    time_shift_magi: xr.DataArray,
+) -> dict:
     """
     Get the time for the grouped data.
 
@@ -144,12 +158,22 @@ def get_time(grouped_data: xr.Dataset, group: int, pkt_counter: xr.DataArray) ->
         Group number.
     pkt_counter : xr.DataArray
         Packet counter.
+    time_shift_mago : xr.DataArray
+        Time shift value mago.
+    time_shift_magi : xr.DataArray
+        Time shift value magi.
 
     Returns
     -------
     time_data : dict
         Coarse and fine time for Primary and Secondary Sensors.
+
+    Notes
+    -----
+    Packet id 0 is course and fine time for the primary sensor PRI.
+    Packet id 2 is the course time for the secondary sensor SEC.
     """
+    # Get the coarse and fine time for the primary and secondary sensors.
     pri_coarsetm = grouped_data["mag_acq_tm_coarse"][
         (grouped_data["group"] == group).values
     ][pkt_counter == 0]
@@ -166,24 +190,111 @@ def get_time(grouped_data: xr.Dataset, group: int, pkt_counter: xr.DataArray) ->
         (grouped_data["group"] == group).values
     ][pkt_counter == 2]
 
-    time_data = {
-        "pri_coarsetm": int(pri_coarsetm),
-        "pri_fintm": int(pri_fintm),
-        "sec_coarsetm": int(sec_coarsetm),
-        "sec_fintm": int(sec_fintm),
+    time_data: dict[str, Union[int, float]] = {
+        "pri_coarsetm": int(pri_coarsetm.item()),
+        "pri_fintm": int(pri_fintm.item()),
+        "sec_coarsetm": int(sec_coarsetm.item()),
+        "sec_fintm": int(sec_fintm.item()),
     }
 
+    primary_time = TimeTuple(int(pri_coarsetm.item()), int(pri_fintm.item()))
+    secondary_time = TimeTuple(int(sec_coarsetm.item()), int(sec_fintm.item()))
+    time_data_pri_met = primary_time.to_seconds()
+    time_data_primary_ttj2000ns = met_to_ttj2000ns(time_data_pri_met)
+    time_data["primary_epoch"] = shift_time(
+        time_data_primary_ttj2000ns, time_shift_mago
+    )
+    time_data_sec_met = secondary_time.to_seconds()
+    time_data_secondary_ttj2000ns = met_to_ttj2000ns(time_data_sec_met)
+    time_data["secondary_epoch"] = shift_time(
+        time_data_secondary_ttj2000ns, time_shift_magi
+    )
+
     return time_data
 
 
-def parse_packet(accumulated_data: xr.Dataset) -> list[dict]:
+def calculate_l1b(
+    grouped_data: xr.Dataset,
+    group: int,
+    pkt_counter: xr.DataArray,
+    science_data: dict,
+    status_data: dict,
+    calibration_dataset: xr.Dataset,
+) -> tuple[np.ndarray, np.ndarray, dict]:
+    """
+    Calculate equivalent of l1b data product.
+
+    Parameters
+    ----------
+    grouped_data : xr.Dataset
+        Grouped data.
+    group : int
+        Group number.
+    pkt_counter : xr.DataArray
+        Packet counter.
+    science_data : dict
+        Science data.
+    status_data : dict
+        Status data.
+    calibration_dataset : xr.Dataset
+        Calibration dataset.
+
+    Returns
+    -------
+    updated_vector_mago : numpy.ndarray
+        Calibrated mago vector.
+    updated_vector_magi : numpy.ndarray
+        Calibrated magi vector.
+    time_data : dict
+        Time data.
+    """
+    calibration_matrix_mago, time_shift_mago = retrieve_matrix_from_l1b_calibration(
+        calibration_dataset, is_mago=True
+    )
+    calibration_matrix_magi, time_shift_magi = retrieve_matrix_from_l1b_calibration(
+        calibration_dataset, is_mago=False
+    )
+
+    # Get time values for each group.
+    time_data = get_time(
+        grouped_data, group, pkt_counter, time_shift_mago, time_shift_magi
+    )
+
+    input_vector_mago = np.array(
+        [
+            science_data["pri_x"],
+            science_data["pri_y"],
+            science_data["pri_z"],
+            status_data["fob_range"],
+        ]
+    )
+    input_vector_magi = np.array(
+        [
+            science_data["sec_x"],
+            science_data["sec_y"],
+            science_data["sec_z"],
+            status_data["fib_range"],
+        ]
+    )
+
+    updated_vector_mago = calibrate_vector(input_vector_mago, calibration_matrix_mago)
+    updated_vector_magi = calibrate_vector(input_vector_magi, calibration_matrix_magi)
+
+    return updated_vector_mago, updated_vector_magi, time_data
+
+
+def process_packet(
+    accumulated_data: xr.Dataset, calibration_dataset: xr.Dataset
+) -> list[dict]:
     """
     Parse the MAG packets.
 
     Parameters
     ----------
     accumulated_data : xr.Dataset
         Packets dataset accumulated over 1 min.
+    calibration_dataset : xr.Dataset
+        Calibration dataset.
 
     Returns
     -------
@@ -237,9 +348,26 @@ def parse_packet(accumulated_data: xr.Dataset) -> list[dict]:
             (grouped_data["group"] == group).values
         ]
         science_data = extract_magnetic_vectors(science_values)
-
-        # Get time values for each group.
-        time_data = get_time(grouped_data, group, pkt_counter)
+        updated_vector_mago, updated_vector_magi, time_data = calculate_l1b(
+            grouped_data,
+            group,
+            pkt_counter,
+            science_data,
+            status_data,
+            calibration_dataset,
+        )
+
+        # Note: primary = MAGo, secondary = MAGi.
+        science_data.update(
+            {
+                "calibrated_pri_x": updated_vector_mago[0],
+                "calibrated_pri_y": updated_vector_mago[1],
+                "calibrated_pri_z": updated_vector_mago[2],
+                "calibrated_sec_x": updated_vector_magi[0],
+                "calibrated_sec_y": updated_vector_magi[1],
+                "calibrated_sec_z": updated_vector_magi[2],
+            }
+        )
 
         mag_data.append({**status_data, **science_data, **time_data})
 

imap_processing/mag/l1b/mag_l1b.py

@@ -369,7 +369,7 @@ def shift_time(epoch_times: xr.DataArray, time_shift: xr.DataArray) -> xr.DataAr
     Parameters
     ----------
     epoch_times : xr.DataArray
-        The input epoch times, in J2000 ns.
+        The input epoch times, in TT J2000 ns.
     time_shift : xr.DataArray
         The time shift to apply for the given sensor. This should be one value and is
         in seconds.

imap_processing/tests/ialirt/unit/test_parse_mag.py

@@ -1,15 +1,22 @@
+"""Tests to support I-ALiRT MAG packet parsing."""
+
 import numpy as np
 import pandas as pd
 import pytest
 import xarray as xr
 
 from imap_processing import imap_module_directory
+from imap_processing.cdf.utils import load_cdf
 from imap_processing.ialirt.l0.parse_mag import (
+    calculate_l1b,
     extract_magnetic_vectors,
     get_pkt_counter,
     get_status_data,
     get_time,
-    parse_packet,
+    process_packet,
+)
+from imap_processing.mag.l1b.mag_l1b import (
+    retrieve_matrix_from_l1b_calibration,
 )
 from imap_processing.utils import packet_file_to_datasets
 
@@ -107,6 +114,15 @@ def grouped_data():
     return grouped_data
 
 
+@pytest.fixture
+def calibration_dataset():
+    """Returns the calibration data."""
+    calibration_dataset = load_cdf(
+        imap_module_directory / "mag" / "l1b" / "imap_calibration_mag_20240229_v01.cdf"
+    )
+    return calibration_dataset
+
+
 def test_get_pkt_counter(xarray_data):
     """Tests the get_pkt_counter function."""
     status_values = xarray_data["mag_status"].values
@@ -127,15 +143,24 @@ def test_get_status_data(xarray_data, mag_test_data):
         assert status_data[key] == matching_row[key.upper()].values[0]
 
 
-def test_get_time(grouped_data):
+def test_get_time(grouped_data, calibration_dataset):
     """Tests the get_time function."""
-    time_data = get_time(grouped_data, 1, np.array([0, 1, 2, 3]))
-    assert time_data == {
-        "pri_coarsetm": 461971386,
-        "pri_fintm": 1500,
-        "sec_coarsetm": 461971386,
-        "sec_fintm": 1503,
-    }
+
+    calibration_matrix_mago, time_shift_mago = retrieve_matrix_from_l1b_calibration(
+        calibration_dataset, is_mago=True
+    )
+    calibration_matrix_magi, time_shift_magi = retrieve_matrix_from_l1b_calibration(
+        calibration_dataset, is_mago=False
+    )
+
+    time_data = get_time(
+        grouped_data, 1, np.array([0, 1, 2, 3]), time_shift_mago, time_shift_magi
+    )
+
+    assert time_data["pri_coarsetm"] == 461971386
+    assert time_data["pri_fintm"] == 1500
+    assert time_data["sec_coarsetm"] == 461971386
+    assert time_data["sec_fintm"] == 1503
 
 
 def test_extract_magnetic_vectors():
@@ -156,13 +181,43 @@ def test_extract_magnetic_vectors():
     }
 
 
-def test_parse_packet(xarray_data, mag_test_data):
+def test_calculate_l1b(grouped_data, xarray_data, calibration_dataset):
+    """Tests the calculate_l1b function."""
+
+    pkt_counter = np.array([0.0, 1.0, 2.0, 3.0])
+
+    science_data = {
+        "pri_x": 1.0,
+        "pri_y": 2.0,
+        "pri_z": 3.0,
+        "sec_x": 4.0,
+        "sec_y": 5.0,
+        "sec_z": 6.0,
+    }
+
+    status_data = {
+        "fob_range": 1,
+        "fib_range": 1,
+    }
+
+    vec_mago, vec_magi, time_data = calculate_l1b(
+        grouped_data, 0, pkt_counter, science_data, status_data, calibration_dataset
+    )
+
+    assert vec_mago.shape == (4,)
+    assert vec_magi.shape == (4,)
+    assert "primary_epoch" in time_data
+    assert "secondary_epoch" in time_data
+
+
+def test_process_packet(xarray_data, mag_test_data, calibration_dataset):
     """Tests the parse_packet function."""
-    parsed_packets = parse_packet(xarray_data)
+    parsed_packets = process_packet(xarray_data, calibration_dataset)
 
     for packet in parsed_packets:
         index = packet["pri_coarsetm"] == mag_test_data["PRI_COARSETM"]
         matching_rows = mag_test_data[index]
 
         for key in packet.keys():
-            assert packet[key] == matching_rows[key.upper()].values[0]
+            if key.upper() in matching_rows.keys():
+                assert packet[key] == matching_rows[key.upper()].values[0]