Source code for yt.fields.cosmology_fields

from .derived_field import ValidateParameter
from .field_exceptions import NeedsConfiguration, NeedsParameter
from .field_plugin_registry import register_field_plugin


[docs] @register_field_plugin def setup_cosmology_fields(registry, ftype="gas", slice_info=None): unit_system = registry.ds.unit_system # 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 _matter_density(field, data): return data[ftype, "density"] + data[ftype, "dark_matter_density"] registry.add_field( (ftype, "matter_density"), sampling_type="local", function=_matter_density, units=unit_system["density"], ) def _matter_mass(field, data): return data[ftype, "matter_density"] * data["index", "cell_volume"] registry.add_field( (ftype, "matter_mass"), sampling_type="local", function=_matter_mass, units=unit_system["mass"], ) # rho_total / rho_cr(z). def _overdensity(field, data): if ( not hasattr(data.ds, "cosmological_simulation") or not data.ds.cosmological_simulation ): raise NeedsConfiguration("cosmological_simulation", 1) co = data.ds.cosmology return data[ftype, "matter_density"] / co.critical_density( data.ds.current_redshift ) registry.add_field( (ftype, "overdensity"), sampling_type="local", function=_overdensity, units="" ) # rho_baryon / <rho_baryon> def _baryon_overdensity(field, data): if ( not hasattr(data.ds, "cosmological_simulation") or not data.ds.cosmological_simulation ): raise NeedsConfiguration("cosmological_simulation", 1) omega_baryon = data.get_field_parameter("omega_baryon") if omega_baryon is None: raise NeedsParameter("omega_baryon") co = data.ds.cosmology # critical_density(z) ~ omega_lambda + omega_matter * (1 + z)^3 # mean matter density(z) ~ omega_matter * (1 + z)^3 return ( data[ftype, "density"] / omega_baryon / co.critical_density(0.0) / (1.0 + data.ds.current_redshift) ** 3 ) registry.add_field( (ftype, "baryon_overdensity"), sampling_type="local", function=_baryon_overdensity, units="", validators=[ValidateParameter("omega_baryon")], ) # rho_matter / <rho_matter> def _matter_overdensity(field, data): if ( not hasattr(data.ds, "cosmological_simulation") or not data.ds.cosmological_simulation ): raise NeedsConfiguration("cosmological_simulation", 1) co = data.ds.cosmology # critical_density(z) ~ omega_lambda + omega_matter * (1 + z)^3 # mean density(z) ~ omega_matter * (1 + z)^3 return ( data[ftype, "matter_density"] / data.ds.omega_matter / co.critical_density(0.0) / (1.0 + data.ds.current_redshift) ** 3 ) registry.add_field( (ftype, "matter_overdensity"), sampling_type="local", function=_matter_overdensity, units="", ) # r / r_vir def _virial_radius_fraction(field, data): virial_radius = data.get_field_parameter("virial_radius") if virial_radius == 0.0: ret = 0.0 else: ret = data[("index", "radius")] / virial_radius return ret registry.add_field( ("index", "virial_radius_fraction"), sampling_type="local", function=_virial_radius_fraction, validators=[ValidateParameter("virial_radius")], units="", ) # Weak lensing convergence. # Eqn 4 of Metzler, White, & Loken (2001, ApJ, 547, 560). # This needs to be checked for accuracy. def _weak_lensing_convergence(field, data): if ( not hasattr(data.ds, "cosmological_simulation") or not data.ds.cosmological_simulation ): raise NeedsConfiguration("cosmological_simulation", 1) co = data.ds.cosmology pc = data.ds.units.physical_constants observer_redshift = data.get_field_parameter("observer_redshift") source_redshift = data.get_field_parameter("source_redshift") # observer to lens dl = co.angular_diameter_distance(observer_redshift, data.ds.current_redshift) # observer to source ds = co.angular_diameter_distance(observer_redshift, source_redshift) # lens to source dls = co.angular_diameter_distance(data.ds.current_redshift, source_redshift) # removed the factor of 1 / a to account for the fact that we are projecting # with a proper distance. return ( 1.5 * (co.hubble_constant / pc.clight) ** 2 * (dl * dls / ds) * data[ftype, "matter_overdensity"] ).in_units("1/cm") registry.add_field( (ftype, "weak_lensing_convergence"), sampling_type="local", function=_weak_lensing_convergence, units=unit_system["length"] ** -1, validators=[ ValidateParameter("observer_redshift"), ValidateParameter("source_redshift"), ], )