Magnetic Intensity map, streamlines and averaged field versus density

plotter.py

@@ -987,6 +987,18 @@ class Plotter(Aggregator, BaseProcessor):
                 "Density slice with magnetic field overlay",
                 dependencies=["rho", "B_h", "B_v"],
             ),
+            "rho_B_vel": PlotRule(
+                self,
+                partial(
+                    self._plot_map,
+                    "rho",
+                    label=r"$\rho$",
+                    unit=cst.Msun_pc3,
+                    overlays=[self._overlay_B, self._overlay_speed],
+                ),
+                "Density slice with magnetic field and velocity  overlay",
+                dependencies=["rho", "B_h", "B_v", "speed_h", "speed_v"],
+            ),
             "jeans_ratio": PlotRule(
                 self,
                 partial(
@@ -1019,6 +1031,12 @@ class Plotter(Aggregator, BaseProcessor):
                 "$\rho$-PDF",
                 dependencies=["rho_pdf"],
             ),
+            "cos_pdf": PlotRule(
+                self,
+                partial(self._plot_hist, "cos_pdf"),
+                "cos-PDF",
+                dependencies=["cos_pdf"],
+            ),
             "avg_coldens_pdf": PlotRule(
                 self,
                 partial(
@@ -1044,6 +1062,20 @@ class Plotter(Aggregator, BaseProcessor):
                 "P-PDF on a 2D slice ",
                 dependencies=["P_pdf"],
             ),
+            "B_int": PlotRule(
+                self,
+                partial(
+                    self._plot_map, "B_int", label=r"$\mid \mathrm{B} \mid$", unit=cst.T
+                ),
+                "Magnetic intensity map",
+                dependencies=["B_int"],
+            ),
+            "Brho": PlotRule(
+                self,
+                partial(self._plot, "/datasets/Brho/rho", "/datasets/Brho/B"),
+                "Brho on a 2D slice ",
+                dependencies=["Brho"],
+            ),
             "rho_prof": PlotRule(
                 self,
                 partial(self._plot, "/profile/axis", "/profile/rho_prof"),

postprocessor.py

@@ -8,6 +8,23 @@ vol_func = lambda dset: dset["dx"] ** 3  # Volume function
 getter_T = lambda dset: dset["P"] / dset["rho"]  # Temperature
 
 
+def _getter_abs_cos_vB(dset):
+    B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1))
+    v_norm = np.sqrt(np.sum(dset["vel"] ** 2, axis=1))
+    # Compute the dot product in each cell
+    dot_prod = np.einsum("ij,ij->i", dset["vel"], dset["Br"])
+    return np.abs(dot_prod) / (v_norm * B_norm)
+
+
+def _getter_B_int(dset):
+    B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1))
+    return B_norm
+
+
+def _getter_rho(dset):
+    return dset["rho"]
+
+
 class PostProcessor(HDF5Container):
     """
     This class enable to compute and save derived quantities from the raw output
@@ -285,6 +302,28 @@ class PostProcessor(HDF5Container):
         centers = 0.5 * (edges[1:] + edges[:-1])
         return (np.stack([values, centers]), {"logbins": logbins})
 
+    def _Brho(self, bins=100):
+        self.load_cells()
+        B = _getter_B_int(self.cells)
+        rho = _getter_rho(self.cells)
+        B = np.log10(B)
+        rho = np.log10(rho)
+        abscisse = np.linspace(np.min(rho), np.max(rho), bins)
+        weights = mass_func(self.cells)
+        bin_number = np.zeros(len(B))
+        for rho_min in abscisse[:-1]:
+            bin_number = bin_number + (rho > rho_min).astype(int)
+
+        B_mean = np.zeros(len(abscisse) - 1)
+        for i in range(len(B_mean)):
+            B_mean[i] = np.mean(B[bin_number == i])
+        values, edges = np.histogram(B, abscisse, weights=weights)
+        centers = 0.5 * (edges[1:] + edges[:-1])
+        print("centers")
+        if self.pp_params.process.unload_cells:
+            self.unload_cells()
+        return {"rho": centers, "B": B_mean}
+
     def _mwa_sigma(self, axes=["x", "y", "z"]):
         mw_speed = self.save.get_node("/globals/mwa_speed").read()
 
@@ -341,6 +380,15 @@ class PostProcessor(HDF5Container):
     def _B_v(self, ax_los, z=0.0):
         return self._vector_v("Br", "unit_mag", ax_los, z)
 
+    def _B_int(self, ax_los, z=0.0):
+        B_op = ScalarOperator(
+            lambda dset: np.sqrt(np.sum(dset["Br"] ** 2, axis=1)),
+            self._ro.info["unit_mag"],
+        )
+        dmap_B = (slicing.SliceMap(self._amr, self._cam[ax_los], B_op, z=z)).map.T
+        dmap_rho = self.save.get_node("/maps/rho_{}".format(ax_los)).read()
+        return dmap_B
+
     def _temperature(self, ax_los, z=0.0):
         P_op = ScalarOperator(lambda dset: dset["P"], self._ro.info["unit_pressure"])
         dmap_P = (slicing.SliceMap(self._amr, self._cam[ax_los], P_op, z=z)).map.T
@@ -749,6 +797,14 @@ class PostProcessor(HDF5Container):
                 "/maps",
                 dependencies=["radial_bins", "rr"],
             ),
+            "B_int": Rule(
+                self,
+                self._B_int,
+                "Magnetic intensity slice",
+                "/maps",
+                dependencies=["rho"],
+                unit=self.info["unit_mag"],
+            ),
             # PDF
             "rho_pdf": Rule(
                 self,
@@ -771,6 +827,20 @@ class PostProcessor(HDF5Container):
                 "/hist",
                 unit=self.info["unit_pressure"],
             ),
+            "cos_pdf": Rule(
+                self,
+                partial(self._vol_pdf, _getter_abs_cos_vB, logbins=False),
+                "Global cos-PDF",
+                "/hist",
+                unit=cst.none,
+            ),
+            "Brho": Rule(
+                self,
+                partial(self._Brho),
+                "Brho average",
+                "/datasets",
+                unit=self.info["unit_mag"],
+            ),
             # Profiles
             "axis": Rule(
                 self,