"""Analysis involving angular momentum.
"""
import logging
import numpy as np
from .. import array, filt, transformation, units, util
from . import halo
logger = logging.getLogger('pynbody.analysis.angmom')
[docs]
def ang_mom_vec(snap):
"""
Calculates the angular momentum vector of the specified snapshot.
Parameters
----------
snap : SimSnap
The snapshot to analyze
Returns
-------
array.SimArray
The angular momentum vector of the snapshot
"""
angmom = (snap['mass'].reshape((len(snap), 1)) *
np.cross(snap['pos'], snap['vel'])).sum(axis=0).view(np.ndarray)
result = angmom.view(array.SimArray)
result.units = snap['mass'].units * snap['pos'].units * snap['vel'].units
return result
[docs]
@util.deprecated("ang_mom_vec_units is deprecated. Use ang_mom_vec instead.")
def ang_mom_vec_units(snap):
"""Deprecated alias for ang_mom_vec"""
return ang_mom_vec(snap)
[docs]
def spin_parameter(snap):
"""Return the spin parameter lambda' of a centered halo
The spin parameter is defined as in eq. (5) of Bullock et al. 2001 (2001MNRAS.321..559B).
Note that the halo has to be aligned such that its center is the coordinate
origin and the velocity must be zeroed. If you are not sure whether this will
be true, calculate the spin parameter of h using:
>>> with pynbody.analysis.angmom.faceon(h):
>>> spin = pynbody.analysis.angmom.spin_parameter(h)
Parameters
----------
snap : SimSnap
The snapshot to analyze
Returns
-------
float
The dimensionless spin parameter lambda' of the halo
"""
m3 = snap['mass'].sum()
m3 = m3 * m3 * m3
l = np.sqrt(((ang_mom_vec_units(snap) ** 2).sum()) /
(2 * units.G * m3 * snap['r'].max()))
return float(l.in_units('1', **snap.conversion_context()))
[docs]
def calc_sideon_matrix(angmom_vec, along = [1.0, 0.0, 0.0]):
"""Calculate the rotation matrix to put the specified angular momentum vector side-on.
The rotation matrix is calculated such that the angular momentum vector will be placed in the y direction
post-transformation.
Parameters
----------
angmom_vec : array_like
The angular momentum vector that will be placed in the y-direction post-transformation.
along : array_like
An additional orientation vector. The components of this vector perpendicular to angmom_vec defines the
direction to transform to 'along', i.e. to the positive x-axis post-transformation.
Returns
-------
array_like
The rotation matrix
"""
vec_in = np.asarray(angmom_vec)
vec_in = vec_in / np.sum(vec_in ** 2).sum() ** 0.5
vec_p1 = np.cross(along, vec_in)
vec_p1 = vec_p1 / np.sum(vec_p1 ** 2).sum() ** 0.5
vec_p2 = np.cross(vec_in, vec_p1)
matr = np.concatenate((vec_p2, vec_in, vec_p1)).reshape((3, 3))
return matr
[docs]
def calc_faceon_matrix(angmom_vec, up=[0.0, 1.0, 0.0]):
"""Calculate the rotation matrix to put the specified angular momentum vector face-on.
The rotation matrix is calculated such that the angular momentum vector will be placed in the z direction
post-transformation.
Parameters
----------
angmom_vec : array_like
The angular momentum vector that will be placed in the z-direction post-transformation.
up : array_like
An additional orientation vector. The components of this vector perpendicular to angmom_vec defines the
direction to transform to 'up', i.e. to the positive y-axis post-transformation.
Returns
-------
array_like
The rotation matrix
"""
vec_in = np.asarray(angmom_vec)
vec_in = vec_in / np.sum(vec_in ** 2).sum() ** 0.5
vec_p1 = np.cross(up, vec_in)
vec_p1 = vec_p1 / np.sum(vec_p1 ** 2).sum() ** 0.5
vec_p2 = np.cross(vec_in, vec_p1)
matr = np.concatenate((vec_p1, vec_p2, vec_in)).reshape((3, 3))
return matr
[docs]
def align(h, vec_to_xform, disk_size="5 kpc", move_all=True, already_centered = False,
center_kwargs = None):
"""Reposition and rotate the ancestor of h to place the angular momentum into a specified orientation.
The routine first calls the center routine to reposition the halo (unless ``already_centered`` is ``True``).
If there are a sufficient number of star particles (more than 100), only the star particles are used for
centering; if there is an offset between dark matter and baryons, it is better to centre on the baryons.
Then, it determines the disk orientation using the angular momentum vector of the gas particles within
a specified radius of the halo center. If there are less than 5 gas particles within this radius, the
routine falls back first to stellar particles, and then to all particles.
Finally, the angular momentum vector is converted into a rotation matrix using the ``vec_to_xform``
function, and the rotation is applied.
.. seealso::
Convenient wrappers around this routine are provided to either align a disk side-on or face-on;
see :func:`sideon` and :func:`faceon`.
.. versionchanged :: 2.3.1
The routine previously checked the total number of gas particles to determine whether to use gas for
angular momentum determination, but now actually checks particles in the specified radius. This
prevents an issue where there could be enough gas particles overall but none in the specified radius,
leading the routine to fail.
Additionally, the docstring regarding choice of particles for centring has been clarified.
Parameters
----------
h : SimSnap
The portion of the simulation from which to extract a centre and orientation. Typically this is a galaxy halo.
vec_to_xform : function
The function to use to calculate the rotation matrix from the measured angular momentum vector.
disk_size : str | float, optional
The size of the disk to use for calculating the angular momentum vector. Default is "5 kpc".
move_all : bool, optional
If True, the ancestor simulation of *h* is transformed. If False, only *h* is moved. Default is True.
already_centered : bool, optional
If True, the simulation is assumed to be already centered. Default is False.
center_kwargs : dict, optional
Dictionary of additional keyword arguments to pass to the centering routine.
"""
if center_kwargs is None:
center_kwargs = {}
center_kwargs.update({'move_all': move_all})
if already_centered:
if move_all:
top = h.ancestor
else:
top = h
tx = transformation.NullTransformation(top)
else:
if len(h.st) > 100:
h_for_centering = h.st
else:
h_for_centering = h
tx = halo.center(h_for_centering, **center_kwargs)
try:
if len(cen_trial := h.gas[filt.Sphere(disk_size)]) > 5:
cen = cen_trial
elif len(cen_trial := h.st[filt.Sphere(disk_size)]) > 5:
cen = cen_trial
else:
cen = h[filt.Sphere(disk_size)]
logger.info("Calculating angular momentum vector...")
trans = vec_to_xform(ang_mom_vec(cen))
logger.info("Transforming simulation...")
tx = tx.rotate(trans)
logger.info("...done!")
return tx
except:
tx.revert()
raise
[docs]
def sideon(h, **kwargs):
"""Reposition and rotate the ancestor of h to place the disk edge-on (i.e. into the x-z plane).
Since pynbody's imaging routines project along the z direction, one can get a side-on view of a disk or
other rotationally-supported structure by calling this routine first.
.. seealso::
For details of how the transformation is calculated, see the documentation for the underlying :func:`align` routine.
A similar routine for aligning face-on is also provided; see :func:`faceon`.
Parameters
----------
h : SimSnap
The portion of the simulation from which to extract a centre and orientation. Typically this is a galaxy halo.
disk_size : str | float, optional
The size of the disk to use for calculating the angular momentum vector. Default is "5 kpc".
move_all : bool, optional
If True, the ancestor simulation of *h* is transformed. If False, only *h* is moved. Default is True.
already_centered : bool, optional
If True, the simulation is assumed to be already centered. Default is False.
center_kwargs : dict, optional
Dictionary of additional keyword arguments to pass to the centering routine.
"""
kwargs.update({'vec_to_xform': calc_sideon_matrix})
return align(h, **kwargs).set_description('sideon')
[docs]
def faceon(h, **kwargs):
"""Reposition and rotate the ancestor of h to place the disk face-on (i.e. into the x-y plane).
Since pynbody's imaging routines project along the z direction, one can get a face-on view of a disk
or other rotationally-supported structure by calling this routine first.
.. seealso::
For details of how the transformation is calculated, see the documentation for the underlying :func:`align` routine.
A similar routine for aligning side-on is also provided; see :func:`sideon`.
Parameters
----------
h : SimSnap
The portion of the simulation from which to extract a centre and orientation. Typically this is a galaxy halo.
disk_size : str | float, optional
The size of the disk to use for calculating the angular momentum vector. Default is "5 kpc".
move_all : bool, optional
If True, the ancestor simulation of *h* is transformed. If False, only *h* is moved. Default is True.
already_centered : bool, optional
If True, the simulation is assumed to be already centered. Default is False.
center_kwargs : dict, optional
Dictionary of additional keyword arguments to pass to the centering routine.
"""
kwargs.update({'vec_to_xform': calc_faceon_matrix})
return align(h, **kwargs).set_description('faceon')
