Harrison 2381 - PMCL/TLAY - Use average of last three spacecraft potentials and halo core values when finding breackpoints. Instead of 4 values except most recent.

Author: Menlo Innovations - CAVA Project

imap_processing/swe/l3/science/pitch_calculations.py

@@ -20,9 +20,8 @@ def piece_wise_model(x: np.ndarray, b0: float, b1: float,
                                ]))
 
 
-def find_breakpoints(energies: np.ndarray, averaged_psd: np.ndarray, initial_spacecraft_potential_guess: float,
-                     initial_core_halo_break_point_guess: float,
-                     latest_spacecraft_potential: float, latest_core_halo_breakpoint: float,
+def find_breakpoints(energies: np.ndarray, averaged_psd: np.ndarray, latest_spacecraft_potentials: list[float],
+                     latest_core_halo_break_points: list[float],
                      config: SweConfiguration) -> tuple[
     float, float]:
     log_psd = np.log(averaged_psd)
@@ -38,7 +37,8 @@ def find_breakpoints(energies: np.ndarray, averaged_psd: np.ndarray, initial_spa
     b3: float = slopes[core_index]
     b5: float = slopes[halo_index]
     b0: float = np.exp(log_psd[0] + b1 * energies[0])
-    initial_guesses = (b0, b1, initial_spacecraft_potential_guess, b3, initial_core_halo_break_point_guess, b5)
+    initial_guesses = (
+        b0, b1, np.average(latest_spacecraft_potentials), b3, np.average(latest_core_halo_break_points), b5)
 
     first_min_index = 0
     for i in range(1, len(slope_ratios) - 1):
@@ -59,8 +59,8 @@ def find_breakpoints(energies: np.ndarray, averaged_psd: np.ndarray, initial_spa
         delta_b2 = -1.0
         delta_b4 = 10
 
-    return try_curve_fit_until_valid(energies, log_psd, initial_guesses, latest_spacecraft_potential,
-                                     latest_core_halo_breakpoint, delta_b2, delta_b4)
+    return try_curve_fit_until_valid(energies, log_psd, initial_guesses, latest_spacecraft_potentials[-1],
+                                     latest_core_halo_break_points[-1], delta_b2, delta_b4)
 
 
 def try_curve_fit_until_valid(energies: np.ndarray, log_psd: np.ndarray, initial_guesses: tuple[float, ...],

imap_processing/swe/swe_processor.py

@@ -50,16 +50,15 @@ def calculate_pitch_angle_products(self, dependencies: SweL3Dependencies) -> Swe
         energy_spectrum_inbound = []
         energy_spectrum_outbound = []
         spacecraft_potential_history = [config["spacecraft_potential_initial_guess"] for _ in
-                                        range(4)]
-        halo_core_history = [config["core_halo_breakpoint_initial_guess"] for _ in range(4)]
+                                        range(3)]
+        halo_core_history = [config["core_halo_breakpoint_initial_guess"] for _ in range(3)]
 
         for i in range(len(swe_epoch)):
             averaged_psd = average_over_look_directions(swe_l2_data.phase_space_density[i],
                                                         np.array(config["geometric_fractions"]))
             spacecraft_potential, halo_core = find_breakpoints(swe_l2_data.energy, averaged_psd,
-                                                               np.average(spacecraft_potential_history[:3]),
-                                                               np.average(halo_core_history[:3]),
-                                                               spacecraft_potential_history[-1], halo_core_history[-1],
+                                                               spacecraft_potential_history,
+                                                               halo_core_history,
                                                                config)
             spacecraft_potential_history = [*spacecraft_potential_history[1:], spacecraft_potential]
             halo_core_history = [*halo_core_history[1:], halo_core]

tests/swe/l3/science/test_pitch_calculations.py

@@ -122,8 +122,7 @@ def test_find_breakpoints_with_noisy_data(self):
                              3.69484446e-01, 2.19359553e-01, 1.19059738e-01, 5.64115725e-02,
                              2.30604686e-02, 9.14406238e-03, 4.24754874e-03, 1.61814681e-03])
         spacecraft_potential, core_halo_breakpoint = find_breakpoints(
-            xs, avg_flux, 10, 80,
-            11, 81, config)
+            xs, avg_flux, [10, 10, 10], [80, 80, 80], config)
         self.assertAlmostEqual(11.1, spacecraft_potential, 1)
         self.assertAlmostEqual(81.1, core_halo_breakpoint, 1)
 
@@ -146,18 +145,18 @@ def test_find_breakpoints_with_synthetic_data(self):
                 noise_floor = 1
                 avg_flux += noise_floor
                 spacecraft_potential, core_halo_breakpoint = find_breakpoints(
-                    xs, avg_flux, 10, 80,
-                    11, 82, config)
+                    xs, avg_flux, [10, 10, 10], [80, 80, 80],
+                    config)
                 self.assertAlmostEqual(expected_potential, spacecraft_potential, 2)
                 self.assertAlmostEqual(expected_core_halo, core_halo_breakpoint, 0)
 
     def test_find_breakpoints_using_initial_guess(self):
         config = build_swe_configuration()
 
         cases = [
-            (4, 40, 4, 50),
-            (10, 80, 12, 100),
-            (12, 60, 10, 80),
+            (4, 40, [4, 4, 4], [50, 50, 50]),
+            (10, 80, [12, 12, 12], [100, 100, 100]),
+            (12, 60, [10, 10, 10], [80, 80, 80]),
         ]
         for case in cases:
             with self.subTest(case):
@@ -172,7 +171,7 @@ def test_find_breakpoints_using_initial_guess(self):
                 avg_flux += noise_floor
 
                 spacecraft_potential, core_halo_breakpoint = find_breakpoints(
-                    xs, avg_flux, guess_potential, guess_halo, 10, 80, config)
+                    xs, avg_flux, guess_potential, guess_halo, config)
                 self.assertAlmostEqual(expected_potential, spacecraft_potential, 2)
                 self.assertAlmostEqual(expected_core_halo, core_halo_breakpoint, 0)
 
@@ -199,9 +198,8 @@ def test_find_breakpoints_determines_b_deltas_correctly(self, mock_try_curve_fit
                 avg_flux = np.exp(log_flux)
 
                 result = find_breakpoints(
-                    xs, avg_flux, 10, 80, 15,
-                    90, config)
-                mock_try_curve_fit_until_valid.assert_called_with(ANY, ANY, ANY, 15, 90, expected_b2_delta,
+                    xs, avg_flux, [10, 10, 10], [80, 80, 80], config)
+                mock_try_curve_fit_until_valid.assert_called_with(ANY, ANY, ANY, 10, 80, expected_b2_delta,
                                                                   expected_b4_delta)
 
                 self.assertEqual(mock_try_curve_fit_until_valid.return_value, result)
@@ -233,8 +231,8 @@ def test_find_breakpoints_uses_config_for_slope_guesses(self, mock_curve_fit):
                 avg_flux = np.exp(log_flux)
 
                 spacecraft_potential, core_halo_breakpoint = find_breakpoints(
-                    xs, avg_flux, 10, 80,
-                    11, 81, config)
+                    xs, avg_flux, [10, 10, 10], [80, 80, 80],
+                    config)
                 expected_guesses = [ANY, b1, 10, b3, 80, b5]
                 rounded_actuals = [round(x, 6) for x in mock_curve_fit.call_args.args[3]]
                 self.assertEqual(expected_guesses, rounded_actuals)
@@ -264,8 +262,7 @@ def test_find_breakpoints_uses_config_for_slope_ratio(self, mock_curve_fit):
                 avg_flux = np.exp(log_flux)
 
                 spacecraft_potential, core_halo_breakpoint = find_breakpoints(
-                    xs, avg_flux, 10, 80,
-                    11, 81, config)
+                    xs, avg_flux, [10, 10, 10], [80, 80, 80], config)
 
                 np.testing.assert_almost_equal(mock_curve_fit.call_args.args[1], xs[:data_length])
                 np.testing.assert_almost_equal(mock_curve_fit.call_args.args[2], log_flux[:data_length])

tests/swe/test_swe_processor.py

@@ -146,14 +146,13 @@ def test_calculate_pitch_angle_products(self, mock_find_closest_neighbor,
             call(NumpyArrayMatcher(swe_l2_data.phase_space_density[1]), NumpyArrayMatcher(geometric_fractions)),
             call(NumpyArrayMatcher(swe_l2_data.phase_space_density[2]), NumpyArrayMatcher(geometric_fractions))])
 
-        spacecraft_potential_initial_guess = swe_config['spacecraft_potential_initial_guess']
-        halo_core_initial_guess = swe_config['core_halo_breakpoint_initial_guess']
         mock_find_breakpoints.assert_has_calls([
-            call(swe_l2_data.energy, mock_average_over_look_directions.return_value, spacecraft_potential_initial_guess,
-                 halo_core_initial_guess, 15, 90, swe_config),
-            call(swe_l2_data.energy, mock_average_over_look_directions.return_value, spacecraft_potential_initial_guess,
-                 halo_core_initial_guess, 12, 96, swe_config),
-            call(swe_l2_data.energy, mock_average_over_look_directions.return_value, 14, 92, 16, 86, swe_config),
+            call(swe_l2_data.energy, mock_average_over_look_directions.return_value, [15, 15, 15],
+                 [90, 90, 90], swe_config),
+            call(swe_l2_data.energy, mock_average_over_look_directions.return_value, [15, 15, 12],
+                 [90, 90, 96], swe_config),
+            call(swe_l2_data.energy, mock_average_over_look_directions.return_value, [15, 12, 16],
+                 [90, 96, 86], swe_config),
         ])
 
         self.assertEqual(UpstreamDataDependency("swe", "l3", datetime(2025, 2, 21),
@@ -174,7 +173,9 @@ def test_calculate_pitch_angle_products(self, mock_find_closest_neighbor,
         def call_with_array_matchers(*args):
             return call(*[NumpyArrayMatcher(x) for x in args])
 
-        self.assertEqual(mock_correct_and_rebin.call_args_list, [
+        actual_calls = mock_correct_and_rebin.call_args_list
+
+        expected_calls = [
             call_with_array_matchers(swe_l2_data.flux[0], swe_l2_data.energy - 12, swe_l2_data.inst_el,
                                      swe_l2_data.inst_az_spin_sector[0],
                                      closest_mag_data[0], closest_swapi_data[0], swe_config),
@@ -193,7 +194,8 @@ def call_with_array_matchers(*args):
             call_with_array_matchers(swe_l2_data.phase_space_density[2], swe_l2_data.energy - 19, swe_l2_data.inst_el,
                                      swe_l2_data.inst_az_spin_sector[2],
                                      closest_mag_data[2], closest_swapi_data[2], swe_config)
-        ])
+        ]
+        self.assertEqual(actual_calls, expected_calls)
         mock_integrate_distribution_to_get_1d_spectrum.assert_has_calls([
             call(rebinned_by_pitch[1], swe_config),
             call(rebinned_by_pitch[3], swe_config),