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 xyplane 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 xyplane 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)[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 xz plane. This gives a sideon 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 xy plane. This gives a faceon 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 + (1omegam)]**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
(f, time)[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.
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):
The parameter eps (scalar, unit or array)
The array f[‘eps’]
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):
The parameter eps (scalar, unit or array)
The array f[‘eps’]
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_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=8, **kwargs)[source]¶ Return the center according to the shrinkingsphere 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')[source]¶ 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)[source]¶ Determine the center of the halo by finding the shrinksphere 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, precentre 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')[source]¶ 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
andlin
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'>)[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 PressSchechter 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 ST fit by the G1 gaussian term and c

pynbody.analysis.hmf.
f_reed_z_evo
(nu, neff, deltac=1.68647)[source]¶ modified ST 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 ColeKaiser (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 ShethTormen (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 friendsoffriends 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)[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
ColeKaiser).
 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)[source]¶ 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
luminosity¶
Calculates luminosities – NEEDS DOCUMENTATION

pynbody.analysis.luminosity.
use_custom_cmd
(path)[source]¶ 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/cgibin/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,x1]=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)[source]¶ Calculating visible magnitudes
Using Padova Simple stellar populations (SSPs) from Girardi http://stev.oapd.inaf.it/cgibin/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')[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', cylindrical=False)[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.
If cylindrical is True compute the half light radius as seen from the zaxis.
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
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)[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 lazyloaded 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 : 3D 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.
Lazyloading arrays:
The Profile class will automatically compute a massweighted profile for any lazyloadable array of its parent SimSnap object.
Dispersions:
To obtain a dispersion profile, attach a
_disp
after the desired quantity name.RMS:
The rootmeansquare of a quantity can be obtained by using a
_rms
suffixDerivatives:
To compute a derivative of a profile, prepend a
d_
to the profile string, as inp['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.
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

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.
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.
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.
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.
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.
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.
j_theta
(self)[source]¶ Angle that the angular momentum vector of the bin makes with respect to the xyplane.

pynbody.analysis.profile.
j_phi
(self)[source]¶ Angle that the angular momentum vector of the bin makes with the xaxis 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)
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)[source]¶ Creates a profile object that uses the absolute value of the zcoordinate 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 edgeon 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)[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 1sigma 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 edgeon 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)[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 cosmologicallyrelevant 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 pkdgravreadable 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/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/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