Source code for pynbody.analysis.angmom
"""
angmom
======
"""
import logging
import numpy as np
from .. import array, filt, units, config, transformation
from . import halo
logger = logging.getLogger('pynbody.analysis.angmom')
[docs]def ang_mom_vec(snap):
"""
Return the angular momentum vector of the specified snapshot.
The return units are [mass]*[dist]*[vel] as per the units of the
snapshot.
"""
angmom = (snap['mass'].reshape((len(snap), 1)) *
np.cross(snap['pos'], snap['vel'])).sum(axis=0).view(np.ndarray)
return angmom
[docs]def ang_mom_vec_units(snap):
"""
Return the angular momentum vector of the specified snapshot
with correct units.
Note that the halo has to be aligned such that the disk
is in the x-y-plane and its center must be the coordinate
origin.
"""
angmom = ang_mom_vec(snap)
return array.SimArray(angmom, snap['mass'].units * snap['pos'].units * snap['vel'].units)
[docs]def spin_parameter(snap):
"""
Return the spin parameter \lambda' of a centered halo
as defined in eq. (5) of Bullock et al. 2001
(2001MNRAS.321..559B).
Note that the halo has to be aligned such that the disk
is in the x-y-plane and its center must be the coordinate
origin.
"""
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()))
def calc_sideon_matrix(angmom_vec):
vec_in = np.asarray(angmom_vec)
vec_in = vec_in / np.sum(vec_in ** 2).sum() ** 0.5
vec_p1 = np.cross([1, 0, 0], 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
def calc_faceon_matrix(angmom_vec, up=[0.0, 1.0, 0.0]):
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 sideon(h, vec_to_xform=calc_sideon_matrix, cen_size="1 kpc",
disk_size="5 kpc", cen=None, vcen=None, move_all=True,
**kwargs):
"""
Reposition and rotate the simulation containing the halo h to see
h's disk edge on.
Given a simulation and a subview of that simulation (probably the
halo of interest), this routine centers the simulation and rotates
it so that the disk lies in the x-z plane. This gives a side-on
view for SPH images, for instance.
"""
global config
if move_all:
top = h.ancestor
else:
top = h
# Top is the top-level view of the simulation, which will be
# transformed
if cen is None:
logger.info("Finding halo center...")
# or h['pos'][h['phi'].argmin()]
cen = halo.center(h, retcen=True, **kwargs)
logger.info("... cen=%s" % cen)
tx = transformation.inverse_translate(top, cen)
if vcen is None:
vcen = halo.vel_center(h, retcen=True, cen_size=cen_size)
tx = transformation.inverse_v_translate(tx, vcen)
# Use gas from inner 10kpc to calculate angular momentum vector
if (len(h.gas) > 0):
cen = h.gas[filt.Sphere(disk_size)]
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 = transformation.transform(tx, trans)
logger.info("...done!")
return tx
[docs]def faceon(h, **kwargs):
"""
Reposition and rotate the simulation containing the halo h to see
h's disk face on.
Given a simulation and a subview of that simulation (probably the
halo of interest), this routine centers the simulation and rotates
it so that the disk lies in the x-y plane. This gives a face-on
view for SPH images, for instance.
"""
return sideon(h, calc_faceon_matrix, **kwargs)