Analysis

angmom

pynbody.analysis.angmom.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.

pynbody.analysis.angmom.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.

pynbody.analysis.angmom.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.

pynbody.analysis.angmom.sideon(h, vec_to_xform=<function 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.

pynbody.analysis.angmom.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.

cosmology

A set of functions for common cosmological calculations.

pynbody.analysis.cosmology.hzoverh0(a, omegam0)

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

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

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')

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')

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(f, time)

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)

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)

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)

Calculate the Hubble parameter of the universe in snapshot f

pynbody.analysis.cosmology.add_hubble(f)

Add the hubble flow to velocities in snapshot f

pynbody.analysis.cosmology.comoving_to_physical(ar)

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')

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)

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)

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)

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)

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)

Return the centre of mass of the SimSnap

pynbody.analysis.halo.center_of_mass_velocity(sim)

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=8, **kwargs)

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, rho_def='matter')

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

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

Input:

sim : a simulation snapshot - this can be any subclass of SimSnap, especially a halo.

Optional Keywords:

cen (default=None): Provides the halo center. If None, assumes that the snapshot is already centered.

*rmax (default=None): Maximum radius to start the search. If None, inferred from the halo particle positions.

*overden (default=178): Overdensity corresponding to the required halo boundary definition. 178 is the virial criterion for spherical collapse in an EdS Universe. Common possible values are 200, 500 etc

*rho_def (default=’matter’): Physical density used to define the overdensity. Default is the matter density at the redshift of the simulation. An other choice is “critical” for the critical density at this redshift.

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

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)

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)

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)

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.

pynbody.analysis.halo.halo_shape(sim, N=100, rin=None, rout=None, bins='equal')

Returns radii in units of sim['pos'], axis ratios b/a and c/a, the alignment angle of axis a in radians, and the rotation matrix for homeoidal shells over a range of N halo radii.

Keyword arguments:

N (100): The number of homeoidal shells to consider. Shells with few particles will take longer to fit.

rin (None): The minimum radial bin in units of sim['pos']. Note that this applies to axis a, so particles within this radius may still be included within homeoidal shells. By default this is taken as rout/1000.

rout (None): The maximum radial bin in units of sim['pos']. By default this is taken as the largest radial value in the halo particle distribution.

bins (equal): The spacing scheme for the homeoidal shell bins. equal initialises radial bins with equal numbers of particles, with the exception of the final bin which will accomodate remainders. This number is not necessarily maintained during fitting. log and lin initialise bins with logarithmic and linear radial spacing.

Halo must be in a centered frame. Caution is advised when assigning large number of bins and radial ranges with many particles, as the algorithm becomes very slow.

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'>)

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)

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)

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

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

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

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

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)

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)

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)

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(context, M, kern=<function cole_kaiser_bias>, pspec=<class 'pynbody.analysis.hmf.PowerSpectrumCAMB'>, delta_crit=1.686)

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)

pynbody.analysis.hmf.simulation_halo_mass_function(snapshot, log_M_min=8.0, log_M_max=15.0, delta_log_M=0.1, masses=None, mass_def='halo_finder', calculate_err=True, subsample_catalogue=None)

Construct the halo mass function from a halo catalogue by binning haloes in mass and counting them. This allows direct plotting of the simulation HMF against the theory HMF.

Args:

snapshot (SimSnap): The snapshot from which to calculate halo masses

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

masses: Provide array of halo masses in Msol h*-1. If none, this is calculated from the snapshot.

*mass_def: Definition of the mass of a halo. Possible extensions could be M200_crit, M200_matter, Mvir etc…

*err: If true, calculates error bars according to Poisson process.

*subsample_catalogue: Sample the halo catalogue every int value. Mostly for speed and debugging issues. WARNING: If your halo catalogue is ordered, this method does not take responsability in randomising the catalogue so you might not recover the true HMF.

Returns:

The binned number density of haloes in this snapshot in comoving Mpc**-3 h**3 per decade of mass

ionfrac

calculates ionization fractions - NEEDS DOCUMENTATION

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

calculate – documentation placeholder

luminosity

Calculates luminosities – NEEDS DOCUMENTATION

pynbody.analysis.luminosity.use_custom_cmd(path)

Use a custom set of stellar populations to calculate magnitudes.

The path is to a numpy archive with a suitable grid of ages/metallicities and corresponding magnitudes.

The following script from Stephanie De Beer should help you make a suitable file starting from ugriz SSPs downloaded from http://stev.oapd.inaf.it/cgi-bin/cmd.

import numpy as np

metals = np.linspace(0.002, 0.05, 25) ages = np.logspace(5.67, 10.13, 25)

mags_bol = np.zeros((len(metals),len(ages))) mags_u = np.zeros((len(metals),len(ages))) mags_g = np.zeros((len(metals),len(ages))) mags_r = np.zeros((len(metals),len(ages))) mags_i = np.zeros((len(metals),len(ages))) mags_z = np.zeros((len(metals),len(ages)))

bands = [‘bol’, ‘u’, ‘g’, ‘r’, ‘i’, ‘z’] k=2 for b in bands:

for x in range(1,26):

with open(‘/users/sdebeer/render_stuff/PGSP_files/’+str(x)+’_output.txt’, ‘r’) as f:

output = f.readlines()

i = 0 for line in output:

magnitudes = line.split() if magnitudes[0]==’#’:

continue

vars()[‘mags_’+b][i,x-1]=magnitudes[k] i +=1

k+=1

np.savez(‘my_cmd.npz’, ages=ages, metals=metals, bol=mags_bol, u=mags_u, g=mags_g, r=mags_r, i=mags_i, z=mags_z)

pynbody.analysis.luminosity.calc_mags(simstars, band='v', cmd_path=None)

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

path (default=None): Path to the CMD grid. If None, use the

default or a path specified by use_custom_cmd. For more information about generating a custom CMD grid, see use_custom_cmd.

pynbody.analysis.luminosity.halo_mag(sim, band='v')

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')

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', cylindrical=False)

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.

If cylindrical is True compute the half light radius as seen from the z-axis.

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)

docs placeholder

Methods

ComoveKickFac(dTime, dDelta)
Exp2Hub(dExp)
ExpDot2(dExp)

ComoveDriftFac
ComoveDriftInt
ComoveKickInt
ComoveLookbackTime2Exp
CosmoTint
Exp2Om
Exp2Time
GrowthFac
GrowthFacDeriv
GrowthFacDot
Time2Exp
Time2Hub

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)

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’)

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

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

nbins (default = 100): number of bins

binsarray like - predefined bin edges in units of binning quantity. If this

keyword is set, the values of the keywords type, nbins, rmin and rmax 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

fourierprovides 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_circcircular 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.

Methods

create_particle_array(profile_name[, …])	Create a particle array with the results of the profile calculation.
derivable_keys()	Returns a list of possible profiles
families()	Returns the family of particles used
keys()	Returns a listing of available profile types
write()	Writes all the vital information of the profile to a file.

profile_property

keys()

Returns a listing of available profile types

derivable_keys()

Returns a list of possible profiles

families()

Returns the family of particles used

create_particle_array(profile_name, particle_name=None, out_sim=None)

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()

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)

Calculate mass in each bin

pynbody.analysis.profile.density(self)

Generate a radial density profile for the current type of profile

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

Generate a profile of fourier coefficients, amplitudes and phases

pynbody.analysis.profile.pattern_frequency(pro)

Estimate the pattern speed from the m=2 Fourier mode

pynbody.analysis.profile.mass_enc(self)

Generate the enclosed mass profile

pynbody.analysis.profile.density_enc(self)

Generate the mean enclosed density profile

pynbody.analysis.profile.dyntime(self)

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

pynbody.analysis.profile.g_spherical(self)

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

pynbody.analysis.profile.rotation_curve_spherical(self)

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

pynbody.analysis.profile.j_circ(p)

Angular momentum of particles on circular orbits.

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

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

pynbody.analysis.profile.E_circ(p)

Energy of particles on circular orbits.

pynbody.analysis.profile.pot(p)

Calculates the potential in the midplane - can be expensive

pynbody.analysis.profile.omega(p)

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)

pynbody.analysis.profile.beta(p)

3D Anisotropy parameter as defined in Binney and Tremiane

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

Calculate magnitudes in each bin

pynbody.analysis.profile.Q(self)

Toomre Q parameter

pynbody.analysis.profile.X(self)

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

pynbody.analysis.profile.jtot(self)

Magnitude of the total angular momentum

pynbody.analysis.profile.j_theta(self)

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

pynbody.analysis.profile.j_phi(self)

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)

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)

Methods

create_particle_array(profile_name[, …])	Create a particle array with the results of the profile calculation.
derivable_keys()	Returns a list of possible profiles
families()	Returns the family of particles used
keys()	Returns a listing of available profile types
write()	Writes all the vital information of the profile to a file.

profile_property

class pynbody.analysis.profile.VerticalProfile(sim, rmin, rmax, zmax, load_from_file=False, ndim=3, type='lin', **kwargs)

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.

Methods

create_particle_array(profile_name[, …])	Create a particle array with the results of the profile calculation.
derivable_keys()	Returns a list of possible profiles
families()	Returns the family of particles used
keys()	Returns a listing of available profile types
write()	Writes all the vital information of the profile to a file.

profile_property

class pynbody.analysis.profile.QuantileProfile(sim, q=0.16, 0.5, 0.84, weights=None, load_from_file=False, ndim=3, type='lin', **kwargs)

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.

Methods

create_particle_array(profile_name[, …])	Create a particle array with the results of the profile calculation.
derivable_keys()	Returns a list of possible profiles
families()	Returns the family of particles used
keys()	Returns a listing of available profile types
write()	Writes all the vital information of the profile to a file.

profile_property

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

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.convert_to_tipsy_simple(output, halo=0, filt=None)

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)

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)

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)

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

pynbody.analysis.ramses_util.write_ahf_input(sim, tipsyfile)

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/utils/f90/part2birth')

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/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)

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

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

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