import numpy as np
from .derived_field import ValidateParameter
from .field_plugin_registry import register_field_plugin
from .vector_operations import create_magnitude_field
[docs]
@register_field_plugin
def setup_astro_fields(registry, ftype="gas", slice_info=None):
unit_system = registry.ds.unit_system
pc = registry.ds.units.physical_constants
# slice_info would be the left, the right, and the factor.
# For example, with the old Enzo-ZEUS fields, this would be:
# slice(None, -2, None)
# slice(1, -1, None)
# 1.0
# Otherwise, we default to a centered difference.
if slice_info is None:
sl_left = slice(None, -2, None)
sl_right = slice(2, None, None)
div_fac = 2.0
else:
sl_left, sl_right, div_fac = slice_info
def _dynamical_time(field, data):
"""
sqrt(3 pi / (16 G rho))
"""
return np.sqrt(3.0 * np.pi / (16.0 * pc.G * data[ftype, "density"]))
registry.add_field(
(ftype, "dynamical_time"),
sampling_type="local",
function=_dynamical_time,
units=unit_system["time"],
)
def _jeans_mass(field, data):
MJ_constant = (((5.0 * pc.kboltz) / (pc.G * pc.mh)) ** 1.5) * (
3.0 / (4.0 * np.pi)
) ** 0.5
u = (
MJ_constant
* (
(data[ftype, "temperature"] / data[ftype, "mean_molecular_weight"])
** 1.5
)
* (data[ftype, "density"] ** (-0.5))
)
return u
registry.add_field(
(ftype, "jeans_mass"),
sampling_type="local",
function=_jeans_mass,
units=unit_system["mass"],
)
def _emission_measure(field, data):
dV = data[ftype, "mass"] / data[ftype, "density"]
nenhdV = data[ftype, "H_nuclei_density"] * dV
nenhdV *= data[ftype, "El_number_density"]
return nenhdV
registry.add_field(
(ftype, "emission_measure"),
sampling_type="local",
function=_emission_measure,
units=unit_system["number_density"],
)
def _mazzotta_weighting(field, data):
# Spectroscopic-like weighting field for galaxy clusters
# Only useful as a weight_field for temperature, metallicity, velocity
ret = data[ftype, "El_number_density"].d ** 2
ret *= data[ftype, "kT"].d ** -0.75
return ret
registry.add_field(
(ftype, "mazzotta_weighting"),
sampling_type="local",
function=_mazzotta_weighting,
units="",
)
def _optical_depth(field, data):
return data[ftype, "El_number_density"] * pc.sigma_thompson
registry.add_field(
(ftype, "optical_depth"),
sampling_type="local",
function=_optical_depth,
units=unit_system["length"] ** -1,
)
def _sz_kinetic(field, data):
# minus sign is because radial velocity is WRT viewer
# See issue #1225
return -data[ftype, "velocity_los"] * data[ftype, "optical_depth"] / pc.clight
registry.add_field(
(ftype, "sz_kinetic"),
sampling_type="local",
function=_sz_kinetic,
units=unit_system["length"] ** -1,
validators=[ValidateParameter("axis", {"axis": [0, 1, 2]})],
)
def _szy(field, data):
kT = data[ftype, "kT"] / (pc.me * pc.clight * pc.clight)
return data[ftype, "optical_depth"] * kT
registry.add_field(
(ftype, "szy"),
sampling_type="local",
function=_szy,
units=unit_system["length"] ** -1,
)
def _entropy(field, data):
mgammam1 = -2.0 / 3.0
tr = data[ftype, "kT"] * data[ftype, "El_number_density"] ** mgammam1
return data.apply_units(tr, field.units)
registry.add_field(
(ftype, "entropy"), sampling_type="local", units="keV*cm**2", function=_entropy
)
def _lorentz_factor(field, data):
b2 = data[ftype, "velocity_magnitude"].to_value("c")
b2 *= b2
return 1.0 / np.sqrt(1.0 - b2)
registry.add_field(
(ftype, "lorentz_factor"),
sampling_type="local",
units="",
function=_lorentz_factor,
)
# 4-velocity spatial components
def four_velocity_xyz(u):
def _four_velocity(field, data):
return data["gas", f"velocity_{u}"] * data["gas", "lorentz_factor"]
return _four_velocity
for u in registry.ds.coordinates.axis_order:
registry.add_field(
("gas", f"four_velocity_{u}"),
sampling_type="local",
function=four_velocity_xyz(u),
units=unit_system["velocity"],
)
# 4-velocity t-component
def _four_velocity_t(field, data):
return data["gas", "lorentz_factor"] * pc.clight
registry.add_field(
("gas", "four_velocity_t"),
sampling_type="local",
function=_four_velocity_t,
units=unit_system["velocity"],
)
create_magnitude_field(
registry,
"four_velocity",
unit_system["velocity"],
ftype=ftype,
)