Pynbody Demo

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.

In [1]:
%matplotlib inline
from matplotlib.pylab import *
rcParams['figure.figsize'] = (10,6)
rcParams['font.size'] = 18

Basic data loading/exploration

loading a simulation output using the pynbody.load() function, which tries to automatically determine which type of code output you have:

In [2]:
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:

In [3]:
s
Out[3]:
<SimSnap "/Users/rokstar/Downloads/testdata/g15784.lr.01024" len=1717156>
In [4]:
len(s)
Out[4]:
1717156
In [5]:
len(s.stars)
Out[5]:
265170

stars, gas, dark also available as s, g, d

In [6]:
len(s.g), len(s.gas), len(s.dark)
Out[6]:
(158755, 158755, 1293231)

the properties attribute of a SimSnap tells us some more basic info

In [7]:
s.properties
Out[7]:
{'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")}
In [8]:
s.properties['time'].in_units('Gyr')
Out[8]:
13.72831064485079

which quantities do we have available?

In [9]:
s.keys()
Out[9]:
[]

None! Because pynbody "lazy-loads" data... so lets see which data is actually on-disk:

In [10]:
s.loadable_keys()
Out[10]:
['phi',
 'FeMassFrac',
 'iord',
 'igasorder',
 'eps',
 'pos',
 'HI',
 'mass',
 'coolontime',
 'HeII',
 'vel',
 'OxMassFrac',
 'HeI']

Accessing data

to access any of these arrays or vectors, you access them like a python dictionary:

In [11]:
s['pos'].in_units('kpc')
TipsySnap: loading data from main file
Out[11]:
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...

In [12]:
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
Out[12]:
(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:

In [13]:
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

In [14]:
s['vt'],s['vr']
Out[14]:
(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

In [15]:
s.derivable_keys()[0:10]
Out[15]:
['nefe', 'hetot', 'mjeans', 'ofe', 'ljeans', 'ne', 'feh', 'c_s', 'mgfe', 'oxh']

Loading and centering/aligning halos

if we have information about substructure (i.e. from the Amiga Halo Finder) we can load in the halos:

In [16]:
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

In [17]:
pynbody.analysis.halo.center(h[1],mode='ssc')
Finding halo velocity center...
vcen= [   1.6295435   -58.21684952 -254.25581364]
Out[17]:
<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:

In [18]:
h1 = h[1]
h1
Out[18]:
<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:

In [19]:
len(h1.g), len(h1.s)
Out[19]:
(79060, 262178)
In [20]:
h[1]['mass'].sum().in_units('1e12 Msol')
Out[20]:
SimArray(1.696392414687047, '1.00e+12 Msol')

Rendering images

lets look at a picture of the dark matter distribution of a part of the whole simulation...

In [21]:
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...

In [22]:
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
Out[22]:
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:

In [23]:
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...
Out[23]:
<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:

In [24]:
s.rotate_x(90)
pynbody.plot.image(h[1].g,width='50 kpc',cmap=cm.Blues)
Rendering image on 4 threads...
Out[24]:
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:

In [25]:
pynbody.plot.image(h1.g,qty='temp',width='40 kpc', cmap=cm.YlOrRd)
s.rotate_x(-90)
Rendering image on 4 threads...
Out[25]:
<pynbody.transformation.GenericRotation at 0x114fcba10>

Profiles

In [26]:
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()
Out[26]:
<matplotlib.text.Text at 0x114c28490>

The Profile class has lots of built-in functionality... including calculating a rotation curve:

In [27]:
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)
Out[27]:
<matplotlib.legend.Legend at 0x10c30e4d0>

Profile can use any loaded (or loadable/derivable!) array:

In [28]:
plot(ps['rbins'],ps['age'])
xlabel('$R$ [kpc]')
ylabel('Age [Gyr]')
Out[28]:
<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:

In [29]:
plot(ps['rbins'],ps['d_v_circ'])
xlabel('$R$ [kpc]')
ylabel('d$v_{c}$/d$R$')
Profile: d_v_circ/dR
Out[29]:
<matplotlib.text.Text at 0x10c30cdd0>
In [30]:
plot(ps['rbins'],ps['vr_rms'])
xlabel('R [kpc]')
ylabel('$v_{rms}$')
Profile: auto-deriving vr_rms
Out[30]:
<matplotlib.text.Text at 0x11b5c6bd0>

Filters

You can define filters to easily select particles based on their quantities. For example, to get stars that are 1-2 Gyr old:

In [31]:
agefilt = pynbody.filt.BandPass('age', '1 Gyr','2 Gyr')
agefilt
Out[31]:
BandPass('age', 'Gyr', '2.00e+00 Gyr')
In [32]:
sphere = pynbody.filt.Sphere('50 kpc')
h1[sphere].s
Out[32]:
<SimSnap "/Users/rokstar/Downloads/testdata/g15784.lr.01024:halo_1:sphere::star" len=242858>
In [33]:
h1.s[sphere]['mass'].sum().in_units('1e12 Msol')
Out[33]:
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

In [34]:
h1.s[sphere&agefilt]['mass'].sum().in_units('1e12 Msol')
Out[34]:
SimArray(0.0023178475592408767, '1.00e+12 Msol')

Units

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:

In [35]:
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:

In [36]:
s['mass']*s['x']
Out[36]:
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.