This notebook highlights some of the basic pynbody
functionality to get you on the way to analyzing your own simulations. It doesn't assume any python knowledge, though familiarity with basic concepts of the language will certainly help.
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 point the load
function below to whatever output you wish to use.
In what follows, we demonstrate the usage of pynbody with Ramses simulation outputs. Support for this particular format is still under active development and although it is becoming increasingly robust, you may encounter strange behavior.
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.
Note that the pynbody.ramses.multiprocess_num
and pynbody.config['number_of_threads']
define how many threads the Ramses reader and the pynbody routines should use. Change these to match your system.
import pynbody
pynbody.ramses.multiprocess_num = 16
pynbody.config['number_of_threads'] = 16
s = pynbody.load('/zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101')
Loading using backend <class 'pynbody.ramses.RamsesSnap'>
A quick look at some basic information about the run we just opened:
s
<SimSnap "/zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101" len=23196808>
len(s)
23196808
len(s.stars)
1243472
stars
, gas
, dark
also available as s
, g
, d
len(s.g), len(s.gas)
(14817745, 14817745)
The properties
attribute of a SimSnap
is a dictionary that tells us some more basic info about the simulation:
s.properties
{'a': 1.00001844594404, 'boxsize': Unit("6.16e+26 cm"), 'h': 0.7040000152587891, 'omegaL0': 0.727999985218048, 'omegaM0': 0.272000014781952, 'time': Unit("1.38e+01 Gyr")}
s.properties['time'].in_units('Gyr')
13.760107374692666
So far, we have only opened the simulation and read its metadata -- none of the actual cell or particle data has been read off the disk. To access any of these arrays or vectors, you access them like a python dictionary:
which quantities do we have available?
s.keys()
[]
None! Because pynbody "lazy-loads" data... we need to try to access an array to read the data off the disk:
s['pos']
RamsesSnap: loading hydro files 135791113151719212325273129/home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) /home/ics/roskar/python/pynbody/ramses.py:258: RuntimeWarning: More hydro variables (8) are in this RAMSES dump than are defined in config.ini (6) warnings.warn("More hydro variables (%d) are in this RAMSES dump than are defined in config.ini (%d)"%(nvar_file,nvar), RuntimeWarning) 2468101214161820222426283230 61 63 37 57 41 39 43 33 35 45 47 53 51 49 55 59 62 64 38 58 42 40 44 34 36 46 48 54 52 50 56 60 91 95 69 93 87 85 75 65 67 79 77 83 81 73 71 89 92 96 70 94 88 86 76 66 68 80 78 84 82 74 72 90 125 127 103 121 113 117 111 97 99 107 115 105 119 109 101 123 126 128 104 122 114 118 112 98 100 108 116 106 120 110 102 124done
SimArray([[ 0.01651346, 0.02221355, 0.05327416], [ 0.16698388, 0.01523324, 0.05411253], [ 0.06062804, 0.05369276, 0.12354728], ..., [ 0.50014782, 0.4998312 , 0.49989605], [ 0.50015163, 0.49983883, 0.49989605], [ 0.50015163, 0.4998312 , 0.49989605]], '6.16e+26 cm')
Note that each array is type SimArray
(an extended numpy
array) and has units attached.
s.keys()
['pos', 'vx', 'vy', 'y', 'x', 'vel', 'z', 'vz']
Note that not all particle families have the same arrays:
s.s.keys(), s.g.keys()
(['pos', 'vx', 'vy', 'y', 'x', 'vel', 'z', 'vz'], ['p', 'metal', 'smooth', 'pos', 'vx', 'rho', 'vy', 'y', 'x', 'vel', 'z', 'vz'])
s['mass'], s['vel']
(SimArray([ 3.97948659e-07, 3.97948659e-07, 3.97948659e-07, ..., 9.86356141e-18, 9.22962765e-18, 9.46693561e-18], '5.92e+50 g'), SimArray([[ 0.00728201, -0.02102769, 0.02122281], [-0.00035387, -0.02282961, 0.02129209], [ 0.0017414 , -0.01511998, 0.00469109], ..., [-0.01208301, -0.00138862, -0.00770179], [-0.01208318, -0.00153632, -0.00770602], [-0.01216255, -0.0014085 , -0.00769588]], '1.41e+09 cm s**-1'))
By default everything is in system units, but you can request an array in different units:
s['mass'].in_units('m_p'), s['vel'].in_units('Mpc yr^-1')
(SimArray([ 1.40943029e+68, 1.40943029e+68, 1.40943029e+68, ..., 3.49341602e+57, 3.26889323e+57, 3.35294151e+57], 'm_p'), SimArray([[ 1.04855120e-10, -3.02781894e-10, 3.05591523e-10], [ -5.09548653e-12, -3.28728237e-10, 3.06589151e-10], [ 2.50748181e-11, -2.17715615e-10, 6.75479899e-11], ..., [ -1.73985636e-10, -1.99949828e-11, -1.10899682e-10], [ -1.73988205e-10, -2.21217356e-11, -1.10960552e-10], [ -1.75131069e-10, -2.02813170e-11, -1.10814482e-10]], 'Mpc yr**-1'))
in_units
produces a "view" of the array in the desired units, but the data is not altered in memory. You can permanently alter it with
s['mass'].convert_units('m_p')
s['mass']
SimArray([ 1.40943029e+68, 1.40943029e+68, 1.40943029e+68, ..., 3.49341602e+57, 3.26889323e+57, 3.35294151e+57], 'm_p')
Or you can ask for everything to be converted to "physical" units (there is a similar convenience function to convert back to original units also):
s.physical_units()
SimSnap: converting pos units from 6.16e+26 cm to kpc SimSnap: converting vel units from 1.41e+09 cm s**-1 to km s**-1 SimSnap: converting mass units from m_p to Msol SimSnap: converting p units from 5.02e-12 g cm**-1 s**-2 to km**2 Msol s**-2 kpc**-3 SimSnap: converting smooth units from 6.16e+26 cm to kpc SimSnap: converting rho units from 2.53e-30 g cm**-3 to Msol kpc**-3
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']
SimSnap: deriving array vt SimSnap: deriving array v2 SimSnap: deriving array vr SimSnap: deriving array r
(SimArray([ 392.0037547 , 435.51170488, 223.91531001, ..., 107.00691314, 105.46517315, 107.53279984], 'km s**-1'), SimArray([ 183.81436575, 59.55845491, -12.07430806, ..., -172.14401735, -173.37977902, -172.9047343 ], 'km s**-1'))
you can get a list of all available derivable quantities
s.derivable_keys()[0:10]
['tform', 'mass', 'temp', 'r_mag', 'i_mag', 'K_lum_den', 'v_mean', 'U_lum_den', 'vcxy', 'j2']
Defining your own derivable quantity requires a single decorator line and can be done inside any module that you import into your python session:
@pynbody.ramses.RamsesSnap.derived_quantity
def my_quantity(self) :
return self['x']*self['y']
s['my_quantity']
SimSnap: deriving array my_quantity
SimArray([ 1.46194737e+07, 1.01377668e+08, 1.29737299e+08, ..., 9.96316553e+09, 9.96339360e+09, 9.96324152e+09], 'kpc**2')
The first thing you probably do after running a cosmological simulation is to run it through a halo finder like hop to get halo information. A few helper functions are defined in pynbody.analysis.ramses_util
to load the halo data from a hop output and center the snapshot. Note that if hop data is not found, it will be automatically generated. Of course, if you are running this on your own machine, both hop and the script need to be in your path, or you must generate the data independently.
Hop output is in system units, so we need to convert back before centering:
s.original_units()
pynbody.analysis.ramses_util.hop_center(s)
SimSnap: converting my_quantity units from kpc**2 to 3.79e+53 cm**2 SimSnap: converting vt units from km s**-1 to 1.41e+09 cm s**-1 SimSnap: converting pos units from kpc to 6.16e+26 cm SimSnap: converting vr units from km s**-1 to 1.41e+09 cm s**-1 SimSnap: converting vel units from km s**-1 to 1.41e+09 cm s**-1 SimSnap: converting v2 units from km**2 s**-2 to 1.98e+18 cm**2 s**-2 SimSnap: converting r units from kpc to 6.16e+26 cm SimSnap: converting mass units from Msol to 5.92e+50 g SimSnap: converting pos units from 6.16e+26 cm to 6.16e+26 cm SimSnap: converting pos units from 6.16e+26 cm to 6.16e+26 cm SimSnap: converting pos units from 6.16e+26 cm to 6.16e+26 cm SimSnap: converting p units from km**2 Msol s**-2 kpc**-3 to 5.02e-12 g cm**-1 s**-2 SimSnap: converting smooth units from kpc to 6.16e+26 cm SimSnap: converting rho units from Msol kpc**-3 to 2.53e-30 g cm**-3
We've recentered the simulation on the most massive halo according to hop, but this usually isn't particularly accurate so we can redo the centering using a "shrinking sphere" method:
pynbody.analysis.halo.center(s[pynbody.filt.Sphere('100 kpc')],mode='ssc')
Finding halo velocity center... vcen= [ 0.00016491 -0.00013521 -0.00032374]
Lets look at the dark matter distribution of a part of the whole simulation box. Note that we can specify the desired units here as well -- by specifying a unit of surface density, we tell pynbody to render a projection:
s.physical_units();
pynbody.plot.image(s.d,width='100 Mpc', units='Msol kpc^-2', cmap=cm.Greys)
SimSnap: converting my_quantity units from 3.79e+53 cm**2 to kpc**2 SimSnap: converting pos units from 6.16e+26 cm to kpc SimSnap: converting vel units from 1.41e+09 cm s**-1 to km s**-1 SimSnap: converting mass units from 5.92e+50 g to Msol SimSnap: converting p units from 5.02e-12 g cm**-1 s**-2 to km**2 Msol s**-2 kpc**-3 SimSnap: converting smooth units from 6.16e+26 cm to kpc SimSnap: converting rho units from 2.53e-30 g cm**-3 to Msol kpc**-3 SimSnap: deriving array rho Building 16 trees with leafsize=16 Tree build done in 5.25 s Calculating SPH density SimSnap: deriving array smooth Smoothing with 32 nearest neighbours Smoothing done in 3.92s Density calculation done in 1.91 s
Rendering image on 16 threads...
SimArray([[ 3823456.75, 3902639.5 , 3981821.25, ..., 2355582. , 2325772. , 2295961.75], [ 3835513.5 , 3908890. , 3982266.25, ..., 2397100.25, 2367628.25, 2338156. ], [ 3847570.25, 3915140.25, 3982710.25, ..., 2438618.5 , 2409484.25, 2380350. ], ..., [ 4023282.5 , 4325098. , 4626914. , ..., 4168680.25, 4159233.75, 4149786.5 ], [ 4068084.5 , 4368396. , 4668708.5 , ..., 4486165. , 4479221. , 4472278.5 ], [ 4112886. , 4411695. , 4710503. , ..., 4803650.5 , 4799210. , 4794770. ]], dtype=float32, 'Msol kpc**-2')
Now that we are centered on the main halo, we can cut out a region of roughly the virial radius:
sph = pynbody.filt.Sphere('180 kpc')
h1 = s[sph]
h1
<SimSnap "/zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101:sphere" len=4007437>
This new object is a SubSnap
of the full simulation, and behaves in exactly the same way:
len(h1.d), len(h1.s), len(h1.g)
(1150790, 1171311, 1685336)
h1['mass'].sum().in_units('1e12 Msol')
SimArray(0.6467235093720013, '1.00e+12 Msol')
Typically we want to make sure the system is aligned in some predictable way, so pynbody allows us to make the disk angular momentum axis the z-axis:
pynbody.analysis.angmom.faceon(h1,cen=(0,0,0),disk_size='5 kpc')
Finding halo velocity center... vcen= [ 4.50101837e-13 -4.09682834e-13 -8.88203760e-13] Calculating angular momentum vector... Transforming simulation...
we use cen=(0,0,0)
because we have already centered it in the previous step... if you don't specify a center, faceon
centers the snapshot before aligning. The disk_size
keyword specifies the region which is used for calculating the angular momentum vector. Now lets look at a slice of the gas disk edge-on:
s.rotate_x(90)
pynbody.plot.image(s.g,width='50 kpc',denoise=True, cmap=cm.Blues);
Rendering image on 16 threads... Rendering image on 16 threads...
By default, rho
(i.e. density) is used to render the image, but any valid array is fair game:
pynbody.plot.image(s.g,qty='temp',width='50 kpc',denoise=True, cmap=cm.YlOrRd)
s.rotate_x(-90)
Rendering image on 16 threads... Rendering image on 16 threads...
SimSnap: deriving array temp
You can also make projections of weighted quantities:
pynbody.plot.image(s.g,qty='vr',width='40 kpc',denoise=True,av_z='rho',log=False, cmap=cm.RdBu);
SimSnap: deriving array vr SimSnap: deriving array
Rendering image on 16 threads... Rendering image on 16 threads... Rendering image on 16 threads... Rendering image on 16 threads...
r
pynbody
includes a Profile
class that implements some basic functionality for creating (you guessed it) profiles. For a more detailed description of its functionality head over to the Profile documentation.
To start with, lets create a 3D spherical profile of the total matter distribution using log-spaced bins, and additional profiles for the stars and gas:
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 0x3d31b50>
The Profile
class has lots of built-in functionality:
p.derivable_keys()
['j_theta', 'v_circ', 'E_circ', 'kappa', 'density', 'mass_enc', 'pot', 'magnitudes', 'Q', 'beta', 'jtot', 'mass', 'j_circ', 'X', 'dyntime', 'g_spherical', 'sb', 'fourier', 'omega', 'rotation_curve_spherical', 'j_phi']
s['eps']=s.g['smooth'].min()
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 10.3 s Profile: v_circ() -- warning, disk must be in the x-y plane Rotation curve calculated in 2.61 s Profile: v_circ() -- warning, disk must be in the x-y plane Rotation curve calculated in 3.81 s
<matplotlib.legend.Legend at 0x2df6910>
In addition, Profile
can use any loaded (or loadable/derivable!) array. If the data is not in-memory already, it will be lazy-loaded or calculated:
plot(ps['rbins'],ps['vr'])
xlabel('$R$ [kpc]')
ylabel('mean $V_R$ [km/s]')
SimSnap: deriving array vr SimSnap: deriving array r
<matplotlib.text.Text at 0x3db1590>
By default, Profile
calculates the means of quantities, but it trivially calculate derivatives, dispersions, 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 0x2ee1390>
plot(ps['rbins'],ps['vr_rms'])
xlabel('R [kpc]')
ylabel('$v_{rms}$')
Profile: auto-deriving vr_rms
<matplotlib.text.Text at 0x3df5290>
You can define filters to select particles based on their quantities. Under the hood, this just calls the numpy.where
function, but it makes your code much more readable. For example, to get stars that are 1-2 Gyr old:
metfilt = pynbody.filt.BandPass('metal', 1e-3,1e-2)
metfilt
BandPass('metal', 1.00e-03, 1.00e-02)
disc = pynbody.filt.Disc('20 kpc', '1 kpc')
h1[disc].s
<SimSnap "/zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101:sphere:disc::star" len=574050>
h1.s[disc]['mass'].sum().in_units('1e10 Msol')
SimArray(4.18613882677458, '1.00e+10 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[disc&metfilt]['mass'].sum().in_units('1e10 Msol')
SimArray(0.47000260266959104, '1.00e+10 Msol')
Often it's desirable to convert a Ramses output into some other format for post processing (running other halo finders, for example). This is quite easily achieved with pynbody, but the unit conversion can be tricky -- so we've created a few helper functions in pynbody.analysis.ramses_util
to make this less painful:
pynbody.analysis.ramses_util.convert_to_tipsy_fullbox(s.filename)
Loading using backend <class 'pynbody.ramses.RamsesSnap'> RamsesSnap: loading hydro files done SimSnap: deriving array 1 3 5 7 9 11 13 15 17 19 21 23 25 29 31 27 2 4 6 8 10 12 14 16 18 20 22 24 26 30 32 28 63 61 39 53 43 37 41 33 35 45 47 55 49 59 57 51 64 62 40 54 44 38 42 34 36 46 48 56 50 60 58 52 95 93 83 81 79 71 85 65 67 73 75 69 77 89 87 91 96 94 84 82 80 72 86 66 68 74 76 70 78 90 88 92 127 121 113 117 109 103 115 97 99 105 107 101 111 123 119 125 128 122 114 118 110 104 116 98 100 106 108 102 112 124 120 126 temp TipsySnap: writing main file as /zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101_fullbox.tipsy TipsySnap: writing auxiliary arrays
This also creates a .param file with requisite units information in order to inform pynbody about the units system.
If you don't want the full box for some reason (e.g. to speed up analysis) and only want to save a part of the simulation for quick-and-dirty post processing, you can do this with another function:
pynbody.analysis.ramses_util.convert_to_tipsy_simple(s.filename,pynbody.filt.Sphere('200 kpc'))
Loading using backend <class 'pynbody.ramses.RamsesSnap'> RamsesSnap: loading hydro files
Finding halo velocity center... vcen= [ 0.00016491 -0.00013521 -0.00032374] Calculating angular momentum vector... Transforming simulation...
done SimSnap: deriving array 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 63 59 39 51 43 37 41 33 35 45 47 55 49 53 61 57 64 60 40 52 44 38 42 34 36 46 48 56 50 54 62 58 95 87 75 89 81 69 79 65 67 73 77 71 83 85 91 93 96 88 76 90 82 70 80 66 68 74 78 72 84 86 92 94 127 117 113 121 119 101 105 97 99 109 111 103 107 115 125 123 128 118 114 122 120 102 106 98 100 110 112 104 108 116 126 124
[ 2213.59814674 2996.23889034 1685.19782238 ..., 268515.38612534 261518.35324757 261736.69514807]
temp SimSnap: deriving array temp TipsySnap: writing main file as output_00101.tipsy RamsesSnap: loading grav files done TipsySnap: writing auxiliary arrays 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 63 59 43 49 39 37 41 33 35 45 47 57 53 51 61 55 64 60 44 50 40 38 42 34 36 46 48 58 54 52 62 56 95 87 73 81 77 69 85 65 67 79 75 89 83 91 93 71 96 88 74 82 78 70 86 66 68 80 76 90 84 92 94 72 127 117 111 115 113 101 119 97 99 105 109 123 107 125 121 103 128 118 112 116 114 102 120 98 100 106 110 124 108 126 122 104
The newly-generated output in tipsy format can now be used in the same way as the RAMSES output:
stipsy = pynbody.load('/zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101_fullbox.tipsy')
Loading using backend <class 'pynbody.tipsy.TipsySnap'> TipsySnap: loading /zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101_fullbox.tipsy
stipsy
<SimSnap "/zbox/data/roskar/cosmo/AMR/rad_fbk/std/output_00101_fullbox.tipsy" len=23196808>
len(stipsy.s),len(stipsy.g)
(1243472, 14817745)
stipsy.properties
{'a': 1.00001844594404, 'boxsize': Unit("2.00e+05 kpc a"), 'h': 0.703969087746729, 'omegaL0': 0.728, 'omegaM0': 0.272, 'time': Unit("1.41e+01 s kpc km**-1")}
We have previously generated an Amiga Halo Finder halo catalogue for this snapshot -- so now we can load the halos and use the halo catalogue framework:
h = stipsy.halos()
AHFCatalogue: loading particles... halos... substructure... done!
pynbody.analysis.halo.center(h[1],mode='ssc')
Finding halo velocity center... vcen= [-0.02096087 -0.0040748 -0.01705168]
TipsySnap: loading data from main file
pynbody.plot.image(stipsy.d,width='100 Mpc',units='Msol kpc^-2',cmap=cm.Greys);
SimSnap: deriving array rho Building 16 trees with leafsize=16 Tree build done in 5.05 s Calculating SPH density SimSnap: deriving array smooth Smoothing with 32 nearest neighbours Smoothing done in 3.98s Density calculation done in 1.89 s
Rendering image on 16 threads...