1148 hi l1c populate despun z (#1253)

* Compute despun_z vector in Hi L1C processing
Add test coverage for pset_geometry function

* Fix docstring in geometry.py cartesian to spherical and spherical to cartesian functions
Calculate despun_z vector and hae lat/lon of spin bins in L1C pset
Add test coverage for hae lat/lon values

* Fix bug where J2000ns were being passed to spiceypy

* Address PR comments

imap_processing/hi/l1c/hi_l1c.py

@@ -8,6 +8,13 @@
 from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes
 from imap_processing.cdf.utils import parse_filename_like
 from imap_processing.hi.utils import create_dataset_variables, full_dataarray
+from imap_processing.spice.geometry import (
+    SpiceFrame,
+    cartesian_to_spherical,
+    frame_transform,
+    spherical_to_cartesian,
+)
+from imap_processing.spice.time import j2000ns_to_j2000s
 
 logger = logging.getLogger(__name__)
 
@@ -72,8 +79,10 @@ def generate_pset_dataset(de_dataset: xr.Dataset) -> xr.Dataset:
     pset_dataset.epoch.data[0] = np.mean(de_dataset.epoch.data[[0, -1]]).astype(
         np.int64
     )
-
-    pset_dataset.update(pset_geometry())
+    pset_et = j2000ns_to_j2000s(pset_dataset.epoch.data[0])
+    # Calculate and add despun_z, hae_latitude, and hae_longitude variables to
+    # the pset_dataset
+    pset_dataset.update(pset_geometry(pset_et, logical_source_parts["sensor"]))
 
     # TODO: The following section will go away as PSET algorithms to populate
     #    these variables are written.
@@ -204,10 +213,17 @@ def empty_pset_dataset(n_esa_steps: int, sensor_str: str) -> xr.Dataset:
     return dataset
 
 
-def pset_geometry() -> dict[str, xr.DataArray]:
+def pset_geometry(pset_et: float, sensor_str: str) -> dict[str, xr.DataArray]:
     """
     Calculate PSET geometry variables.
 
+    Parameters
+    ----------
+    pset_et : float
+        Pointing set ephemeris time for which to calculate PSET geometry.
+    sensor_str : str
+        '45sensor' or '90sensor'.
+
     Returns
     -------
     geometry_vars : dict[str, xarray.DataArray]
@@ -216,13 +232,43 @@ def pset_geometry() -> dict[str, xr.DataArray]:
     geometry_vars = create_dataset_variables(
         ["despun_z"], (1, 3), att_manager_lookup_str="hi_pset_{0}"
     )
-    # TODO: Calculate despun_z
+    despun_z = frame_transform(
+        pset_et,
+        np.array([0, 0, 1]),
+        SpiceFrame.IMAP_DPS,
+        SpiceFrame.ECLIPJ2000,
+    )
+    geometry_vars["despun_z"].values = despun_z[np.newaxis, :].astype(np.float32)
+
+    # Calculate hae_latitude and hae_longitude of the spin bins
+    # define the azimuth/elevation coordinates in the Pointing Frame (DPS)
+    # TODO: get the sensor's true elevation using SPICE?
+    el = 0 if "90" in sensor_str else -45
+    dps_az_el = np.array(
+        [
+            np.ones(3600),
+            np.deg2rad(np.arange(0.05, 360, 0.1)),
+            np.full(3600, np.deg2rad(el)),
+        ]
+    ).T
+    dps_cartesian = spherical_to_cartesian(dps_az_el)
+    # Transform DPS Cartesian coords into HAE Ecliptic
+    hae_eclip_cartesian = frame_transform(
+        pset_et, dps_cartesian, SpiceFrame.IMAP_DPS, SpiceFrame.ECLIPJ2000
+    )
+    hae_az_el = cartesian_to_spherical(hae_eclip_cartesian, degrees=True)
+
     geometry_vars.update(
         create_dataset_variables(
             ["hae_latitude", "hae_longitude"],
             (1, 3600),
             att_manager_lookup_str="hi_pset_{0}",
         )
     )
-    # TODO: Calculate HAE Lat/Lon
+    geometry_vars["hae_longitude"].values = hae_az_el[:, 1].astype(np.float32)[
+        np.newaxis, :
+    ]
+    geometry_vars["hae_latitude"].values = hae_az_el[:, 2].astype(np.float32)[
+        np.newaxis, :
+    ]
     return geometry_vars

imap_processing/spice/geometry.py

@@ -524,7 +524,7 @@ def cartesian_to_spherical(
           output range=[0, 360],
           otherwise in radians if degrees parameter is False:
           output range=[0, 2*pi].
-        - elevation : angle from the z-axis
+        - elevation : angle from the xy-plane
           In degrees if degrees parameter is True (by default):
           output range=[0, 180],
           otherwise in radians if degrees parameter is False:
@@ -535,7 +535,7 @@ def cartesian_to_spherical(
 
     vhat = v / magnitude_v
 
-    # Elevation angle (angle from the z-axis, range: [-pi/2, pi/2])
+    # Elevation angle (angle from the xy-plane, range: [-pi/2, pi/2])
     el = np.arcsin(vhat[..., 2])
 
     # Azimuth angle (angle in the xy-plane, range: [0, 2*pi])
@@ -565,7 +565,7 @@ def spherical_to_cartesian(spherical_coords: NDArray) -> NDArray:
 
         - r : Distance of the point from the origin.
         - azimuth : angle in the xy-plane in radians [0, 2*pi].
-        - elevation : angle from the z-axis in radians [-pi/2, pi/2].
+        - elevation : angle from the xy-plane in radians [-pi/2, pi/2].
 
     Returns
     -------

imap_processing/tests/hi/test_hi_l1c.py

@@ -1,12 +1,17 @@
 """Test coverage for imap_processing.hi.l1c.hi_l1c.py"""
 
+from unittest import mock
+
 import numpy as np
+import pytest
 
 from imap_processing.cdf.utils import load_cdf, write_cdf
 from imap_processing.hi.l1c import hi_l1c
 from imap_processing.hi.utils import HIAPID
 
 
+@pytest.mark.external_kernel()
+@pytest.mark.use_test_metakernel("imap_ena_sim_metakernel.template")
 def test_generate_pset_dataset(hi_l1_test_data_path):
     """Test coverage for generate_pset_dataset function"""
     l1b_de_path = hi_l1_test_data_path / "imap_hi_l1b_45sensor-de_20250415_v999.cdf"
@@ -44,3 +49,31 @@ def test_empty_pset_dataset():
     assert dataset.spin_angle_bin.size == 3600
     assert dataset.esa_energy_step.size == n_esa_steps
     assert dataset.calibration_prod.size == n_calibration_prods
+
+
+@pytest.mark.parametrize("sensor_str", ["90sensor", "45sensor"])
+@mock.patch("imap_processing.hi.l1c.hi_l1c.frame_transform")
+def test_pset_geometry(mock_frame_transform, sensor_str):
+    """Test coverage for pset_geometry function"""
+    # Mock frame transform to simply return the input position vectors (no transform)
+    mock_frame_transform.side_effect = lambda et, pos, from_frame, to_frame: pos
+
+    geometry_vars = hi_l1c.pset_geometry(0, sensor_str)
+
+    assert "despun_z" in geometry_vars
+    np.testing.assert_array_equal(geometry_vars["despun_z"].data, [[0, 0, 1]])
+
+    assert "hae_latitude" in geometry_vars
+    assert "hae_longitude" in geometry_vars
+    # frame_transform is mocked to return the input vectors. For Hi-90, we
+    # expect hae_latitude to be 0, and for Hi-45 we expect -45. Both sensors
+    # have an expected longitude to be 0.1 degree steps starting at 0.05
+    expected_latitude = 0 if sensor_str == "90sensor" else -45
+    np.testing.assert_array_equal(
+        geometry_vars["hae_latitude"].data, np.full((1, 3600), expected_latitude)
+    )
+    np.testing.assert_allclose(
+        geometry_vars["hae_longitude"].data,
+        np.arange(0.05, 360, 0.1, dtype=np.float32).reshape((1, 3600)),
+        atol=4e-05,
+    )