"""
Calculations relevant to galactic morphology
.. versionadded :: 2.0
This module was added in version 2.0.
Previously, only the :func:`decomp` function was available, in its own module within :mod:`pynbody.analysis`.
"""
import logging
import numpy as np
from .. import filt, transformation, util
from . import angmom, profile
logger = logging.getLogger('pynbody.analysis.morphology')
[docs]
def estimate_jcirc_from_energy(h, particles_per_bin=1000, quantile=0.99):
"""Estimate the circular angular momentum as a function of energy.
This routine calculates the circular angular momentum as a function of energy
for the stars in the simulation, using a profile with a fixed number of particles
per bin. It then estimates a circular angular momentum for each individual particle
by interpolating the profile.
Arguments
---------
h : SimSnap
The simulation snapshot to analyze
quantile : float
The circular angular momentum will be estimated as the specified quantile of the scalar angular momentum
particles_per_bin : int
The approximate number of particles per bin in the profile. Default is 1000.
"""
nbins = len(h) // particles_per_bin
pro_d = profile.QuantileProfile(h, q=(quantile,), nbins=nbins, type='equaln', calc_x = lambda sim : sim['te'])
pro_d.create_particle_array("j2", particle_name='j_circ2', target_simulation=h)
h['j_circ'] = np.sqrt(h['j_circ2'])
del h.ancestor['j_circ2']
return pro_d
[docs]
def estimate_jcirc_from_rotation_curve(h, particles_per_bin=1000):
"""Estimate the circular angular momentum as a function of radius in the disk (x-y) plane.
This routine calculates the circular velocity as a function of radius for the disk, using a profile with a fixed
number of particles per radial bin. It then estimates a circular angular momentum for each individual particle by
interpolating the profile.
.. warning::
This routine is only valid for simulations where all the stars are anyway in quite a narrow disc. Otherwise
the interpolation back onto the individual particles carries limited meaning.
For more general cases, use :func:`estimate_jcirc_from_energy` instead.
Arguments
---------
h : SimSnap
The simulation snapshot to analyze
particles_per_bin : int
The approximate number of particles per bin in the profile. Default is 1000.
"""
d = h[filt.Disc('1 Mpc', h['eps'].min() * 3)]
nbins = len(d) // particles_per_bin
pro_d = profile.Profile(d, nbins=nbins, type='equaln')
pro_d.create_particle_array("j_circ", target_simulation=h)
[docs]
def decomp(h, aligned=False, j_disk_min=0.8, j_disk_max=1.1, E_cut=None, j_circ_from_r=False,
angmom_size="3 kpc", particles_per_bin = 500):
"""Creates an array 'decomp' for star particles in the simulation, with an integer specifying components.
The possible values of the components are:
1. thin disk
2. halo
3. bulge
4. thick disk
5. pseudo bulge
First, the simulation is aligned so that the disk is in the xy plane. Then, the maximum angular momentum of
stars as a function of energy is estimated, by default using the routine :func:`estimate_jcirc_from_energy`.
The stars are then classified into the components based on their energy and disk-aligned angular momentum component.
* The thin disk is defined as stars with angular momentum between j_disk_min and j_disk_max.
* From the remaining stars, a critical angular momentum j_crit is then calculated that separates rotating from
non-rotating components. By definition, this is chosen such that the mean rotaiton velocity of the non-rotating
part is zero
* The most tightly bound rotating component is labelled as the pseudo bulge, while less tightly bound rotating
stars are labelled as the thick disk, based on E_cut.
* The non-rotating stars are then separated into bulge and halo components based on their binding energy, again
based on E_cut.
This routine is based on an original IDL procedure by Chris Brook.
.. versionchanged :: 2.0
The routine now defaults to using a new method to estimate the circular angular momentum as a function of energy;
see :func:`estimate_jcirc_from_energy`.
The critical angular momentum for the spheroid is now determined by insisting the mean angular momentum of the
spheroid should be zero.
The above changes lead to different, but probably better, classifications compared with pynbody version 1.
Additionally, this routine is now inside the :mod:`pynbody.analysis.morphology` module.
Arguments
---------
h : SimSnap
The simulation snapshot to analyze
aligned : bool
If True, the simulation is assumed to be already aligned so that the disk is in the xy plane.
Otherwise, the simulation is recentered and aligned into the xy plane.
j_disk_min : float
The minimum angular momentum as a proportion of the circular angular momentum which a particle must have to be
part of the 'disk'.
j_disk_max : float
The maximum angular momentum as a proportion of the circular angular momentum which a particle can have to be
part of the 'disk'.
E_cut : float
The energy boundary between bulge and halo. If None, this is taken to be the median energy of the stars. Note
that the distinction between bulge and halo is somewhat arbitrary and may not be physically meaningful.
j_circ_from_r : bool
If True, the maximum angular momentum is determined as a function of disc radius, rather than as a function of
energy. Default False (determine as function of energy). This option is only valid for simulations where
all the stars are anyway in quite a narrow disc.
angmom_size : str
The size of the disk to use for calculating the angular momentum vector. Default is "3 kpc".
particles_per_bin : int
The approximate number of particles per bin in the profile. Default is 500.
"""
if aligned:
tx = transformation.NullTransformation(h)
else:
tx = angmom.faceon(h, disk_size=angmom_size)
with tx:
if j_circ_from_r:
estimate_jcirc_from_rotation_curve(h, particles_per_bin=particles_per_bin)
else:
estimate_jcirc_from_energy(h, particles_per_bin=particles_per_bin)
h['jz_by_jzcirc'] = h['j'][:, 2] / h['j_circ']
h_star = h.star
if 'decomp' not in h_star:
h_star._create_array('decomp', dtype=int)
disk = np.where(
(h_star['jz_by_jzcirc'] > j_disk_min) * (h_star['jz_by_jzcirc'] < j_disk_max))
h_star['decomp', disk[0]] = 1
# h_star = h_star[np.where(h_star['decomp']!=1)]
# Find disk/spheroid angular momentum cut-off to make spheroid
# angular momentum exactly zero
JzJcirc = h_star['jz_by_jzcirc']
te = h_star['te']
logger.info("Finding spheroid/disk angular momentum boundary...")
j_crit = util.bisect(-2.0, 2.0, lambda c: np.mean(JzJcirc[JzJcirc < c]))
logger.info("j_crit = %.2e" % j_crit)
if j_crit > j_disk_min:
logger.warning(
"!! j_crit exceeds j_disk_min. This is usually a sign that something is going wrong (train-wreck galaxy?)")
logger.warning("!! j_crit will be reset to j_disk_min=%.2e" % j_disk_min)
j_crit = j_disk_min
sphere = np.where(h_star['jz_by_jzcirc'] < j_crit)
if E_cut is None:
E_cut = np.median(h_star['te'])
logger.info("E_cut = %.2e" % E_cut)
halo = np.where((te > E_cut) * (JzJcirc < j_crit))
bulge = np.where((te <= E_cut) * (JzJcirc < j_crit))
pbulge = np.where((te <= E_cut) * (JzJcirc > j_crit)
* ((JzJcirc < j_disk_min) + (JzJcirc > j_disk_max)))
thick = np.where((te > E_cut) * (JzJcirc > j_crit)
* ((JzJcirc < j_disk_min) + (JzJcirc > j_disk_max)))
h_star['decomp', halo] = 2
h_star['decomp', bulge] = 3
h_star['decomp', thick] = 4
h_star['decomp', pbulge] = 5