diff --git a/plotter.py b/plotter.py index e741397..e074459 100644 --- a/plotter.py +++ b/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), - dependencies=["radial_centers", "rad_avg_" + name], + "Azimuthal mass weighted average of {}".format(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"), diff --git a/postprocessor.py b/postprocessor.py index 4017e9e..8009518 100644 --- a/postprocessor.py +++ b/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"): - radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() - bins_on_map = self.save.get_node("/maps/bins_on_map_" + ax_los).read() - dmap = self.save.get_node("/maps/" + name + "_" + ax_los).read() + def _rad_avg(self, name, ax_los="z", mass_weighted=False): + radial_bins = self.get_value("/radial/radial_bins_" + ax_los) + bins_on_map = self.get_value("/maps/bins_on_map_" + ax_los) + 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() - bins_on_map = self.save.get_node("/maps/bins_on_map_" + ax_los).read() - rr = self.save.get_node("/maps/rr_" + ax_los).read() - mean_bin = self.save.get_node("/radial/rad_avg_" + name + "_" + ax_los).read() + radial_bins = self.get_value("/radial/radial_bins_" + ax_los) + bins_on_map = self.get_value("/maps/bins_on_map_" + ax_los) + rr = self.get_value("/maps/rr_" + ax_los) + + 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"): - fluct_map = self.save.get_node("/maps/fluct_" + name + "_" + ax_los).read() - rr = self.save.get_node("/maps/rr_" + ax_los).read() + def _rad_fluct_pdf(self, name, ax_los="z", mass_weighted=False): + if mass_weighted: + 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),