Return the angular momentum vector of the specified snapshot.
The return units are [mass]*[dist]*[vel] as per the units of the snapshot.
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.
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.
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.
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.
A set of functions for common cosmological calculations.
returns: H(a) / H0 = [omegam/a**3 + (1-omegam)]**0.5
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.
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.
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
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
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.
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.
Tools for bulge/disk/halo decomposition
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_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).
Functions for calculating the gravitational potential and accelerations.
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):
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):
Calculates the midplane rotation curve at a specified set of points. Used primarily by pynbody.analysis.profile.v_circ()
Calculates the midplane potential at a specified set of points. Used primarily by pynbody.analysis.profile.pot()
Functions for dealing with and manipulating halos in simulations.
Return the center of mass velocity of the SimSnap
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:
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.
Determine the center of the halo by finding the shrink-sphere -center inside the specified distance of the potential minimum
Determine the center of mass based on specific particles.
Supply a list of indices using the ind keyword.
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.
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:
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.
Various halo mass function routines.
Calculate the linear density field correlation function.
Args:
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.
The Press-Schechter kernel used by halo_mass_function
Sheth & Tormen (1999) fit (see also Sheth Mo & Tormen 2001)
modified S-T fit by the G1 gaussian term and c
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
The Cole-Kaiser (1989) bias function. Also in Mo & White 1996.
The Sheth-Tormen (1999) bias function [eq 8]
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)
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)
calculates ionization fractions - NEEDS DOCUMENTATION
Calculates luminosities – NEEDS DOCUMENTATION
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
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
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
Cosmological module from PKDGRAV.
N.B. This code is being shared with skid and the I.C. generator.
NEEDS DOCUMENTATION
A set of classes and functions for making profiles of simulation properties.
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:
min (default = min(x)): minimum value to consider
max (default = max(x)): maximum value to consider
nbins (default = 100): number of bins
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
dyntime : dynamical time
g_spherical: GM_enc/r^2
j_circ : angular momentum of particles on circular orbits
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'])
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.
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.
Generate a radial density profile for the current type of profile
Generate a profile of fourier coefficients, amplitudes and phases
Estimate the pattern speed from the m=2 Fourier mode
The naive gravitational acceleration assuming spherical symmetry = GM_enc/r^2
The naive rotation curve assuming spherical symmetry: vc = sqrt(G M_enc/r)
Circular velocity, i.e. rotation curve. Calculated by computing the gravity in the midplane - can be expensive
Calculates the potential in the midplane - can be expensive
Circular frequency Omega = v_circ/radius (see Binney & Tremaine Sect. 3.2)
X parameter defined as kappa^2*R/(2*pi*G*sigma*m) See Binney & Tremaine 2008, eq. 6.77
Angle that the angular momentum vector of the bin makes with respect to the xy-plane.
Angle that the angular momentum vector of the bin makes with the x-axis in the xy plane.
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)
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:
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
Optional Keywords:
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
>>> 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
>>> 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.
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.
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
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)
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)
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
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
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)
Write a pkdgrav-readable parameter file for RAMSES snapshot sim with the prefix filename
Write an input file that can be used by the Amiga Halo Finder with the corresponding tipsyfile which is the sim in tipsy format.
Use part2birth to calculate the formation time of stars in Gyr and replaces the original tform array.
Input:
sim: RAMSES snapshot
Optional Keywords:
$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