From c423dc8ac3ce51ec1700aaae5f572dfd98cb2b13 Mon Sep 17 00:00:00 2001
From: daralynnrhode <daralynn.rhode@lasp.colorado.edu>
Date: Mon, 19 Oct 2026 10:02:58 -0600
Subject: [PATCH] Use geometric mean for HIT energy bin centers

---
 imap_processing/hit/hit_utils.py            | 10 +++++--
 imap_processing/tests/hit/test_hit_utils.py | 31 +++++++++++++--------
 2 files changed, 28 insertions(+), 13 deletions(-)

diff --git a/imap_processing/hit/hit_utils.py b/imap_processing/hit/hit_utils.py
index 5fb0c56c1..3930c9335 100644
--- a/imap_processing/hit/hit_utils.py
+++ b/imap_processing/hit/hit_utils.py
@@ -356,7 +356,12 @@ def add_energy_variables(
     attr_mgr: ImapCdfAttributes | None = None,
 ) -> xr.Dataset:
     """
-    Add energy min and max variables to the dataset.
+    Add energy mean and delta variables to the dataset.
+
+    The energy mean for each energy range is the geometric mean of the
+    minimum and maximum energies, sqrt(energy_min * energy_max), rounded to
+    three decimal places. The deltas are the distances from that mean to the
+    minimum (delta minus) and maximum (delta plus) energies.
 
     Parameters
     ----------
@@ -385,8 +390,9 @@ def add_energy_variables(
     """
     updated_ds = dataset.copy()
 
+    # Geometric mean of the minimum and maximum energy of each range
     energy_mean: np.ndarray = np.round(
-        np.mean(np.array([energy_min_values, energy_max_values]), axis=0), 3
+        np.sqrt(np.array(energy_min_values) * np.array(energy_max_values)), 3
     ).astype(np.float32)
 
     def _attrs(var_name: str) -> dict:
diff --git a/imap_processing/tests/hit/test_hit_utils.py b/imap_processing/tests/hit/test_hit_utils.py
index df948bace..a3d6647d4 100644
--- a/imap_processing/tests/hit/test_hit_utils.py
+++ b/imap_processing/tests/hit/test_hit_utils.py
@@ -259,7 +259,7 @@ def test_add_energy_variables():
     particle = "test_particle"
     energy_min = np.array([1.8, 4.0, 6.0], dtype=np.float32)
     energy_max = np.array([2.2, 6.0, 10.0], dtype=np.float32)
-    energy_mean = np.mean([energy_min, energy_max], axis=0)
+    energy_mean = np.array([1.99, 4.899, 7.746], dtype=np.float32)
 
     # Call the function
     dataset = add_energy_variables(dataset, particle, energy_min, energy_max)
@@ -268,15 +268,21 @@ def test_add_energy_variables():
     assert f"{particle}_energy_delta_minus" in dataset.data_vars
     assert f"{particle}_energy_delta_plus" in dataset.data_vars
     assert f"{particle}_energy_mean" in dataset.coords
-    assert np.all(
-        dataset[f"{particle}_energy_delta_minus"].values
-        == np.array(energy_mean - np.array(energy_min), dtype=np.float32)
+    np.testing.assert_allclose(
+        dataset[f"{particle}_energy_mean"].values, energy_mean, rtol=0, atol=1e-6
     )
-    assert np.all(
-        dataset[f"{particle}_energy_delta_plus"].values
-        == np.array(energy_max - energy_mean, dtype=np.float32)
+    np.testing.assert_allclose(
+        dataset[f"{particle}_energy_delta_minus"].values,
+        energy_mean - energy_min,
+        rtol=0,
+        atol=1e-6,
+    )
+    np.testing.assert_allclose(
+        dataset[f"{particle}_energy_delta_plus"].values,
+        energy_max - energy_mean,
+        rtol=0,
+        atol=1e-6,
     )
-    assert np.all(dataset[f"{particle}_energy_mean"].values == energy_mean)
 
 
 def test_sum_particle_data(sample_dataset):
@@ -364,9 +370,12 @@ def test_add_summed_particle_data_to_dataset(sample_dataset):
         len(energy_ranges),
     )
 
-    assert np.all(
-        dataset_to_update[f"{particle}_energy_mean"].values
-        == np.mean([[1.8, 4.0], [2.2, 6.0]], axis=0)
+    # Geometric means: sqrt(1.8 * 2.2) = 1.98997, sqrt(4.0 * 6.0) = 4.89898
+    np.testing.assert_allclose(
+        dataset_to_update[f"{particle}_energy_mean"].values,
+        np.array([1.99, 4.899], dtype=np.float32),
+        rtol=0,
+        atol=1e-6,
     )