diff --git a/comparator.py b/comparator.py index e74faf8..0410110 100644 --- a/comparator.py +++ b/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, diff --git a/plotter.py b/plotter.py index 2b2a5e5..9087c12 100644 --- a/plotter.py +++ b/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( diff --git a/postprocessor.py b/postprocessor.py index 2d82c93..e6a986d 100644 --- a/postprocessor.py +++ b/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: diff --git a/pp_params.yml b/pp_params.yml index bf95340..b81239c 100644 --- a/pp_params.yml +++ b/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