Ultra l1c pset updates (#1059)

* Ultra l1c pset updates

imap_processing/spice/geometry.py

@@ -19,6 +19,7 @@
 import numpy.typing as npt
 import pandas as pd
 import spiceypy as spice
+from numpy.typing import NDArray
 
 from imap_processing.spice.kernels import ensure_spice
 
@@ -489,3 +490,78 @@ def basis_vectors(
     ... spacecraft_z = basis_vectors[:, 2]
     """
     return np.moveaxis(get_rotation_matrix(et, from_frame, to_frame), -1, -2)
+
+
+def cartesian_to_spherical(
+    v: NDArray,
+) -> NDArray:
+    """
+    Convert cartesian coordinates to spherical coordinates.
+
+    Parameters
+    ----------
+    v : np.ndarray
+        A NumPy array with shape (n, 3) where each
+        row represents a vector
+        with x, y, z-components.
+
+    Returns
+    -------
+    spherical_coords : np.ndarray
+        A NumPy array with shape (n, 3), where each row contains
+        the spherical coordinates (r, azimuth, elevation):
+
+        - 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].
+    """
+    # Magnitude of the velocity vector
+    magnitude_v = np.linalg.norm(v, axis=-1, keepdims=True)
+
+    vhat = v / magnitude_v
+
+    # Elevation angle (angle from the z-axis, range: [-pi/2, pi/2])
+    el = np.arcsin(vhat[..., 2])
+
+    # Azimuth angle (angle in the xy-plane, range: [0, 2*pi])
+    az = np.arctan2(vhat[..., 1], vhat[..., 0])
+
+    # Ensure azimuth is from 0 to 2PI
+    az = az % (2 * np.pi)
+    spherical_coords = np.stack(
+        (np.squeeze(magnitude_v), np.degrees(az), np.degrees(el)), axis=-1
+    )
+
+    return spherical_coords
+
+
+def spherical_to_cartesian(spherical_coords: NDArray) -> NDArray:
+    """
+    Convert spherical coordinates to Cartesian coordinates.
+
+    Parameters
+    ----------
+    spherical_coords : np.ndarray
+        A NumPy array with shape (n, 3), where each row contains
+        the spherical coordinates (r, azimuth, elevation):
+
+        - 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].
+
+    Returns
+    -------
+    cartesian_coords : np.ndarray
+        Cartesian coordinates.
+    """
+    r = spherical_coords[..., 0]
+    azimuth = spherical_coords[..., 1]
+    elevation = spherical_coords[..., 2]
+
+    x = r * np.cos(elevation) * np.cos(azimuth)
+    y = r * np.cos(elevation) * np.sin(azimuth)
+    z = r * np.sin(elevation)
+
+    cartesian_coords = np.stack((x, y, z), axis=-1)
+
+    return cartesian_coords

imap_processing/tests/spice/test_data/imap_ena_sim_metakernel.template

@@ -3,4 +3,6 @@
 {SPICE_TEST_DATA_PATH}/imap_spk_demo.bsp
 {SPICE_TEST_DATA_PATH}/sim_1yr_imap_attitude.bc
 {SPICE_TEST_DATA_PATH}/imap_wkcp.tf
-{SPICE_TEST_DATA_PATH}/de440s.bsp
+{SPICE_TEST_DATA_PATH}/de440s.bsp
+{SPICE_TEST_DATA_PATH}/imap_science_0001.tf
+{SPICE_TEST_DATA_PATH}/sim_1yr_imap_pointing_frame.bc

imap_processing/tests/spice/test_geometry.py

@@ -9,6 +9,7 @@
     SpiceBody,
     SpiceFrame,
     basis_vectors,
+    cartesian_to_spherical,
     frame_transform,
     get_instrument_spin_phase,
     get_rotation_matrix,
@@ -17,6 +18,7 @@
     get_spin_data,
     imap_state,
     instrument_pointing,
+    spherical_to_cartesian,
 )
 from imap_processing.spice.kernels import ensure_spice
 
@@ -335,3 +337,52 @@ def test_basis_vectors():
                 SpiceFrame.ECLIPJ2000,
             ),
         )
+
+
+def test_cartesian_to_spherical():
+    """Tests cartesian_to_spherical function."""
+
+    step = 0.05
+    x = np.arange(-1, 1 + step, step)
+    y = np.arange(-1, 1 + step, step)
+    z = np.arange(-1, 1 + step, step)
+    x, y, z = np.meshgrid(x, y, z)
+
+    cartesian_points = np.stack((x.ravel(), y.ravel(), z.ravel()), axis=-1)
+
+    for point in cartesian_points:
+        r, az, el = cartesian_to_spherical(point)
+        r_spice, colat_spice, slong_spice = spice.recsph(point)
+
+        # Convert SPICE co-latitude to elevation
+        el_spice = 90 - np.degrees(colat_spice)
+        az_spice = np.degrees(slong_spice)
+
+        # Normalize azimuth to [0, 360]
+        az_spice = az_spice % 360
+
+        np.testing.assert_allclose(r, r_spice, atol=1e-5)
+        np.testing.assert_allclose(az, az_spice, atol=1e-5)
+        np.testing.assert_allclose(el, el_spice, atol=1e-5)
+
+
+def test_spherical_to_cartesian():
+    """Tests spherical_to_cartesian function."""
+
+    azimuth = np.linspace(0, 2 * np.pi, 50)
+    elevation = np.linspace(-np.pi / 2, np.pi / 2, 50)
+    theta, elev = np.meshgrid(azimuth, elevation)
+    r = 1.0
+
+    spherical_points = np.stack(
+        (r * np.ones_like(theta).ravel(), theta.ravel(), elev.ravel()), axis=-1
+    )
+
+    # Convert elevation to colatitude for SPICE
+    colat = np.pi / 2 - spherical_points[:, 2]
+
+    for i in range(len(colat)):
+        cartesian_coords = spherical_to_cartesian(np.array([spherical_points[i]]))
+        spice_coords = spice.sphrec(r, colat[i], spherical_points[i, 1])
+
+        np.testing.assert_allclose(cartesian_coords[0], spice_coords, atol=1e-5)

imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py

@@ -3,14 +3,12 @@
 import cdflib
 import numpy as np
 import pytest
-import spiceypy as spice
 from cdflib import CDF
 
 from imap_processing import imap_module_directory
 from imap_processing.ultra.l1c.ultra_l1c_pset_bins import (
     build_energy_bins,
     build_spatial_bins,
-    cartesian_to_spherical,
     get_helio_exposure_times,
     get_histogram,
     get_pointing_frame_exposure_times,
@@ -32,33 +30,15 @@ def test_data():
     return v, energy
 
 
-@pytest.fixture()
-def kernels(spice_test_data_path):
-    """List SPICE kernels."""
-    required_kernels = [
-        "imap_science_0001.tf",
-        "imap_sclk_0000.tsc",
-        "sim_1yr_imap_attitude.bc",
-        "imap_wkcp.tf",
-        "naif0012.tls",
-        "sim_1yr_imap_pointing_frame.bc",
-        "de440s.bsp",
-        "imap_spk_demo.bsp",
-    ]
-    kernels = [str(spice_test_data_path / kernel) for kernel in required_kernels]
-
-    return kernels
-
-
 def test_build_energy_bins():
     """Tests build_energy_bins function."""
     energy_bin_edges, energy_midpoints = build_energy_bins()
-    energy_bin_start = energy_bin_edges[:-1]
-    energy_bin_end = energy_bin_edges[1:]
+    energy_bin_start = [interval[0] for interval in energy_bin_edges]
+    energy_bin_end = [interval[1] for interval in energy_bin_edges]
 
     assert energy_bin_start[0] == 0
     assert energy_bin_start[1] == 3.385
-    assert len(energy_bin_edges) == 25
+    assert len(energy_bin_edges) == 24
     assert energy_midpoints[0] == (energy_bin_start[0] + energy_bin_end[0]) / 2
 
     # Comparison to expected values.
@@ -90,21 +70,6 @@ def test_build_spatial_bins():
     np.testing.assert_allclose(el_bin_midpoints[-1], 89.75, atol=1e-4)
 
 
-def test_cartesian_to_spherical(test_data):
-    """Tests cartesian_to_spherical function."""
-    v, _ = test_data
-
-    az_sc, el_sc, r = cartesian_to_spherical(v)
-
-    # MATLAB code outputs:
-    np.testing.assert_allclose(
-        np.unique(np.radians(az_sc)), np.array([1.31300, 2.34891]), atol=1e-05, rtol=0
-    )
-    np.testing.assert_allclose(
-        np.unique(np.radians(el_sc)), np.array([-0.88901, -0.70136]), atol=1e-05, rtol=0
-    )
-
-
 def test_get_histogram(test_data):
     """Tests get_histogram function."""
     v, energy = test_data
@@ -119,7 +84,7 @@ def test_get_histogram(test_data):
     assert hist.shape == (
         len(az_bin_edges) - 1,
         len(el_bin_edges) - 1,
-        len(energy_bin_edges) - 1,
+        len(energy_bin_edges),
     )
 
 
@@ -143,10 +108,10 @@ def test_get_pointing_frame_exposure_times():
 
 
 @pytest.mark.external_kernel()
-def test_et_helio_exposure_times(kernels):
+@pytest.mark.use_test_metakernel("imap_ena_sim_metakernel.template")
+def test_et_helio_exposure_times():
     """Tests get_helio_exposure_times function."""
 
-    spice.furnsh(kernels)
     constant_exposure = BASE_PATH / "dps_grid45_compressed.cdf"
     start_time = 829485054.185627
     end_time = 829567884.185627
@@ -185,9 +150,9 @@ def test_et_helio_exposure_times(kernels):
             transposed_exposure = np.transpose(exposure_data, (2, 1, 0))
             exposures.append(transposed_exposure)
 
-    np.array_equal(exposures[0], exposure_3d[:, :, 0])
-    np.array_equal(exposures[1], exposure_3d[:, :, 1])
-    np.array_equal(exposures[2], exposure_3d[:, :, 2])
+    assert np.array_equal(np.squeeze(exposures[0]), exposure_3d[:, :, 0])
+    assert np.array_equal(np.squeeze(exposures[1]), exposure_3d[:, :, 11])
+    assert np.array_equal(np.squeeze(exposures[2]), exposure_3d[:, :, 23])
 
 
 def test_get_pointing_frame_sensitivity():

imap_processing/ultra/constants.py

@@ -62,3 +62,8 @@ class UltraConstants:
     # Conversion factors
     KEV_J = 1.602180000000000e-16  # 1.6021766339999e-16 # keV to joules
     MASS_H = 1.6735575e-27  # Mass of a hydrogen atom in kilograms.
+
+    # Energy bin constants
+    ALPHA = 0.2  # deltaE/E
+    ENERGY_START = 3.385  # energy start for the Ultra grids
+    N_BINS = 23  # number of energy bins