[postprocessor] add a mode to compute Q with generic units

This commit is contained in:
Noe Brucy
2021-06-23 13:56:06 +02:00
parent f0bea238c1
commit 5e08d1e751
+16 -9
View File
@@ -944,26 +944,34 @@ class PostProcessor(HDF5Container):
dmap_omega = rt_omega.process(self._cam[ax_los]).map.T dmap_omega = rt_omega.process(self._cam[ax_los]).map.T
return dmap_omega return dmap_omega
def _toomreQ_disk(self, ax_los, omega_approx=True): def _toomreQ_disk(self, ax_los, omega_approx=True, G1_units=True):
""" """
Compute the Toomre Q parameter Compute the Toomre Q parameter
""" """
# Get maps # Get maps
cs2 = self.get_value("/maps/T_mwavg_z") if G1_units:
coldens = self.get_value("/maps/coldens_z") G = 1.
omega = self.get_value("/maps/omega_mwavg_z") cs2 = self.get_value("/maps/T_mwavg_z")
coldens = self.get_value("/maps/coldens_z")
omega = self.get_value("/maps/omega_mwavg_z")
space_coeff = self.lbox
else:
G = U.G.express(U.kg**-1 * U.m ** 3 * U.s ** -2)
cs2 = self.get_value("/maps/T_mwavg_z", unit=U.m ** 2 * U.s ** -2)
coldens = self.get_value("/maps/coldens_z", unit=U.kg * U.m ** -2)
omega = self.get_value("/maps/omega_mwavg_z", unit=U.s ** -1)
space_coeff = self.info["unit_length"].express(U.m)
# Compute Q # Compute Q
if omega_approx: if omega_approx:
# Use angular frequency as epiciclic frequency (true if the disk is Keplerian) # Use angular frequency as epiciclic frequency (true if the disk is Keplerian)
kappa = omega kappa = omega
else: else:
# Get coordinates # Get coordinates
im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * self.lbox im_extent = np.array(self.save.root.maps._v_attrs.im_extent) * space_coeff
map_size = self.pp_params.pymses.map_size map_size = self.pp_params.pymses.map_size
center = np.array(self.pp_params.disk.center) * self.lbox center = np.array(self.pp_params.disk.center) * space_coeff
# Physical size of cells # Physical size of cells
dx = (im_extent[1] - im_extent[0]) / map_size dx = (im_extent[1] - im_extent[0]) / map_size
@@ -982,8 +990,7 @@ class PostProcessor(HDF5Container):
kappa_square = (2 * omega / rr) * Gr kappa_square = (2 * omega / rr) * Gr
kappa = np.sqrt(kappa_square) kappa = np.sqrt(kappa_square)
map_Q = (np.sqrt(cs2) * kappa) / (np.pi * self.G * coldens) map_Q = (np.sqrt(cs2) * kappa) / (np.pi * G * coldens)
return map_Q return map_Q
def _radial_bins(self, ax_los="z"): def _radial_bins(self, ax_los="z"):