2342 hi l2 update how statistical uncertainty is combined from calibration products (#2364)

* Change equation for how we combine statistical uncertainties across calibration products

* Fix incorrect calculation of stat_unc_sc in interpolate_map_flux_to_helio_frame

* Add check for interpolated statistical uncertianty

Author: tmplummer

imap_processing/ena_maps/utils/corrections.py

@@ -707,10 +707,11 @@ def interpolate_map_flux_to_helio_frame(
         )
 
         # Statistical uncertainty propagation (Equation 75):
-        # δJ = J * sqrt((δJ_left/J_left)^2 * (1 + unc_factor^2) + (δJ_right/J_right)^2)
+        # δJ = J * sqrt((δJ_left/J_left)^2 * (1 + unc_factor^2)
+        #               + unc_factor^2 * (δJ_right/J_right)^2)
         stat_unc_sc = flux_sc * np.sqrt(
             (stat_unc_left / flux_left) ** 2 * (1.0 + unc_factor**2)
-            + (stat_unc_right / flux_right) ** 2
+            + unc_factor**2 * (stat_unc_right / flux_right) ** 2
         )
 
         # Systematic uncertainty propagation (Equation 76):

imap_processing/hi/hi_l2.py

@@ -379,8 +379,9 @@ def combine_calibration_products(
         combined_flux = weighted_flux_sum / flux_weights.sum(dim="calibration_prod")
 
     map_ds["ena_intensity"] = combined_flux
+    # Statistical uncertainty
     map_ds["ena_intensity_stat_uncert"] = np.sqrt(
-        (map_ds["ena_intensity_stat_uncert"] ** 2).sum(dim="calibration_prod")
+        1 / (1 / (map_ds["ena_intensity_stat_uncert"] ** 2)).sum(dim="calibration_prod")
     )
     # For systematic error, just do quadrature sum over the systematic error for
     # each calibration product.

imap_processing/tests/ena_maps/test_corrections.py

@@ -712,6 +712,10 @@ def test_power_law_interpolation_accuracy(self):
         expected_flux_middle = (
             (e_sc**power_law_slope) * (e_helio / e_sc) * spatial_factor
         )
+        unc_factor = np.log(e_sc / 500) / np.log(1000 / 500)
+        expected_stat_uncert = expected_flux_middle * np.sqrt(
+            0.1**2 * (1 + unc_factor**2) + unc_factor**2 * 0.1**2
+        )
 
         # Compare interpolated result to expected value
         # (should be very close for a perfect power-law)
@@ -721,6 +725,13 @@ def test_power_law_interpolation_accuracy(self):
             rtol=1e-10,
         )
 
+        # Check expected stat. unc.
+        np.testing.assert_allclose(
+            result_ds["ena_intensity_stat_uncert"].values[1, 0],
+            expected_stat_uncert,
+            rtol=1e-10,
+        )
+
         # The flux should be finite and positive
         assert np.all(np.isfinite(result_ds["ena_intensity"].values))
         assert np.all(result_ds["ena_intensity"].values > 0)

imap_processing/tests/hi/test_hi_l2.py

@@ -632,9 +632,8 @@ def test_statistical_uncertainty_combination_correctness():
     result_ds = combine_calibration_products(test_ds, geom_factors, esa_energies)
 
     # Manual calculation of expected statistical uncertainty combination
-    # stat_weights = 1/σ² = [1/151, 1/151]
-    # combined_stat_unc = sqrt(1/sum(stat_weights)) = sqrt(2/151) = sqrt(20) ≈ 4.47
-    expected_combined_stat_unc = np.sqrt(np.sum(stat_unc_values**2))
+    # combined_stat_unc = sqrt(1/sum(1 / stat_unc**2))
+    expected_combined_stat_unc = np.sqrt(1 / np.sum(1 / stat_unc_values**2))
     flux_weights = 1.0 / (np.array([101, 101]) + np.array([4, 16]))
     expected_flux = np.sum(flux_values.squeeze() * flux_weights) / np.sum(flux_weights)