Analysis

angmom

pynbody.analysis.angmom.ang_mom_vec(snap)[source]

Return the angular momentum vector of the specified snapshot.

The return units are [mass]*[dist]*[vel] as per the units of the snapshot.

pynbody.analysis.angmom.ang_mom_vec_units(snap)[source]

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.

pynbody.analysis.angmom.spin_parameter(snap)[source]

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.

pynbody.analysis.angmom.sideon(h, vec_to_xform=<function calc_sideon_matrix at 0x1082a37d0>, cen_size='1 kpc', disk_size='5 kpc', cen=None, vcen=None, move_all=True, **kwargs)[source]

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.

pynbody.analysis.angmom.faceon(h, **kwargs)[source]

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.

cosmology

A set of functions for common cosmological calculations.

pynbody.analysis.cosmology.hzoverh0(a, omegam0)[source]

returns: H(a) / H0 = [omegam/a**3 + (1-omegam)]**0.5

pynbody.analysis.cosmology.linear_growth_factor(f, z=None)[source]

Calculate the linear growth factor b(a), normalized to 1 at z=0, for the cosmology of snapshot f.

The output is dimensionless. If a redshift z is specified, it is used in place of the redshift in output f.

pynbody.analysis.cosmology.rate_linear_growth(f, z=None, unit='h Gyr^-1')[source]

Calculate the linear growth rate b’(a), normalized to 1 at z=0, for the cosmology of snapshot f.

The output is in ‘h Gyr^-1’ by default. If a redshift z is specified, it is used in place of the redshift in output f.

pynbody.analysis.cosmology.age(f, z=None, unit='Gyr')[source]

Calculate the age of the universe in the snapshot f by integrating the Friedmann equation.

The output is given in the specified units. If a redshift z is specified, it is used in place of the redshift in the output f.

Input:

f: SimSnap

Optional Keywords:

z (None): desired redshift. Can be a single number, a list, or a numpy.ndarray.

unit (‘Gyr’): desired units for age output

pynbody.analysis.cosmology.redshift(*fn_args, **fn_kwargs)[source]

Calculate the redshift given a snapshot and a time since Big Bang in Gyr.

Uses scipy.optimize.newton to do the root finding if number of elements in the time array is less than 1000, otherwise uses a linear interpolation.

Input:

f: SimSnap with cosmological parameters defined

time: time since the Big Bang in Gyr for which a redshift should
be returned. float, list, or numpy.ndarray
pynbody.analysis.cosmology.rho_crit(f, z=None, unit=None)[source]

Calculate the critical density of the universe in the snapshot f.

z specifies the redshift. If z is none, the redshift of the provided snapshot is used.

unit specifies the units of the returned density. If unit is None, the returned density will be in the units of f[“mass”].units/f[“pos”].units**3. If that unit cannot be calculated, the returned units are Msol kpc^-3 comoving.

Note that you can get slightly confusing results if your simulation is in comoving units and you specify a different redshift z. Specifically, the physical density for the redshift you specify is calulated, but expressed as a comoving density at the redshift of the snapshot. This is intentional behaviour.

pynbody.analysis.cosmology.rho_M(f, z=None, unit=None)[source]

Calculate the matter density of the universe in snapshot f.

unit and z are used if not None, as by rho_crit. See also the note in rho_crit about confusion over comoving units in this case.

pynbody.analysis.cosmology.H(f)[source]

Calculate the Hubble parameter of the universe in snapshot f

pynbody.analysis.cosmology.add_hubble(f)[source]

Add the hubble flow to velocities in snapshot f

pynbody.analysis.cosmology.comoving_to_physical(ar)[source]

Given an array, modify it to be in physical units (remove any dependencies on a or aform).

decomp

Tools for bulge/disk/halo decomposition

pynbody.analysis.decomp.decomp(h, aligned=False, j_disk_min=0.8, j_disk_max=1.1, E_cut=None, j_circ_from_r=False, cen=None, vcen=None, log_interp=False, angmom_size='3 kpc')[source]

Creates an array ‘decomp’ for star particles in the simulation, with an integer specifying a kinematic decomposition. The possible values are:

1 – thin disk

2 – halo

3 – bulge

4 – thick disk

5 – pseudo bulge

This routine is based on an original IDL procedure by Chris Brook.

Parameters:

h – the halo to work on

j_disk_min – 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 – 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 – the energy boundary between bulge and spheroid. If
None, this is taken to be the median energy of the stars.
aligned – if False, the simulation is recenterd and aligned so
the disk is in the xy plane as required for the rest of the analysis.
cen – if not None, specifies the center of the halo. Otherwise
it is found. This has no effect if aligned=True
vcen – if not None, specifies the velocity center of the
halo. Otherwise it is found. This has no effect if aligned=True

j_circ_from_r – if True, the maximum angular momentum is determined as a function of radius, rather than as a function of orbital energy. Default False (determine as function of energy).

angmom_size – the size of the gas sphere used to determine the
plane of the disk

gravity

Functions for calculating the gravitational potential and accelerations.

pynbody.analysis.gravity.potential(f, pos_vec, eps=None, unit=None)[source]

Calculates the gravitational potential at the specified position using all particles in the specified snapshot.

The gravitational softening length is determined by (in order of preference):

  1. The parameter eps (scalar, unit or array)
  2. The array f[‘eps’]
  3. f.properties[‘eps’] (scalar or unit)
pynbody.analysis.gravity.accel(f, pos_vec, eps=None)[source]

Calculates the gravitational acceleration vector at the specified position using all particles in the specified snapshot.

The gravitational softening length is determined by (in order of preference):

  1. The parameter eps (scalar, unit or array)
  2. The array f[‘eps’]
  3. f.properties[‘eps’] (scalar or unit)
pynbody.analysis.gravity.midplane_rot_curve(f, rxy_points, eps=None)[source]

Calculates the midplane rotation curve at a specified set of points. Used primarily by pynbody.analysis.profile.v_circ()

pynbody.analysis.gravity.midplane_potential(f, rxy_points, eps=None)[source]

Calculates the midplane potential at a specified set of points. Used primarily by pynbody.analysis.profile.pot()

halo

Functions for dealing with and manipulating halos in simulations.

pynbody.analysis.halo.center_of_mass(sim)[source]

Return the centre of mass of the SimSnap

pynbody.analysis.halo.center_of_mass_velocity(sim)[source]

Return the center of mass velocity of the SimSnap

pynbody.analysis.halo.shrink_sphere_center(sim, r=None, shrink_factor=0.7, min_particles=100, verbose=False, num_threads=4, **kwargs)[source]

Return the center according to the shrinking-sphere method of Power et al (2003)

Input:

sim : a simulation snapshot - this can be any subclass of SimSnap

Optional Keywords:

r (default=None): initial search radius. This can be a string
indicating the unit, i.e. “200 kpc”, or an instance of Unit().
shrink_factor (default=0.7): the amount to shrink the search
radius by on each iteration
min_particles (default=100): minimum number of particles within
the search radius. When this number is reached, the search is complete.
verbose (default=False): if True, prints out the diagnostics at
each iteration. Useful to determine whether the centering is zeroing in on the wrong part of the simulation.
pynbody.analysis.halo.virial_radius(sim, cen=None, overden=178, r_max=None)[source]

Calculate the virial radius of the halo centered on the given coordinates.

This is here defined by the sphere centered on cen which contains a mean density of overden * rho_M_0 * (1+z)^3.

pynbody.analysis.halo.hybrid_center(sim, r='3 kpc', **kwargs)[source]

Determine the center of the halo by finding the shrink-sphere -center inside the specified distance of the potential minimum

pynbody.analysis.halo.index_center(sim, **kwargs)[source]

Determine the center of mass based on specific particles.

Supply a list of indices using the ind keyword.

pynbody.analysis.halo.vel_center(sim, mode=None, cen_size='1 kpc', retcen=False, move_all=True, **kwargs)[source]

Use stars from a sphere to calculate center of velocity. The size of the sphere is given by the cen_size keyword and defaults to 1 kpc.

Keyword arguments:

mode: reserved for future use; currently ignored

move_all: if True (default), move the entire snapshot. Otherwise only move the particles in the halo passed in.

retcen: if True only return the velocity center without moving the
snapshot (default = False)
pynbody.analysis.halo.center(sim, mode=None, retcen=False, vel=True, cen_size='1 kpc', move_all=True, wrap=False, **kwargs)[source]

Determine the center of mass of the given particles using the specified mode, then recenter the particles (of the entire ancestor snapshot) accordingly

Accepted values for mode are

pot: potential minimum

com: center of mass

ssc: shrink sphere center

ind: center on specific particles; supply the list of particles using the ind keyword.

hyb: for sane halos, returns the same as ssc, but works faster by
starting iteration near potential minimum

or a function returning the COM.

Other keywords:

retcen: if True only return the center without centering the
snapshot (default = False)
ind: only used when mode=ind – specifies the indices of
particles to be used for centering

vel: if True, translate velocities so that the velocity of the central 1kpc (default) is zeroed. Other values can be passed with cen_size.

move_all: if True (default), move the entire snapshot. Otherwise only move the particles in the halo passed in.

wrap: if True, pre-centre and wrap the simulation so that halos on the edge of the box are handled correctly. Default False.

halo mass function (hmf)

Various halo mass function routines.

pynbody.analysis.hmf.correlation_func(context, log_r_min=-3, log_r_max=2, delta_log_r=0.2, pspec=<class 'pynbody.analysis.hmf.PowerSpectrumCAMB'>)[source]

Calculate the linear density field correlation function.

Args:

context (SimSnap): The snapshot from which to pull the
cosmological context (includes sigma8 normalization and growth function integrations, but does not currently affect transfer function)

Kwargs:

log_r_min: log10 of the minimum separation (Mpc h^-1) to
consider
log_r_max: log10 of the maximum separation (Mpc h^-1) to
consider

delta_log_r: The value spacing in dex

pspec: A power spectrum object; default is a WMAP7 cosmology
calculated by CAMB.

Returns:

r: Array of the r values (Mpc h^-1) for which the correlation
function was evaluated.
Xi: Array of the dimensionless correlation for each
separation.
pynbody.analysis.hmf.f_press_schechter(nu)[source]

The Press-Schechter kernel used by halo_mass_function

pynbody.analysis.hmf.f_sheth_tormen(nu, Anorm=0.3222, a=0.707, p=0.3)[source]

Sheth & Tormen (1999) fit (see also Sheth Mo & Tormen 2001)

pynbody.analysis.hmf.f_reed_no_z(nu, deltac=1.68647)[source]

modified S-T fit by the G1 gaussian term and c

pynbody.analysis.hmf.f_reed_z_evo(nu, neff, deltac=1.68647)[source]

modified S-T fit by the n_eff dependence and the G1 and G2 gaussian terms and c where P(k) proportional to k_halo**(n_eff) and k_halo = Mhalo / r_halo_precollapse. eqn 13 of Reed et al 2007 estimtes neff = 6 d ln(1/sigma(M))/ d ln M - 3

pynbody.analysis.hmf.cole_kaiser_bias(nu, delta_c)[source]

The Cole-Kaiser (1989) bias function. Also in Mo & White 1996.

pynbody.analysis.hmf.sheth_tormen_bias(nu, delta_c, a=0.707, b=0.5, c=0.6)[source]

The Sheth-Tormen (1999) bias function [eq 8]

pynbody.analysis.hmf.halo_mass_function(context, log_M_min=8.0, log_M_max=15.0, delta_log_M=0.1, kern='ST', pspec=<class 'pynbody.analysis.hmf.PowerSpectrumCAMB'>, delta_crit=1.686, no_h=False)[source]

Returns the halo mass function, dN/d log_{10} M in units of Mpc^-3 h^3.

Args:

context (SimSnap): The snapshot from which to pull the cosmological context (includes sigma8 normalization and growth function integrations, but does not currently affect transfer function)

Kwargs:

log_M_min: The minimum halo mass (Msol h^-1) to consider

log_M_max: The maximum halo mass (Msol h^-1) to consider

delta_log_M: The bin spacing of halo masses (see warning below)

kern: The kernel function which dictates what type of mass function to calculate; or a string (“PS” or “ST”) for one of the defaults

pspec: A power spectrum object (which also defines the window function); default is a WMAP7 cosmology calculated by CAMB, and a top hat window

delta_crit: The critical overdensity for collapse

Returns:

M: The centre of the mass bins, in Msol h^-1

sigma: The linear variance of the corresponding sphere

N: The abundance of halos of that mass (Mpc^-3 h^3 comoving, per decade of mass)

Because numerical derivatives are involved, the value of delta_log_M affects the accuracy. Numerical experiments suggest that delta_log_M=0.1 gives more than enough accuracy, but you should check for your own use case.

Recommended m.f. for friends-of-friends linking length 0.2 particle sep.: z <~ 2 : bhattacharya z >~ 5 : reed_universal (no redshift dependence) : or reed_evolving (w/redshift dependence for additional accuracy)

pynbody.analysis.hmf.halo_bias(*fn_args, **fn_kwargs)[source]

Return the halo bias for the given halo mass.

Args:

context (SimSnap): The snapshot from which to pull the
cosmological context
M: float, unit or string describing the halo mass. If a
float, units are Msol h^-1.

Kwargs:

kern: The kernel function describing the halo bias (default
Cole-Kaiser).
pspec: A power spectrum object (which also defines the window
function); default is a WMAP7 cosmology calculated by CAMB, and a top hat window

delta_crit: The critical overdensity for collapse

Returns:

The halo bias (single float)

ionfrac

calculates ionization fractions - NEEDS DOCUMENTATION

pynbody.analysis.ionfrac.calculate(sim, ion='ovi', mode='old')[source]

calculate – documentation placeholder

luminosity

Calculates luminosities – NEEDS DOCUMENTATION

pynbody.analysis.luminosity.calc_mags(simstars, band='v')[source]

Calculating visible magnitudes

Using Padova Simple stellar populations (SSPs) from Girardi http://stev.oapd.inaf.it/cgi-bin/cmd Marigo+ (2008), Girardi+ (2010)

pynbody includes a grid of SSP luminosities for many bandpasses for various stellar ages and metallicities. This function linearly interpolates to the desired value and returns the value as a magnitude.

Usage:

>>> import pynbody
>>> pynbody.analysis.luminosity.calc_mags(h[1].s)

Optional keyword arguments:

band (default=’v’): Which observed bandpass magnitude in which
magnitude should be calculated
pynbody.analysis.luminosity.halo_mag(sim, band='v')[source]

Calculating halo magnitude

Calls pynbody.analysis.luminosity.calc_mags for ever star in passed in simulation, converts those magnitudes back to luminosities, adds those luminosities, then converts that luminosity back to magnitudes, which are returned.

Usage:

>>> import pynbody
>>> pynbody.analysis.luminosity.halo_mag(h[1].s)

Optional keyword arguments:

band (default=’v’): Which observed bandpass magnitude in which
magnitude should be calculated
pynbody.analysis.luminosity.halo_lum(sim, band='v')[source]

Calculating halo luminosiy

Calls pynbody.analysis.luminosity.calc_mags for every star in passed in simulation, converts those magnitudes back to luminosities, adds those luminosities, which are returned. Uses solar magnitudes from http://www.ucolick.org/~cnaw/sun.html.

Usage:

>>> import pynbody
>>> pynbody.analysis.luminosity.halo_mag(h[1].s)

Optional keyword arguments:

band (default=’v’): Which observed bandpass magnitude in which
magnitude should be calculated
pynbody.analysis.luminosity.half_light_r(sim, band='v')[source]

Calculate half light radius

Calculates entire luminosity of simulation, finds half that, sorts stars by distance from halo center, and finds out inside which radius the half luminosity is reached.

pkdgrav_cosmo

Cosmological module from PKDGRAV.

N.B. This code is being shared with skid and the I.C. generator.

NEEDS DOCUMENTATION

class pynbody.analysis.pkdgrav_cosmo.Cosmology(sim=None, H0=2.8944050182330705, Om=0.272, L=0.728, Ob=0.0456, Or=0.0, Quin=0.0, Ok=0.0)[source]

docs placeholder

profile

A set of classes and functions for making profiles of simulation properties.

class pynbody.analysis.profile.Profile(sim, load_from_file=False, ndim=2, type='lin', calc_x=None, weight_by='mass', **kwargs)[source]

A basic profile class for arbitrary profiles. Stores information about bins etc.

Made to work with the pynbody SimSnap instances. The constructor only generates the bins and figures out which particles belong in which bin. Profiles are generated lazy-loaded when a given property is requested.

Input:

sim : a simulation snapshot - this can be any subclass of SimSnap

Optional Keywords:

ndim (default = 2): specifies whether it’s a 2D or 3D profile - in the
2D case, the bins are generated in the xy plane
type (default = ‘lin’): specifies whether bins should be spaced linearly (‘lin’),
logarithmically (‘log’) or contain equal numbers of particles (‘equaln’)

min (default = min(x)): minimum value to consider

max (default = max(x)): maximum value to consider

nbins (default = 100): number of bins

bins : array like - predefined bin edges in units of binning quantity. If this
keyword is set, the values of the keywords type, nbins, min and max will be ignored
calc_x (default = None): function to use to calculate the value
for binning. If None it defaults to the radial distance from origin (in either 2 or 3 dimensions), ut you can specify this function to return any value you want for making profiles along arbitrary axes. Depening on your function, the units of certain profiles (such as density) might not make sense.
weight_by (default = ‘mass’): name of the array to use for weighting
averages across particles in each bin

Output:

a Profile object. To find out which profiles are available, use keys().

Implemented profile functions:

density : density

mass : mass in each bin

mass_enc : enclosed mass

fourier : provides fourier coefficients, amplitude and phase for
m=0 to m=6. To access the amplitude profile of m=2 mode, do p['fourier']['amp'][2,:]

dyntime : dynamical time

g_spherical: GM_enc/r^2

rotation_curve_spherical: rotation curve from vc = sqrt(GM/R) -
can be very wrong!

j_circ : angular momentum of particles on circular orbits

v_circ : circular velocity, aka rotation curve - calculated from
the midplane gravity, so this can be expensive

E_circ : energy of particles on circular orbits in the midplane

omega : circular orbital frequency

kappa : radial orbital frequency

beta : 3-D velocity anisotropy parameter

magnitudes : magnitudes in each bin - default band = ‘v’

sb : surface brightness - default band = ‘v’

Additional functions should use the profile_property to yield the desired profile.

Lazy-loading arrays:

The Profile class will automatically compute a mass-weighted profile for any lazy-loadable array of its parent SimSnap object.

Dispersions:

To obtain a dispersion profile, attach a _disp after the desired quantity name.

RMS:

The root-mean-square of a quantity can be obtained by using a _rms suffix

Derivatives:

To compute a derivative of a profile, prepend a d_ to the profile string, as in p['d_temp'] to get a temperature gradient.

Saving and loading previously generated profiles:

Use the write() function to write the current profiles with all the necessary information to a file. Initialize a profile with the load_from_file=True keyword to automatically load a previously saved profile. The filename is chosen automatically and corresponds to a hash generated from the positions of the particles used in the profile. This is to ensure that you are always looking at the same set of particles, centered in the same way. It also means you must use the same centering method if you want to reuse a saved profile.

Examples:

Density profile of the entire simulation:

>>> s = pynbody.load('mysim')
>>> import pynbody.profile as profile
>>> p = profile.Profile(s) # 2D profile of the whole simulation - note
                           # that this only makes the bins etc. but
                           # doesn't generate the density
>>> p['density'] # now we have a density profile
>>> p.keys()
['mass', 'n', 'density']
>>> p.families()
[<Family dm>, <Family star>, <Family gas>]

Density profile of the stars:

>>> ps = profile.Profile(s.s) # xy profile of the stars
>>> ps = profile.Profile(s.s, type='log') # same, but with log bins
>>> ps.families()
[<Family star>]
>>> import matplotlib.pyplot as plt
>>> plt.plot(ps['rbins'], ps['density'], 'o')
>>> plt.semilogy()

Metallicity profile of the gas in spherical shells (requires appropriate auxilliary files):

>>> pg = profile.Profile(s.g, ndim=3)
>>> pg['feh']
SimSnap: deriving array feh
TipsySnap: attempting to load auxiliary array 10/12M_hr.01000.FeMassFrac
SimSnap: deriving array hydrogen
SimSnap: deriving array hetot
SimArray([  1.83251940e-01,   1.48439968e-02,  -4.09390892e-01,
-1.82734736e+01])

Radial velocity dispersion profile and its gradient:

>>> ps = profile.Profile(s.s, max=15)
>>> ps['vr_disp']
SimSnap: deriving array vr
SimSnap: deriving array r
SimArray([ 118.80420996,  122.06102431,  131.13872886,  144.74447697,
35.89904165,   37.59565128,   35.21608633,   35.03373379], '1.01e+00 km s**-1')
>>> p['d_vr_disp']
SimArray([  21.71365764,   41.11802081,   75.61694836,  105.28110255,
2.664999  ,   -2.27668148,   -8.54033931,   -1.21577105], '1.01e+00 km s**-1 kpc**-1')

Using another quantity for binning:

>>> ps = profile.Profile(s.s, calc_x = lambda x: x.s['rform'])
keys()[source]

Returns a listing of available profile types

derivable_keys()[source]

Returns a list of possible profiles

families()[source]

Returns the family of particles used

create_particle_array(profile_name, particle_name=None, out_sim=None)[source]

Create a particle array with the results of the profile calculation.

After calling this function, sim[particle_name][i] == profile[profile_name][bin_in_which_particle_i_sits]

If particle_name is not specified, it defaults to the same as profile_name.

write()[source]

Writes all the vital information of the profile to a file.

To recover the profile, initialize a profile with the load_from_file=True keyword to automatically load a previously saved profile. The filename is chosen automatically and corresponds to a hash generated from the positions of the particles used in the profile. This is to ensure that you are always looking at the same set of particles, centered in the same way. It also means you must use the same centering method if you want to reuse a saved profile.

pynbody.analysis.profile.weight_fn(self, weight_by=None)[source]

Calculate mass in each bin

pynbody.analysis.profile.density(self)[source]

Generate a radial density profile for the current type of profile

pynbody.analysis.profile.fourier(self, delta_t='0.1 Myr', phi_bins=100)[source]

Generate a profile of fourier coefficients, amplitudes and phases

pynbody.analysis.profile.pattern_frequency(pro)[source]

Estimate the pattern speed from the m=2 Fourier mode

pynbody.analysis.profile.mass_enc(self)[source]

Generate the enclosed mass profile

pynbody.analysis.profile.density_enc(self)[source]

Generate the mean enclosed density profile

pynbody.analysis.profile.dyntime(self)[source]

The dynamical time of the bin, sqrt(R^3/2GM).

pynbody.analysis.profile.g_spherical(self)[source]

The naive gravitational acceleration assuming spherical symmetry = GM_enc/r^2

pynbody.analysis.profile.rotation_curve_spherical(self)[source]

The naive rotation curve assuming spherical symmetry: vc = sqrt(G M_enc/r)

pynbody.analysis.profile.j_circ(p)[source]

Angular momentum of particles on circular orbits.

pynbody.analysis.profile.v_circ(p, grav_sim=None)[source]

Circular velocity, i.e. rotation curve. Calculated by computing the gravity in the midplane - can be expensive

pynbody.analysis.profile.E_circ(p)[source]

Energy of particles on circular orbits.

pynbody.analysis.profile.pot(p)[source]

Calculates the potential in the midplane - can be expensive

pynbody.analysis.profile.omega(p)[source]

Circular frequency Omega = v_circ/radius (see Binney & Tremaine Sect. 3.2)

pynbody.analysis.profile.kappa(R dOmega^2/dR + 4 Omega^2) (see Binney & Tremaine Sect. 3.2)[source]
pynbody.analysis.profile.beta(p)[source]

3D Anisotropy parameter as defined in Binney and Tremiane

pynbody.analysis.profile.magnitudes(self, band='v')[source]

Calculate magnitudes in each bin

pynbody.analysis.profile.Q(self)[source]

Toomre Q parameter

pynbody.analysis.profile.X(self)[source]

X parameter defined as kappa^2*R/(2*pi*G*sigma*m) See Binney & Tremaine 2008, eq. 6.77

pynbody.analysis.profile.jtot(self)[source]

Magnitude of the total angular momentum

pynbody.analysis.profile.j_theta(self)[source]

Angle that the angular momentum vector of the bin makes with respect to the xy-plane.

pynbody.analysis.profile.j_phi(self)[source]

Angle that the angular momentum vector of the bin makes with the x-axis in the xy plane.

class pynbody.analysis.profile.InclinedProfile(sim, angle, load_from_file=False, ndim=2, type='lin', **kwargs)[source]

Creates a profile object to be used with a snapshot inclined by some known angle to the xy plane.

In addition to the SimSnap object, it also requires the angle to initialize.

Example:

>>> s = pynbody.load('sim')
>>> pynbody.analysis.angmom.faceon(s)
>>> s.rotate_x(60)
>>> p = pynbody.profile.InclinedProfile(s, 60)
class pynbody.analysis.profile.VerticalProfile(sim, rmin, rmax, zmax, load_from_file=False, ndim=3, type='lin', **kwargs)[source]

Creates a profile object that uses the absolute value of the z-coordinate for binning.

Input:

sim: snapshot to make a profile from

rmin: minimum radius for particle selection in kpc

rmax: maximum radius for particle selection in kpc

zmax: maximum height to consider in kpc

Optional Keywords:

ndim: if ndim=2, an edge-on projected profile is produced,
i.e. density is in units of mass/pc^2. If ndim=3 a volume profile is made, i.e. density is in units of mass/pc^3.
class pynbody.analysis.profile.QuantileProfile(sim, q=(0.16, 0.5, 0.84), weights=None, load_from_file=False, ndim=3, type='lin', **kwargs)[source]

Creates a profile object that returns the requested quantiles for a given array in a given bin. The quantiles may be mass weighted.

Input:

sim: snapshot to make a profile from

q (default: (0.16,0.5,0.84)):
The quantiles that will be returned. Default is median with 1-sigma on either side. q can be of arbitrary length allowing the user to select any quantiles they desire.
weights (default:None):
What should be used to weight the quantile. You will usually want to use particle mass: sim[‘mass’]. The default is to not weight by anything, weights=None.

Optional Keywords:

ndim: if ndim=2, an edge-on projected profile is produced,
i.e. density is in units of mass/pc^2. If ndim=3 a volume profile is made, i.e. density is in units of mass/pc^3.

ramses_util

Handy utilities for using RAMSES outputs in pynbody. For a complete demo on how to use RAMSES outputs with pynbody, have a look at the ipython notebook demo

Loading and centering

>>> s = pynbody.analysis.ramses_util.load_center('output_00101', align=False) # centered on halo 0 
>>> pynbody.analysis.ramses_util.hop_center(s,10) # centered on the halo 10

File Conversion

>>> pynbody.analysis.ramses_util.convert_to_tipsy_fullbox('output_00101') # will convert the whole output

Now you can run AHF or pkdgrav using the file named output_00101_fullbox.tipsy as an input or

>>> s_tipsy = pynbody.load('output_00101_fullbox.tipsy')

You can also just output a part of the simulation :

>>> s = pynbody.analysis.ramses_util.load_center('output_00101', align=False) # centered on halo 0 
>>> pynbody.analysis.ramses_util.convert_to_tipsy_simple('output_00101', file = pynbody.filt.Sphere('200 kpc')

Now we’ve got a file called output_00101.tipsy which holds only the 200 kpc sphere centered on halo 0.

Generating tform

A problem with RAMSES outputs in pynbody is that the tform array is in funny units that aren’t easily usable. To generate a new tform array (in Gyr) you can use the get_tform() defined here. It’s very easy:

>>> s = pynbody.load('output_00101')
>>> pynbody.analysis.ramses_util.get_tform(s)

This now generated a directory called birth in the parent directory of your output. It then calls the routine part2birth located in the RAMSES utils (see the bitbucket repository. get_tform() also deletes the previous tform array (not from disk, just from the currently loaded snapshot). The next time you call get_tform(), the data will be loaded from the disk and part2birth won’t need to be run again.

Feb 2016 - RS Note that in this version of ramses_util I have moved the ‘birth’ files into the output directory to which they pertain. The old code wasn’t placing them in a subdir of ‘birth’ in the top output and since I had to fix it I thought it better not to have a single directory with potentially 1000’s of files for all RAMSES outputs. The get_tform routine now creates the files under the approrpiate “output_” dir and reads them back from there.

pynbody.analysis.ramses_util.load_hop(s, hop='$HOME/ramses/trunk/ramses/utils/scripts/script_hop.sh')[source]

Loads the hop catalog for the given RAMSES snapshot. If the catalog doesn’t exist, it tries to run hop to create one via the ‘script_hop.sh’ script found in the RAMSES distribution. The hop output should be in a ‘hop’ directory in the base directory of the simulation.

Input:

s : loaded RAMSES snapshot

Optional Keywords:

hop : path to script_hop.sh

pynbody.analysis.ramses_util.hop_center(s, halo=0)[source]

Center the simulation snapshot on the specified halo using the halo data from hop.

Input:

s : RAMSES snapshot

Optional Keywords:

halo : halo ID to use for centering (default = 0)

pynbody.analysis.ramses_util.load_center(output, align=True, halo=0)[source]

Loads a RAMSES output and centers it on the desired halo. The hop center is used for an initial estimate, but for more precise centering, a shrinking-sphere center is calculated.

Inputs:

output : path to RAMSES output directory

Optional Keywords:

align : whether to align the snapshot based on the angular momentum in the central region (default = True)

halo : halo to center on (default = 0)

pynbody.analysis.ramses_util.convert_to_tipsy_simple(output, halo=0, filt=None)[source]

Convert RAMSES output to tipsy format readable by e.g. pkdgrav. This is a quick and dirty conversion, meant to be used for quick visualization or other simple post processing. Importantly, none of the cosmologically-relevant information is carried forward. For a more complete conversion for e.g. running through pkdgrav or Amiga Halo Finder, see convert_to_tipsy_fullbox().

The snapshot is put into units where G=1, time unit = 1 Gyr and mass unit = 2.222286e5 Msol.

Input:

output : path to RAMSES output directory

Optional Keywords:

filt : a filter to apply to the box before writing out the tipsy file

halo : which hop halo to center on – default = 0

pynbody.analysis.ramses_util.get_tipsy_units(sim)[source]

Returns snapshot sim units in the pkdgrav/gasoline unit system. This is probably not a function to be called by users, but it is used instead by other routines for file conversion.

Input:

sim: RAMSES simulation snapshot

Return values:

lenunit, massunit, timeunit : tuple specifying the units in kpc, Msol, and Gyr

pynbody.analysis.ramses_util.convert_to_tipsy_fullbox(output, write_param=True)[source]

Convert RAMSES file output to tipsy format readable by pkdgrav and Amiga Halo Finder. Does all unit conversions etc. into the pkdgrav unit system. Creates a file called output_fullbox.tipsy.

Input:

output: name of RAMSES output

Optional Keywords:

write_param: whether or not to write the parameter file (default = True)

pynbody.analysis.ramses_util.write_tipsy_param(sim, tipsyfile)[source]

Write a pkdgrav-readable parameter file for RAMSES snapshot sim with the prefix filename

pynbody.analysis.ramses_util.write_ahf_input(sim, tipsyfile)[source]

Write an input file that can be used by the Amiga Halo Finder with the corresponding tipsyfile which is the sim in tipsy format.

pynbody.analysis.ramses_util.get_tform(sim, part2birth_path='$HOME/ramses/trunk/ramses/utils/f90/part2birth')[source]

Use part2birth to calculate the formation time of stars in Gyr and replaces the original tform array.

Input:

sim: RAMSES snapshot

Optional Keywords:

part2birth_path: by default, this is
$HOME/ramses/trunk/ramses/utils/f90/part2birth, as specified in default_config.ini in your pynbody install directory. You can override this like so – make a file called ”.pynbodyrc” in your home directory, and include

[ramses]

ramses_utils = /path/to/your/ramses/utils/directory

interpolate

2D and 3D Interpolation routines written in cython

pynbody.analysis.interpolate.interpolate3d(x, y, z, x_vals, y_vals, z_vals, vals)[source]

Interpolate on a 3D regular grid. Yields results identical to scipy.interpolate.interpn.

pynbody.analysis.interpolate.interpolate2d(x, y, x_vals, y_vals, vals)[source]

Interpolate on a 2D regular grid. Yields results identical to scipy.interpolate.interpn.