Harrison 2382 - SWYN/PSCH - Now performing the main fitting logic on halo electrons, regress, find_density_temp

Author: Menlo Innovations - CAVA Project

imap_processing/swe/swe_processor.py

@@ -58,33 +58,60 @@ def calculate_moment_products(self, dependencies: SweL3Dependencies):
 
             weights: np.ndarray[float] = compute_maxwellian_weight_factors(ccounts)
 
-            halo_core_breakpoint_index: int = next(
-                i - 1 for i, energy in enumerate(swe_l2_data.energy) if energy > halo_core)
             spacecraft_potential_core_breakpoint_index: int = next(
                 i for i, energy in enumerate(swe_l2_data.energy) if energy >= spacecraft_potential)
+            halo_core_breakpoint_index: int = next(
+                i - 1 for i, energy in enumerate(swe_l2_data.energy) if energy > halo_core)
+
+            core_end_index = halo_core_breakpoint_index
+            while True:
+                filtered_velocity_vectors, filtered_weights, filtered_yreg = filter_and_flatten_regress_parameters(
+                    corrected_energy_bins,
+                    velocity_vectors,
+                    swe_l2_data.phase_space_density[i],
+                    weights,
+                    spacecraft_potential_core_breakpoint_index, core_end_index)
+
+                fit_function, chisq = regress(filtered_velocity_vectors,
+                                              filtered_weights, filtered_yreg)
+                core_moments = calculate_fit_temperature_density_velocity(fit_function)
+
+                if 0 < core_moments.density < 185 or (
+                        core_end_index - spacecraft_potential_core_breakpoint_index) <= 3:
+                    break
+                else:
+                    core_end_index -= 1
+
+            rtn_velocity = rotate_dps_vector_to_rtn(swe_epoch[i],
+                                                    np.array(
+                                                        [core_moments.velocity_x, core_moments.velocity_y,
+                                                         core_moments.velocity_z]))
+            rotate_temperature(swe_epoch[i], core_moments.alpha, core_moments.beta)
 
+            halo_end_index = len(swe_l2_data.energy)
             while True:
                 filtered_velocity_vectors, filtered_weights, filtered_yreg = filter_and_flatten_regress_parameters(
                     corrected_energy_bins,
                     velocity_vectors,
                     swe_l2_data.phase_space_density[i],
                     weights,
-                    spacecraft_potential_core_breakpoint_index, halo_core_breakpoint_index)
+                    halo_core_breakpoint_index, halo_end_index)
 
                 fit_function, chisq = regress(filtered_velocity_vectors,
                                               filtered_weights, filtered_yreg)
-                moments = calculate_fit_temperature_density_velocity(fit_function)
+                halo_moments = calculate_fit_temperature_density_velocity(fit_function)
 
-                if 0 < moments.density < 185 or (
-                        halo_core_breakpoint_index - spacecraft_potential_core_breakpoint_index) <= 3:
+                if 0 < halo_moments.density < 185 or (
+                        halo_end_index - halo_core_breakpoint_index) <= 3:
                     break
                 else:
-                    halo_core_breakpoint_index -= 1
+                    halo_end_index -= 1
 
             rtn_velocity = rotate_dps_vector_to_rtn(swe_epoch[i],
                                                     np.array(
-                                                        [moments.velocity_x, moments.velocity_y, moments.velocity_z]))
-            rotate_temperature(swe_epoch[i], moments.alpha, moments.beta)
+                                                        [halo_moments.velocity_x, halo_moments.velocity_y,
+                                                         halo_moments.velocity_z]))
+            rotate_temperature(swe_epoch[i], halo_moments.alpha, halo_moments.beta)
 
     def calculate_pitch_angle_products(self, dependencies: SweL3Dependencies) -> SweL3Data:
         swe_l2_data = dependencies.swe_l2_data

run_local.py

@@ -298,9 +298,9 @@ def create_hit_direct_event_cdf():
 
     if "swe_moments" in sys.argv:
         dependencies = SweL3Dependencies.from_file_paths(
-            get_test_data_path("swe/imap_swe_l2_sci-with-ace-data_19990609_v002.cdf"),
-            get_test_data_path("swe/imap_mag_l1d_mago-normal_19990609_v001.cdf"),
-            get_test_data_path("swe/imap_swapi_l3a_proton-sw_19990609_v001.cdf"),
+            get_test_data_path("swe/imap_swe_l2_sci-with-ace-data_19990906_v002.cdf"),
+            get_test_data_path("swe/imap_mag_l1d_mago-normal_19990906_v001.cdf"),
+            get_test_data_path("swe/imap_swapi_l3a_proton-sw_19990906_v001.cdf"),
             get_test_data_path("swe/example_swe_config.json"),
         )
         print(create_swe_moments_cdf(dependencies))