Averaged magnetic field versus time

comparator.py

@@ -574,6 +574,7 @@ class Comparator(Aggregator, HDF5Container):
 
         self._gen_rule_time_global("mwa_sigma", "time_sigma", unit="unit_velocity")
         self._gen_rule_time_global("max_fluct_coldens")
+        self._gen_rule_time_global("mwa_B_int", unit="unit_mag")
 
         for name in [
             "issfr",

plotter.py

@@ -78,12 +78,14 @@ class Plotter(Aggregator, BaseProcessor):
         "coldens0": "$\Sigma_0$",
         "sfr_avg_window": "window",
         "comp_frac": "$\\zeta$",
+        "bx_bound": "$B_x$",
     }
 
     # Conversion table from namelist values (from amses config file)  into LaTex strings
     value_convert = {
         "sfr_avg_window": lambda x: "${:g}$ Myr".format(80 * x),
         "comp_frac": lambda x: "${:g}$".format(1 - x),
+        "bx_bound": lambda x: "${:g}$".format(5.267 * 10 ** (-10) * x),
     }
 
     def __init__(
@@ -1066,7 +1068,7 @@ class Plotter(Aggregator, BaseProcessor):
                     unit=cst.T,
                     overlays=[self._overlay_B, self._overlay_speed],
                 ),
-                "Density slice with magnetic field and velocity  overlay",
+                "Magnetic slice with magnetic field and velocity  overlay",
                 dependencies=["B_int", "B_h", "B_v", "speed_h", "speed_v"],
             ),
             "jeans_ratio": PlotRule(
@@ -1142,10 +1144,28 @@ class Plotter(Aggregator, BaseProcessor):
             ),
             "Brho": PlotRule(
                 self,
-                partial(self._plot, "/datasets/Brho/rho", "/datasets/Brho/B"),
+                partial(
+                    self._plot,
+                    "/datasets/Brho/rho",
+                    "/datasets/Brho/B",
+                    label=r"$\mathrm{B} $",
+                    put_time=True,
+                ),
                 "Brho on a 2D slice ",
                 dependencies=["Brho"],
             ),
+            "Ek_Eb_rho": PlotRule(
+                self,
+                partial(
+                    self._plot,
+                    "/datasets/Ek_Eb_rho/rho",
+                    "/datasets/Ek_Eb_rho/Ek_Eb_rho",
+                    label=r"Ek/Eb",
+                    put_time=True,
+                ),
+                "Ek/Eb on a 2D slice ",
+                dependencies=["Ek_Eb_rho", "mwa_speed"],
+            ),
             "rho_prof": PlotRule(
                 self,
                 partial(self._plot, "/profile/axis", "/profile/rho_prof"),
@@ -1256,6 +1276,15 @@ class Plotter(Aggregator, BaseProcessor):
                 kind="series",
                 dependencies=["time_sigma"],
             ),
+            "mwa_B_int": PlotRule(
+                self,
+                partial(
+                    self._plot, "/series/time", "/series/time_mwa_B_int", xunit=cst.Myr
+                ),
+                "Magnetic intensity average",
+                kind="series",
+                dependencies=["time_mwa_B_int"],
+            ),
             "max_fluct_coldens": PlotRule(
                 self,
                 partial(

postprocessor.py

@@ -17,6 +17,11 @@ def getter_rho(dset):
     return dset["rho"]
 
 
+def getter_v_norm(dset):
+    v_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1))
+    return v_norm
+
+
 class PostProcessor(HDF5Container):
     """
     This class enable to compute and save derived quantities from the raw output
@@ -333,6 +338,52 @@ class PostProcessor(HDF5Container):
             self.unload_cells()
         return ({"rho": centers, "B": B_mean}, {"logbins": logbins})
 
+    def _Ek_Eb_rho(self, bins=100, logbins=True):
+        """
+        Mean of Ek/Eb in rho bins
+        """
+        self.load_cells()
+        mean_speed = self.save.get_node("/globals/mwa_speed").read()
+        vel_fluct = (self.cells)["vel"] - mean_speed
+        B_norm = getter_B_int(self.cells)
+        v_norm = np.sqrt(np.sum((vel_fluct * 10 ** (3)) ** 2, axis=1))  # v [km/s]
+        print(v_norm)
+        rho = getter_rho(self.cells)
+        eb = 0.5 * (B_norm) ** 2 / (4 * np.pi * 10 ** (-7))  # mettre le bon mu
+        ek = 0.5 * v_norm ** 2 * rho  # mettre la masse de la cellule
+        rapport = ek / eb
+
+        if logbins:
+            rho_bins = np.logspace(
+                np.log10(np.min(rho)), np.log10(np.max(rho)), bins, base=10
+            )
+        else:
+            rho_bins = np.linspace(np.min(rho), np.max(rho), bins)
+
+        weights = mass_func(self.cells)
+
+        # For each cell, bin_number contains the number of the bins it belongs to
+        bin_number = np.zeros(len(B_norm))
+
+        # Go through the min value of rho of each bin
+        for rho_min in rho_bins[:-1]:
+            bin_number = bin_number + (rho > rho_min).astype(int)
+
+        # Compute the mean in each bin
+        ek_eb = np.zeros(len(rho_bins) - 1)
+        for i in range(len(ek_eb)):
+            ek_eb[i] = np.mean(rapport[bin_number == i])
+
+        # Get the center of each bin
+        if logbins:
+            centers = 10 ** (0.5 * (np.log10(rho_bins[1:]) + np.log10(rho_bins[:-1])))
+        else:
+            centers = 0.5 * (rho_bins[1:] + rho_bins[:-1])
+
+        if self.pp_params.process.unload_cells:
+            self.unload_cells()
+        return ({"rho": centers, "Ek_Eb_rho": ek_eb}, {"logbins": logbins})
+
     def cos_vfluct_B(self):
 
         mean_speed = self.save.get_node("/globals/mwa_speed").read()
@@ -836,6 +887,14 @@ class PostProcessor(HDF5Container):
                 "/datasets",
                 unit={"rho": self.info["unit_density"], "B": self.info["unit_mag"]},
             ),
+            "Ek_Eb_rho": Rule(
+                self,
+                self._Ek_Eb_rho,
+                "Average of Ek/Eb as a function of rho",
+                "/datasets",
+                dependencies=["mwa_speed"],
+                unit={"rho": self.info["unit_density"], "Ek_Eb_rho": cst.none},
+            ),
             # Profiles
             "axis": Rule(
                 self,
@@ -875,6 +934,13 @@ class PostProcessor(HDF5Container):
                 dependencies={"mwa_speed": None},
                 unit=self.info["unit_velocity"],
             ),
+            "mwa_B_int": Rule(
+                self,
+                partial(self._vol_avg, getter_B_int),
+                "Mass weighted Magnetic intensity average",
+                "/globals",
+                unit=self.info["unit_mag"],
+            ),
         }
 
         # Average and other