This notebook highlights some of the basic pynbody
functionality to get you on the way to analyzing your own simulations. The pynbody
webpage is here and you can find the documentation with additional tutorials here. If you find that things are broken please please please let us know by submitting a bug report on github. If you want to run the notebook on your own machine, you need to get the testdata tarball and change the path in the load
function below accordingly.
%matplotlib inline
from matplotlib.pylab import *
rcParams['figure.figsize'] = (10,6)
rcParams['font.size'] = 18
loading a simulation output using the pynbody.load()
function, which tries to automatically determine which type of code output you have:
import pynbody
s = pynbody.load('/Users/rokstar/Downloads/testdata/g15784.lr.01024.gz')
Loading using backend <class 'pynbody.tipsy.TipsySnap'> TipsySnap: loading /Users/rokstar/Downloads/testdata/g15784.lr.01024.gz
a quick look at some basic information about the run we just opened:
s
<SimSnap "/Users/rokstar/Downloads/testdata/g15784.lr.01024" len=1717156>
len(s)
1717156
len(s.stars)
265170
stars
, gas
, dark
also available as s
, g
, d
len(s.g), len(s.gas), len(s.dark)
(158755, 158755, 1293231)
the properties
attribute of a SimSnap
tells us some more basic info
s.properties
{'a': 0.9999999999999911, 'boxsize': Unit("6.85e+04 kpc a"), 'h': 0.7301145776501103, 'omegaL0': 0.76, 'omegaM0': 0.24, 'time': Unit("1.40e+01 s kpc km**-1")}
s.properties['time'].in_units('Gyr')
13.72831064485079
which quantities do we have available?
s.keys()
[]
None! Because pynbody "lazy-loads" data... so lets see which data is actually on-disk:
s.loadable_keys()
['phi', 'FeMassFrac', 'iord', 'igasorder', 'eps', 'pos', 'HI', 'mass', 'coolontime', 'HeII', 'vel', 'OxMassFrac', 'HeI']
to access any of these arrays or vectors, you access them like a python dictionary:
s['pos'].in_units('kpc')
TipsySnap: loading data from main file
SimArray([[ 733.51516346, -2479.31184326, -11394.52608302], [ 781.21390409, -2209.05638147, -11392.16434933], [ 794.4949563 , -2232.21137376, -11432.30055384], ..., [ 1672.78285106, -2332.06734179, -8386.4632496 ], [ 1675.12072757, -2335.55559387, -8386.09837501], [ 1682.45100698, -2338.79940549, -8386.07898308]], 'kpc')
Note that each array has units attached...
s['mass'], s['vel'], s.g['HeI'], s.s['posform']
Reading /Users/rokstar/Downloads/testdata/g15784.lr.starlog Bridging starlog into SimSnap
TipsySnap: attempting to load auxiliary array /Users/rokstar/Downloads/testdata/g15784.lr.01024.HeI TipsySnap: attempting to load auxiliary array
/Users/rokstar/Downloads/testdata/g15784.lr.01024.iord
(SimArray([ 3.19310585e-11, 3.19310585e-11, 3.19310585e-11, ..., 1.24175999e-11, 1.24175999e-11, 1.24175999e-11], '4.75e+16 Msol'), SimArray([[ 0.01673503, -0.01547551, -0.16838804], [ 0.05142973, -0.00647544, -0.16129912], [ 0.02738659, -0.00184103, -0.1736342 ], ..., [-0.13633673, -0.13716312, -0.12805983], [-0.13267632, -0.09228273, -0.1247804 ], [ 0.06591902, 0.13406622, -0.15172358]], '1.73e+03 km a s**-1'), SimArray([ 3.51669648e-12, 1.60111191e-09, 1.88118882e-10, ..., 1.94256286e-12, 2.30497435e-12, 1.81767830e-11], dtype=float32, '1.00e+00'), SimArray([[ 0.0178572 , -0.0131114 , -0.00092187], [ 0.01790398, -0.01335163, -0.00124265], [ 0.01781229, -0.01359182, -0.00148088], ..., [ 0.02442261, -0.03404816, -0.12244228], [ 0.02445674, -0.03409909, -0.12243695], [ 0.02456377, -0.03414645, -0.12243667]], '6.85e+04 kpc aform'))
by default everything is in system units, but most of the time thinking in physical units is easier:
s.physical_units()
We have defined many useful quantities that are automatically calculated for you. For example, the radial and tangential velocities are simply obtained by
s['vt'],s['vr']
(SimArray([ 37.28290565, 81.78390831, 60.75947459, ..., 330.09191502, 288.75911925, 303.1394157 ], 'km s**-1'), SimArray([ 291.20028176, 281.05583586, 297.57331009, ..., 227.22613725, 202.62874018, 208.42590924], 'km s**-1'))
you can get a list of all available derivable quantities
s.derivable_keys()[0:10]
['nefe', 'hetot', 'mjeans', 'ofe', 'ljeans', 'ne', 'feh', 'c_s', 'mgfe', 'oxh']
if we have information about substructure (i.e. from the Amiga Halo Finder) we can load in the halos:
h = s.halos()
AHFCatalogue: loading particles... halos... done!
We are not really near the main halo, so lets center the simulation there using the "shrinking sphere" method
pynbody.analysis.halo.center(h[1],mode='ssc')
Finding halo velocity center... vcen= [ 1.6295435 -58.21684952 -254.25581364]
<pynbody.transformation.GenericTranslate at 0x1090d4f50>
we're going to work with the main halo a lot, so it's useful to save a reference to it in a variable to save some typing:
h1 = h[1]
h1
<SimSnap "/Users/rokstar/Downloads/testdata/g15784.lr.01024:halo_1" len=502300>
this new object is a SubSnap
and behaves exactly like any other SimSnap
:
len(h1.g), len(h1.s)
(79060, 262178)
h[1]['mass'].sum().in_units('1e12 Msol')
SimArray(1.696392414687047, '1.00e+12 Msol')
lets look at a picture of the dark matter distribution of a part of the whole simulation...
im = pynbody.plot.image(s.d,width='50 Mpc', units='Msol kpc^-2', cmap=cm.Greys)
Building 4 trees with leafsize=16 Tree build done in 1.55 s Calculating SPH density Smoothing with 32 nearest neighbours Smoothing done in 3.36s Density calculation done in 1.48 s
Rendering image on 4 threads...
pynbody.plot.image(h[1].g,width='200 kpc', av_z=True, cmap=cm.Blues)
Rendering image on 4 threads... Rendering image on 4 threads...
Building 4 trees with leafsize=16 Tree build done in 0.177 s Smoothing with 32 nearest neighbours Smoothing done in 0.451s
SimArray([[ 773.51574707, 775.52825928, 777.54034424, ..., 619.25726318, 617.1003418 , 614.94384766], [ 776.91687012, 778.97802734, 781.03839111, ..., 621.38391113, 619.21942139, 617.05560303], [ 780.31646729, 782.42590332, 784.53491211, ..., 623.51269531, 621.3404541 , 619.16906738], ..., [ 741.20367432, 743.32147217, 745.43774414, ..., 882.82843018, 879.76422119, 876.706604 ], [ 737.50848389, 739.64691162, 741.78356934, ..., 880.90118408, 877.8538208 , 874.81311035], [ 733.80169678, 735.96105957, 738.11846924, ..., 878.97247314, 875.94213867, 872.9185791 ]], dtype=float32, 'Msol kpc**-3')
Typically we want to make sure the system is aligned in some predictable way, so pynbody allows us to easily make the disk angular momentum axis the z-axis:
pynbody.analysis.angmom.faceon(h[1],cen=(0,0,0))
Finding halo velocity center... vcen= [ -1.28863327e-14 4.35041217e-13 1.22950847e-12] Calculating angular momentum vector... Transforming simulation...
<pynbody.transformation.GenericRotation at 0x114f23c10>
we use cen=(0,0,0)
because we have already centered it in the previous step... if you don't specify a center, it is first centered and then aligned. Now lets look at a slice of the gas disk edge-on:
s.rotate_x(90)
pynbody.plot.image(h[1].g,width='50 kpc',cmap=cm.Blues)
Rendering image on 4 threads...
SimArray([[ 20220.50976562, 20200.53125 , 20180.5546875 , ..., 19654.78320312, 19650.26367188, 19645.74414062], [ 20184.61523438, 20167.7109375 , 20150.80664062, ..., 19755.99609375, 19749.9296875 , 19743.86328125], [ 20148.72265625, 20134.88867188, 20121.05664062, ..., 19857.20703125, 19849.59570312, 19841.98632812], ..., [ 24377.94921875, 24315.953125 , 24253.95507812, ..., 19957.59765625, 19975.203125 , 19992.80664062], [ 24303.73828125, 24237.01171875, 24170.28515625, ..., 19852.88867188, 19864.32617188, 19875.76367188], [ 24229.52539062, 24158.0703125 , 24086.6171875 , ..., 19748.1796875 , 19753.44921875, 19758.72070312]], dtype=float32, 'Msol kpc**-3')
By default, rho
(i.e. density) is used to render the image, but you can any valid array is fair game:
pynbody.plot.image(h1.g,qty='temp',width='40 kpc', cmap=cm.YlOrRd)
s.rotate_x(-90)
Rendering image on 4 threads...
<pynbody.transformation.GenericRotation at 0x114fcba10>
p = pynbody.analysis.profile.Profile(h1,max='200 kpc',min='.1 kpc',type='log',nbins=100)
ps = pynbody.analysis.profile.Profile(h1.s,max='50 kpc',nbins=50)
pg = pynbody.analysis.profile.Profile(h1.g,max='50 kpc',nbins=50)
plot(p['rbins'], p['density'], label = 'all')
plot(ps['rbins'], ps['density'], label = 'stars')
plot(pg['rbins'], pg['density'], label = 'gas')
semilogy()
xlim(0,50)
legend()
xlabel('$R$ [kpc]')
ylabel('$\Sigma$ [M$_{\odot}$ kpc$^{-2}$]')
Profile: density() Profile: mass() Profile: density() Profile: mass() Profile: density() Profile: mass()
<matplotlib.text.Text at 0x114c28490>
The Profile
class has lots of built-in functionality... including calculating a rotation curve:
plot(p['rbins'], p['v_circ'], label = 'all')
plot(ps['rbins'], ps['v_circ'], label = 'stars')
plot(pg['rbins'], pg['v_circ'], label = 'gas')
xlim(0,50)
xlabel('$R$ [kpc]')
ylabel('$V_{circ}$')
legend()
Profile: v_circ() -- warning, disk must be in the x-y plane Rotation curve calculated in 4.13 s Profile: v_circ() -- warning, disk must be in the x-y plane Rotation curve calculated in 1.33 s Profile: v_circ() -- warning, disk must be in the x-y plane Rotation curve calculated in 0.45 s
/Users/rokstar/python/pynbody/gravity/calc.py:105: RuntimeWarning: OpenMP gravity nt able to load -- using single cpu warnings.warn("OpenMP gravity nt able to load -- using single cpu", RuntimeWarning)
<matplotlib.legend.Legend at 0x10c30e4d0>
Profile
can use any loaded (or loadable/derivable!) array:
plot(ps['rbins'],ps['age'])
xlabel('$R$ [kpc]')
ylabel('Age [Gyr]')
<matplotlib.text.Text at 0x10c30cd10>
By default, Profile
calculates the means of quantities, but it trivially calculate derivatives and root-mean-square values as well:
plot(ps['rbins'],ps['d_v_circ'])
xlabel('$R$ [kpc]')
ylabel('d$v_{c}$/d$R$')
Profile: d_v_circ/dR
<matplotlib.text.Text at 0x10c30cdd0>
plot(ps['rbins'],ps['vr_rms'])
xlabel('R [kpc]')
ylabel('$v_{rms}$')
Profile: auto-deriving vr_rms
<matplotlib.text.Text at 0x11b5c6bd0>
You can define filters to easily select particles based on their quantities. For example, to get stars that are 1-2 Gyr old:
agefilt = pynbody.filt.BandPass('age', '1 Gyr','2 Gyr')
agefilt
BandPass('age', 'Gyr', '2.00e+00 Gyr')
sphere = pynbody.filt.Sphere('50 kpc')
h1[sphere].s
<SimSnap "/Users/rokstar/Downloads/testdata/g15784.lr.01024:halo_1:sphere::star" len=242858>
h1.s[sphere]['mass'].sum().in_units('1e12 Msol')
SimArray(0.10313311567845437, '1.00e+12 Msol')
you can combine filters by using simple logic operators -- so to get the total mass of all 1-2 Gyr old stars within 50 kpc, you can do
h1.s[sphere&agefilt]['mass'].sum().in_units('1e12 Msol')
SimArray(0.0023178475592408767, '1.00e+12 Msol')
You've noticed by now that we've been asking for things in various kinds of units. Pynbody
has a unit-handling system built in that automatically ensures that arrays in different units are compatible and performs the necessary conversions for calculations. Adding mass and position together makes no sense, and pynbody knows this:
s['mass'] + s['x']
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-35-3fb92e8fe716> in <module>() ----> 1 s['mass'] + s['x'] /Users/rokstar/python/pynbody/array.pyc in __add__(self, x) 421 return x+self 422 else : --> 423 return self._generic_add(x) 424 425 def __sub__(self, x) : /Users/rokstar/python/pynbody/array.pyc in _generic_add(self, x, add_op) 393 **context) 394 except units.UnitsException : --> 395 raise ValueError("Incompatible physical dimensions %r and %r, context %r"%(str(self.units),str(x.units), str(self.conversion_context()))) 396 397 ValueError: Incompatible physical dimensions 'Msol' and 'kpc', context "{'a': 0.9999999999999911, 'h': 0.7301145776501103}"
but, it handles this just fine and returns an array with the right units:
s['mass']*s['x']
SimArray([ -1.11405446e+09, -1.04208223e+09, -1.01781848e+09, ..., -1.34612686e+06, 1.56832943e+04, 4.33084960e+06], 'Msol kpc')
This concludes the basic overview of what pynbody
can do -- for more tutorials and in-depth explanation of the inner-workings of the framework, head over to the walkthroughs and cookbook examples. If you're stuck implementing a particular analysis workflow, stop by our Github issues and post your question, we'll be happy to help.