In [1]:
# Author: Volker Hoffmann <volmker@cheleb.net>, <volker@physik.uzh.ch>
# Date: 02 December 2014
In [2]:
%matplotlib inline
In [3]:
import matplotlib as mpl; mpl.rcParams['savefig.dpi'] = 144
import numpy as np
import matplotlib.pyplot as plt

Using Pynbody & Pkdgrav2 Data

In this notebook, we will use Pynbody [1] to load particle data from Pkdgrav2 [2]. Pkdgrav2 is a dark-matter only simulation, so we will only be using a small subset of Pynbody can do. Using different backends, we can, for example, load SPH particle data and use Pynbody to construct density, velocity, and pressure fields. We will expore this in a forthcoming article. Note that this notebook is somewhat compressed, and we refer to the official Pynbody introduction [3] for more details.

In [4]:
import pynbody

Load File Data

To load Pynbody simulation data, all we need to do is point to the simulation output. Note that Pynbody lazy-loads, i.e. it doesn't actually load the data unless they are specifically requested.

In [5]:
# Load Data, Show Some Stuff
s = pynbody.load('/zbox/data/dpotter/sr8.00100')
print s
INFO:pynbody.snapshot.tipsy:Loading /zbox/data/dpotter/sr8.00100
/zbox/opt/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/pynbody/snapshot/tipsy.py:130: RuntimeWarning: No readable param file in the run directory or parent directory: using defaults.
  "No readable param file in the run directory or parent directory: using defaults.", RuntimeWarning)
INFO:pynbody:Loading using backend <class 'pynbody.snapshot.tipsy.TipsySnap'>
INFO:pynbody.snapshot.tipsy:Loading /zbox/data/dpotter/sr8.00100

<SimSnap "/zbox/data/dpotter/sr8.00100" len=83840>

This simulation has only dark matter particles. Let's first check how many, and then check some properties, such as units.

In [6]:
print len(s)
print len(s.dark)
print s.properties
83840
83840
{'time': Unit("1.00e-05 kpc**3/2 G**-1/2 Msol**-1/2")}

We mentioned that Pynbody lazy-loads. So if we check what data are available to us with keys(), we will get an empty list because nothing has actually been loaded. We can use loadable_keys() to check what data we can load. We can also call all_keys() to check all the stuff that could be avaible through some post-processing routines.

In [7]:
print s.keys()
print s.loadable_keys()
print s.all_keys()
[]
['pos', 'phi', 'mass', 'vel', 'eps']
['nefe', 'hetot', 'mjeans', 'ofe', 'ljeans', 'ne', 'feh', 'c_s', 'mgfe', 'oxh', 'HeIII', 'sife', 'mjeans_turb', 'HII', 'ljeans_turb', 'hydrogen', 'c_s_turb', 'r_mag', 'i_mag', 'K_lum_den', 'vtheta', 'U_lum_den', 'vcxy', 'j2', 'I_lum_den', 'u_mag', 'cs', 'vr', 'vt', 'H_lum_den', 'V_lum_den', 'alt', 'i_lum_den', 'u', 'mu', 'az', 'I_mag', 'vrxy', 'u_lum_den', 'J_lum_den', 'k_mag', 'v_mag', 'U_mag', 'v_mean', 'h_lum_den', 'v_lum_den', 'theta', 'b_mag', 'zeldovich_offset', 'j_mag', 'B_lum_den', 'jz', 'j_lum_den', 'K_mag', 'V_mag', 'v2', 'J_mag', 'rho', 'H_mag', 'h_mag', 'B_mag', 'aform', 'k_lum_den', 'te', 'b_lum_den', 'temp', 'ke', 'age', 'j', 'smooth', 'rxy', 'p', 'r', 'R_lum_den', 'vphi', 'r_lum_den', 'R_mag', 'v_disp', 'pos', 'phi', 'mass', 'vel', 'eps']

Directly Accessing Particle Data

For example, to get a arrays with the mass, positions, and velcoity components of all dark matter particles, we can do the following. We can also access the gravitational softening of a particular particle ("eps") or the local gravitational potential ("phi").

In [8]:
print s.dark["mass"]
print s.dark["pos"]
print s.dark["vel"]
INFO:pynbody.snapshot.tipsy:Loading data from main file /zbox/data/dpotter/sr8.00100

[  2.53319754e-10   2.53319754e-10   2.53319754e-10 ...,   8.30078170e-06
   8.30078170e-06   8.30078170e-06]
[[-0.05146891 -0.02104948  0.02989316]
 [-0.05146838 -0.02105897  0.02991315]
 [-0.04741035 -0.02474593  0.02890897]
 ..., 
 [ 0.43263644 -0.47343358  0.39332488]
 [ 0.45491138 -0.46382084  0.3929171 ]
 [-0.49250564 -0.45727074  0.4013322 ]]
[[-0.05842622 -0.02269436  0.07257371]
 [-0.05656721 -0.02469489  0.07708147]
 [-0.05491808 -0.02834519  0.0724081 ]
 ..., 
 [ 0.03300809  0.08756977 -0.09323223]
 [ 0.00304865  0.10864146 -0.14196818]
 [ 0.03959431  0.10282841 -0.13245329]]

Custom Plots

With have access to the raw data, we can build any analysis and plots we want. For example, a histogram in mass could be made like follows. Note that this isn't a particularly good example though.

In [9]:
fig, ax = plt.subplots(1,1)

# Pretty Custom Color
teal = ( 119./255., \
         151./255., \
         161./255. )

# ax.hist(np.log10(i[a>0]), bins=128, log=True, color=c_amazing, histtype="stepfilled", alpha=0.5)
ax.hist(np.log10(s.dark["mass"]), bins=8, log=True, \
        color=teal, histtype="bar", \
        alpha=0.5)

ax.set_title("Mass Function")
ax.set_xlabel("Log10(Mass)")
ax.set_ylabel("dN/dM")
pass

Pynbody Plots

We can also use Pynbody's extensive plotting and analysis routines to make much prettier plots.

In [10]:
im = pynbody.plot.image(s.dark, width='1 kpc', \
                        units="Msol kpc^-2", cmap=mpl.cm.bone)
INFO:pynbody.snapshot:Deriving array rho
INFO:pynbody.sph:Building 4 trees with leafsize=16
INFO:pynbody.sph:Tree build done in  0.04 s
INFO:pynbody.snapshot:Deriving array smooth
INFO:pynbody.sph:Smoothing with 32 nearest neighbours
INFO:pynbody.sph:Smoothing done in 0.172s
INFO:pynbody.sph:Calculating SPH density
INFO:pynbody.sph:Density calculation done in 0.0725 s
INFO:pynbody.sph:Rendering image on 4 threads...