[postprocessor] [plotter] add mass weighted radial average

plotter.py

@@ -1643,10 +1643,10 @@ class Plotter(Aggregator, BaseProcessor):
                 partial(
                     self._plot,
                     "/radial/radial_centers",
-                    "/radial/rad_avg_" + name,
+                    "/radial/rad_mwavg_" + name,
                 ),
-                "Azimuthal average of {}".format(name),
+                "Azimuthal mass weighted average of {}".format(name),
-                dependencies=["radial_centers", "rad_avg_" + name],
+                dependencies=["radial_centers", "rad_mwavg_" + name],
             )
 
             self.rules["avg_map_" + name] = PlotRule(
@@ -1656,6 +1656,13 @@ class Plotter(Aggregator, BaseProcessor):
                 dependencies=["avg_map_" + name],
             )
 
+            self.rules["mwavg_map_" + name] = PlotRule(
+                self,
+                partial(self._plot_map, "mwavg_map_" + name),
+                "Map of the mass weighted radial average of {}".format(name),
+                dependencies=["avg_map_" + name],
+            )
+
             self.rules["fluct_" + name] = PlotRule(
                 self,
                 partial(self._plot_map, "fluct_" + name, cmap="RdBu_r"),

postprocessor.py

@@ -876,7 +876,7 @@ class PostProcessor(HDF5Container):
         dmap_omega = rt_omega.process(self._cam[ax_los]).map.T
         return dmap_omega
 
-    def _toomreQ_disk(self, ax_los, omega_approx=False):
+    def _toomreQ_disk(self, ax_los, omega_approx=True):
         """
         Compute the Toomre Q parameter
         """
@@ -986,23 +986,35 @@ class PostProcessor(HDF5Container):
             bins = bins + (rr >= r).astype(int)
         return bins
 
-    def _rad_avg(self, name, ax_los="z"):
+    def _rad_avg(self, name, ax_los="z", mass_weighted=False):
-        radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read()
+        radial_bins = self.get_value("/radial/radial_bins_" + ax_los)
-        bins_on_map = self.save.get_node("/maps/bins_on_map_" + ax_los).read()
+        bins_on_map = self.get_value("/maps/bins_on_map_" + ax_los)
-        dmap = self.save.get_node("/maps/" + name + "_" + ax_los).read()
+        dmap = self.get_value("/maps/" + name + "_" + ax_los)
+
+        if mass_weighted:
+            coldens = self.get_value("/maps/coldens_" + ax_los)
 
         # mean of all the cells in the bin
         mean_bin = np.zeros(len(radial_bins) - 1)
         for j in range(len(radial_bins) - 1):
-            mean_bin[j] = np.mean(dmap[bins_on_map == j])
+            if mass_weighted:
+                weight = coldens[bins_on_map == j]
+                mean_bin[j] = np.mean(dmap[bins_on_map == j] * weight)
+                mean_bin[j] = mean_bin[j] / np.mean(weight)
+            else:
+                mean_bin[j] = np.mean(dmap[bins_on_map == j])
         return mean_bin
 
-    def _rad_avg_map(self, name, ax_los="z"):
+    def _rad_avg_map(self, name, ax_los="z", mass_weighted=False):
 
-        radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read()
+        radial_bins = self.get_value("/radial/radial_bins_" + ax_los)
-        bins_on_map = self.save.get_node("/maps/bins_on_map_" + ax_los).read()
+        bins_on_map = self.get_value("/maps/bins_on_map_" + ax_los)
-        rr = self.save.get_node("/maps/rr_" + ax_los).read()
+        rr = self.get_value("/maps/rr_" + ax_los)
-        mean_bin = self.save.get_node("/radial/rad_avg_" + name + "_" + ax_los).read()
+
+        if mass_weighted:
+            mean_bin = self.get_value("/radial/rad_mwavg_" + name + "_" + ax_los)
+        else:
+            mean_bin = self.get_value("/radial/rad_avg_" + name + "_" + ax_los)
 
         # Add value for border
         mean_bin = np.concatenate((mean_bin, [mean_bin[-1]]))
@@ -1026,16 +1038,22 @@ class PostProcessor(HDF5Container):
 
         return avg_map
 
-    def _fluct_map(self, name, ax_los="z"):
+    def _fluct_map(self, name, ax_los="z", mass_weighted=False):
 
         dmap = self.save.get_node("/maps/" + name + "_" + ax_los).read()
-        avg_map = self.save.get_node("/maps/avg_map_" + name + "_" + ax_los).read()
+        if mass_weighted:
-
+            avg_map = self.get_value("/maps/mwavg_map_" + name + "_" + ax_los)
+        else:
+            avg_map = self.get_value("/maps/avg_map_" + name + "_" + ax_los)
         return dmap / avg_map
 
-    def _rad_fluct_pdf(self, name, ax_los="z"):
+    def _rad_fluct_pdf(self, name, ax_los="z", mass_weighted=False):
-        fluct_map = self.save.get_node("/maps/fluct_" + name + "_" + ax_los).read()
+        if mass_weighted:
-        rr = self.save.get_node("/maps/rr_" + ax_los).read()
+            fluct_map = self.get_value("/maps/mwfluct_" + name + "_" + ax_los)
+        else:
+            fluct_map = self.get_value("/maps/fluct_" + name + "_" + ax_los)
+
+        rr = self.get_value("/maps/rr_" + ax_los)
 
         mask_pdf = (
             (rr > self.pp_params.disk.rmin_pdf)
@@ -1684,7 +1702,13 @@ class PostProcessor(HDF5Container):
                 "/radial",
                 dependencies=["radial_centers", "bins_on_map", name],
             )
-
+            self.rules["rad_mwavg_" + name] = Rule(
+                self,
+                partial(self._rad_avg, name, mass_weighted=True),
+                "Mass weighted azimuthal average of {}".format(name),
+                "/radial",
+                dependencies=["coldens", "radial_centers", "bins_on_map", name],
+            )
             self.rules["avg_map_" + name] = Rule(
                 self,
                 partial(self._rad_avg_map, name),
@@ -1692,6 +1716,13 @@ class PostProcessor(HDF5Container):
                 "/maps",
                 dependencies=["radial_centers", "bins_on_map", "rad_avg_" + name],
             )
+            self.rules["mwavg_map_" + name] = Rule(
+                self,
+                partial(self._rad_avg_map, name, mass_weighted=True),
+                "Interpolated map of azimuthal average of {}".format(name),
+                "/maps",
+                dependencies=["radial_centers", "bins_on_map", "rad_mwavg_" + name],
+            )
             self.rules["fluct_" + name] = Rule(
                 self,
                 partial(self._fluct_map, name),
@@ -1699,6 +1730,13 @@ class PostProcessor(HDF5Container):
                 "/maps",
                 dependencies=[name, "avg_map_" + name],
             )
+            self.rules["mwfluct_" + name] = Rule(
+                self,
+                partial(self._fluct_map, name, mass_weigthed=True),
+                "Fluctuation wrt to mass weighted average of {}".format(name),
+                "/maps",
+                dependencies=[name, "mwavg_map_" + name],
+            )
             self.rules["pdf_" + name] = Rule(
                 self,
                 partial(self._rad_fluct_pdf, name),
@@ -1706,7 +1744,15 @@ class PostProcessor(HDF5Container):
                 "/hist",
                 dependencies=["rr", "fluct_" + name],
             )
-
+            self.rules["mwpdf_" + name] = Rule(
+                self,
+                partial(self._rad_fluct_pdf, name, mass_weigthed=True),
+                "Probability density function of {} mass weighted fluctuations".format(
+                    name
+                ),
+                "/hist",
+                dependencies=["rr", "mwfluct_" + name],
+            )
             self.rules["fit_pdf_" + name] = Rule(
                 self,
                 partial(self._fit_pdf, name),