"""
Provides the :class:`~martini.sources.eagle_source.EAGLESource` class.
Facilitates working with EAGLE simulations as input.
"""
import numpy as np
from scipy.spatial.transform import Rotation
from typing import TYPE_CHECKING
from .sph_source import SPHSource
from ..sph_kernels import _WendlandC2Kernel, find_fwhm
from ..L_coords import L_coords
from os.path import join, normpath, sep
import astropy.units as U
from astropy.coordinates import ICRS
if TYPE_CHECKING:
from astropy.coordinates.builtin_frames.baseradec import BaseRADecFrame
[docs]
class EAGLESource(SPHSource):
"""
Class abstracting HI sources from publicly available EAGLE snapshot and group data.
For file access, see http://icc.dur.ac.uk/Eagle/database.php.
Parameters
----------
snapPath : str
Directory containing snapshot files. The directory structure unpacked
from the publicly available tarballs is expected; removing/renaming
files or directories below this will cause errors.
snapBase : str
Filename of snapshot files, omitting portion ``'.X.hdf5'``.
fof : int
FOF group number of the target object. Note that all particles in the
FOF group to which the subhalo belongs are used to construct the data
cube. This avoids strange "holes" at the locations of other subhaloes
in the same group, and gives a more realistic treatment of foreground
and background emission local to the source. In the EAGLE database,
this is the 'GroupNumber'.
sub : int
Subfind subhalo number of the target object. For centrals the subhalo
number is 0, for satellites >0. In the EAGLE database, this is then
'SubGroupNumber'.
db_user : str
Database username.
db_key : str, optional
Database password, or omit for a prompt at runtime.
subBoxSize : ~astropy.units.Quantity
:class:`~astropy.units.Quantity`, with dimensions of length.
Box half-side length of a region to load around the object of interest,
in physical (not comoving, no little h) units. Using larger values
will include more foreground/background, which may be desirable, but
will also slow down execution and impair the automatic routine used
to find a disc plane.
distance : ~astropy.units.Quantity
:class:`~astropy.units.Quantity`, with dimensions of length.
Source distance, also used to set the velocity offset via Hubble's law.
vpeculiar : ~astropy.units.Quantity, optional
:class:`~astropy.units.Quantity`, with dimensions of velocity.
Source peculiar velocity along the direction to the source centre.
rotation : ~scipy.spatial.transform.Rotation, optional
A rotation to apply to the source particles, specified using the
:class:`~scipy.spatial.transform.Rotation` class. That class supports many ways to
specify a rotation (Euler angle, rotation matrices, quaternions, etc.). Refer to
the :mod:`scipy` documentation for details. Note that the ``y-z`` plane will be
the one eventually placed in the plane of the "sky". Cannot be used at the same
time as ``L_coords``.
L_coords : ~martini.L_coords.L_coords, optional
A named tuple specifying 3 angles. Import it as ``from martini import L_coords``.
The angles are used to orient the galaxy relative to its angular momentum vector,
"L". The routine will first identify a preferred plane based on the angular
momenta of the central 1/3 of HI gas. This plane will then be rotated to lie in
the plane of the "sky" (``y-z`` plane), rotated by an angle ``az_rot`` around the
angular momentum vector (rotation around ``x``), then inclined by ``incl`` towards
or away from the line of sight (rotation around ``y``) and finally rotated on the
sky to set the position angle ``pa`` (second rotation around ``x``). All rotations
are extrinsic. The position angle refers to the receding side of the galaxy
measured East of North. The angles should be specified using syntax like:
``L_coords=L_coords(incl=0 * U.deg, pa=270 * U.deg, az_rot=0 * U.deg)``. These
example values are the defaults. Cannot be used at the same time as ``rotation``.
ra : ~astropy.units.Quantity, optional
:class:`~astropy.units.Quantity`, with dimensions of angle.
Right ascension for the source centroid.
dec : ~astropy.units.Quantity, optional
:class:`~astropy.units.Quantity`, with dimensions of angle.
Declination for the source centroid.
coordinate_frame : ~astropy.coordinates.builtin_frames.baseradec.BaseRADecFrame, \
optional
The coordinate frame assumed in converting particle coordinates to RA and Dec, and
for transforming coordinates and velocities to the data cube frame. The frame
needs to have a well-defined velocity as well as spatial origin. Recommended
frames are :class:`~astropy.coordinates.GCRS`, :class:`~astropy.coordinates.ICRS`,
:class:`~astropy.coordinates.HCRS`, :class:`~astropy.coordinates.LSRK`,
:class:`~astropy.coordinates.LSRD` or :class:`~astropy.coordinates.LSR`. The frame
should be passed initialized, e.g. ``ICRS()`` (not just ``ICRS``).
print_query : bool, optional
If True, the SQL query submitted to the EAGLE database is printed.
"""
def __init__(
self,
*,
snapPath: str,
snapBase: str,
fof: int,
sub: int,
db_user: str,
db_key: str | None = None,
subBoxSize: U.Quantity[U.kpc] = 50.0 * U.kpc,
distance: U.Quantity[U.Mpc],
vpeculiar: U.Quantity[U.km / U.s] = 0 * U.km / U.s,
rotation: Rotation | None = None,
L_coords: L_coords | None = None,
ra: U.Quantity[U.deg] = 0.0 * U.deg,
dec: U.Quantity[U.deg] = 0.0 * U.deg,
coordinate_frame: "BaseRADecFrame" = ICRS(),
print_query: bool = False,
) -> None:
# optional dependencies for this source class
from eagleSqlTools import connect, execute_query
from pyread_eagle import EagleSnapshot
from Hdecompose.atomic_frac import atomic_frac
import h5py
snapNum = int(snapBase.split("_")[1])
volCode = normpath(snapPath).split(sep)[-2]
query = (
"SELECT "
" sh.redshift as redshift, "
" sh.CentreOfPotential_x as x, "
" sh.CentreOfPotential_y as y, "
" sh.CentreOfPotential_z as z, "
" sh.Velocity_x as vx, "
" sh.Velocity_y as vy, "
" sh.Velocity_z as vz "
"FROM "
" {:s}_SubHalo as sh ".format(volCode)
+ "WHERE sh.Snapnum = {:d} ".format(snapNum)
+ " and sh.GroupNumber = {:d} ".format(fof)
+ " and sh.SubGroupNumber = {:d}".format(sub)
)
if print_query:
print("-----EAGLE-DB-QUERY-----")
print(query)
print("-------QUERY-ENDS-------")
if db_key is None:
from getpass import getpass
db_key = getpass("EAGLE database password:")
q = execute_query(connect(db_user, db_key), query)
redshift = q["redshift"]
a = np.power(1 + redshift, -1)
cop = np.array([q[coord] for coord in "xyz"]) * a * U.Mpc
vcent = np.array([q["v" + coord] for coord in "xyz"]) * U.km / U.s
snapFile = join(snapPath, snapBase + ".0.hdf5")
with h5py.File(snapFile, "r") as f:
h = f["RuntimePars"].attrs["HubbleParam"]
subBoxSize = (subBoxSize * h / a).to(U.Mpc).value
centre = (cop * h / a).to(U.Mpc).value
eagle_data = EagleSnapshot(snapFile)
region = np.vstack((centre - subBoxSize, centre + subBoxSize)).T.flatten()
eagle_data.select_region(*region)
lbox = f["/Header"].attrs["BoxSize"] * U.Mpc / h
fH = f["/RuntimePars"].attrs["InitAbundance_Hydrogen"]
fHe = f["/RuntimePars"].attrs["InitAbundance_Helium"]
proton_mass = f["/Constants"].attrs["PROTONMASS"] * U.g
mu = 1 / (fH + 0.25 * fHe)
gamma = f["/RuntimePars"].attrs["EOS_Jeans_GammaEffective"]
T0 = f["/RuntimePars"].attrs["EOS_Jeans_TempNorm_K"] * U.K
def fetch(att: str, ptype: int = 0) -> np.ndarray:
# gas is type 0, only need gas properties
tmp = eagle_data.read_dataset(ptype, att)
dset = f["/PartType{:d}/{:s}".format(ptype, att)]
aexp = dset.attrs.get("aexp-scale-exponent")
hexp = dset.attrs.get("h-scale-exponent")
return tmp[()] * np.power(a, aexp) * np.power(h, hexp)
code_to_g = f["/Units"].attrs["UnitMass_in_g"] * U.g
code_to_cm = f["/Units"].attrs["UnitLength_in_cm"] * U.cm
code_to_cm_s = f["/Units"].attrs["UnitVelocity_in_cm_per_s"] * U.cm / U.s
ng_g = fetch("GroupNumber")
particles = {
"xyz_g": (fetch("Coordinates") * code_to_cm).to(U.kpc),
"vxyz_g": (fetch("Velocity") * code_to_cm_s).to(U.km / U.s),
"T_g": fetch("Temperature") * U.K,
"hsm_g": (fetch("SmoothingLength") * code_to_cm).to(U.kpc)
* find_fwhm(_WendlandC2Kernel().kernel),
}
rho_g = fetch("Density") * U.g * U.cm**-3
SFR_g = fetch("StarFormationRate")
Habundance_g = fetch("ElementAbundance/Hydrogen")
m_g = fetch("Mass") # raw, converted in mHI_g below
particles["mHI_g"] = (
m_g
* atomic_frac(
redshift,
rho_g * Habundance_g / (mu * proton_mass),
particles["T_g"],
rho_g,
Habundance_g,
onlyA1=True,
EAGLE_corrections=True,
SFR=SFR_g,
mu=mu,
gamma=gamma,
fH=fH,
T0=T0,
)
* code_to_g
).to(U.Msun)
mask = ng_g == fof
for k, v in particles.items():
particles[k] = v[mask]
particles["xyz_g"] -= cop
particles["xyz_g"][particles["xyz_g"] > lbox / 2.0] -= lbox.to(U.kpc)
particles["xyz_g"][particles["xyz_g"] < -lbox / 2.0] += lbox.to(U.kpc)
particles["vxyz_g"] -= vcent
super().__init__(
distance=distance,
vpeculiar=vpeculiar,
rotation=rotation,
L_coords=L_coords,
ra=ra,
dec=dec,
h=h,
coordinate_frame=coordinate_frame,
**particles,
)
return