Added sbeta_onavg + Switched to densty PDF

comparator.py

@@ -367,6 +367,29 @@ class Comparator(Aggregator, HDF5Container):
             turb_power[run] = energy / dt
         return turb_power
 
+    def _sbeta_onavg(self):
+        """
+        [ismfeed] Compute the slope of the Sigma pdf as a function of the value of beta
+        """
+        col_pdf = self.get_value("/comp/avg_time_coldens_pdf_z")
+        beta = self.get_value("/comp/nml_cloud_params/beta_cool")
+
+        slope = np.zeros(len(col_pdf["runs"]))
+        origin = np.zeros(len(col_pdf["runs"]))
+        stderr = np.zeros(len(col_pdf["runs"]))
+
+        for i, run in enumerate(col_pdf["runs"]):
+            values, centers = col_pdf["mean"][i]
+            mask_fit = (
+                (centers > self.pp_params.pdf.xmin_fit)
+                & (centers < self.pp_params.pdf.xmax_fit)
+                & (values > np.max(values) * self.pp_params.pdf.fit_cut)
+            )
+            (slope[i], origin[i], correlation, _, stderr[i]) = linregress(
+                centers[mask_fit], np.log10(values[mask_fit])
+            )
+        return {"beta": beta, "slope": slope, "origin": origin, "stderr": stderr}
+
     def _gen_rule_time_global(
         self,
         glob_name,
@@ -584,6 +607,15 @@ class Comparator(Aggregator, HDF5Container):
                 group="/series",
                 dependencies={"time": None, "fit_pdf_coldens": "z"},
             ),
+            "sbeta_onavg": Rule(
+                self,
+                partial(self._sbeta_onavg),
+                group="/comp",
+                dependencies={
+                    "avg_time_coldens_pdf": "z",
+                    "nml": "cloud_params/beta_cool",
+                },
+            ),
             # namelist
             "nml": Rule(
                 self,

plotter.py

@@ -708,8 +708,6 @@ class Plotter(Aggregator, BaseProcessor):
 
         P.streamplot(hh, vv, map_Bh_red, map_Bv_red, **kwargs)
 
-        P.streamplot(hh, vv, map_Bh_red, map_Bv_red)
-
     def _plot_radial(
         self,
         name,
@@ -719,6 +717,7 @@ class Plotter(Aggregator, BaseProcessor):
         xlog=False,
         ylog=False,
         ytransform=None,
+        label=None,
         title=None,
         nml_key=None,
         put_title=True,
@@ -763,6 +762,10 @@ class Plotter(Aggregator, BaseProcessor):
 
             P.title(title)
 
+        if label is None:
+            label = title
+        P.plot(bin_centers, mean_bin, label=label, **kwargs)
+
         P.plot(bin_centers, mean_bin, label=title, **kwargs)
 
     def _plot_hist(
@@ -863,7 +866,7 @@ class Plotter(Aggregator, BaseProcessor):
         if not ylabel is None:
             P.ylabel(ylabel)
 
-        if not ax_los is None and "/hist/fit_" + name + "_" + ax_los in self.save:
+        if ax_los is not None and "/hist/fit_" + name + "_" + ax_los in self.save:
             slope = node.attrs.slope
             origin = node.attrs.origin
             P.plot(
@@ -900,6 +903,7 @@ class Plotter(Aggregator, BaseProcessor):
         nml_key=None,
         run=None,
         runs=None,
+        yerr=None,
         yerr_kind="std",
         sigma_err=2.0,
         grid=False,
@@ -955,7 +959,6 @@ class Plotter(Aggregator, BaseProcessor):
             else:
                 label = time_str
 
-        yerr = None
         if node_y._v_attrs.CLASS == "ARRAY":
             x = node_x.read() * xunit_old.express(xunit)
             y = node_y.read() * yunit_old.express(yunit)
@@ -981,7 +984,14 @@ class Plotter(Aggregator, BaseProcessor):
                         color = colors[nml]
                     except:
                         color = colors(nml)
-                (base_line,) = P.plot(x, y, label=label, color=color, **kwargs)
+                if yerr is None:
+                    (base_line,) = P.plot(x, y, label=label, color=color, **kwargs)
+                else:
+                    if isinstance(yerr, str):
+                        yerr = self.save.get_node(yerr).read()
+
+                    base_line, _, _ = P.errorbar(x, y, yerr=yerr, label=label, **kwargs)
+
         elif "mean" in node_y:
             x = node_x.read() * xunit_old.express(xunit)
             y = node_y.mean.read() * yunit_old.express(yunit)
@@ -1401,7 +1411,7 @@ class Plotter(Aggregator, BaseProcessor):
                 dependencies=["axis", "rho_prof"],
             ),
             "pspec": PlotRule(self, self._pspec, dependencies={"pspec": None}),
-            "kappa_beta": PlotRule(
+            "sbeta": PlotRule(
                 self,
                 partial(
                     self._plot,
@@ -1415,6 +1425,18 @@ class Plotter(Aggregator, BaseProcessor):
                     "avg_time_pdf_slope_coldens": None,
                 },
             ),
+            "sbeta_onavg": PlotRule(
+                self,
+                partial(
+                    self._plot,
+                    "/comp/sbeta_onavg/beta",
+                    "/comp/sbeta_onavg/slope",
+                    yerr="/comp/sbeta_onavg/stderr",
+                ),
+                "Slope of the time averaged Sigma-PDF against cooling beta factor",
+                kind="comp",
+                dependencies=["sbeta_onavg"],
+            ),
             "sink_mass": PlotRule(
                 self,
                 partial(

postprocessor.py

@@ -911,12 +911,11 @@ class PostProcessor(HDF5Container):
             & (fluct_map > 0)
         )
 
-        nb_cells = np.sum(mask_pdf.flatten())
         values, edges = np.histogram(
             np.log10(fluct_map[mask_pdf].flatten()),
             self.pp_params.pdf.nb_bin,
             range=self.pp_params.pdf.range,
-            weights=np.ones(nb_cells) / nb_cells,
+            density=True,
         )
         centers = 0.5 * (edges[1:] + edges[:-1])
         return np.stack([values, centers])
@@ -1518,7 +1517,8 @@ class PostProcessor(HDF5Container):
             "P_avg",
             "T_avg",
             "P_mwavg",
-            "T_mwavg" "alpha_disk",
+            "T_mwavg",
+            "alpha_disk",
             "alpha_grav",
         ]
         for name in averageables:

pp_params.yml

@@ -22,10 +22,11 @@ disk: # Disk speficic parameters
 
 
 pdf: # parameters for probability density functions
-    nb_bin    : 100                   # Number of bins for the PDF
+    nb_bin    : 200                   # Number of bins for the PDF
     range     : [-1.5, 2.5]           # Range of the PDF (log of fluctuation)
-    xmin_fit  : 0.                    # Lower boundary of the fit (log of fluctuation)
-    xmax_fit  : 1.25                  # Upper boundary of the fit (log of fluctuation)
+    xmin_fit  : 0.1                   # Lower boundary of the fit (log of fluctuation)
+    xmax_fit  : 1.5                   # Upper boundary of the fit (log of fluctuation)
+    fit_cut   : 1e-4                  # Exclude value that are < fit_cut * maximum
 
 
 filaments:  # parameters for FilFinder