Making profiles of all kinds of quantities is easy -- here's a simple example that shows how to plot a density profile:
tutorials/example_code/density_profile.py
Below is a more extended description of the ~pynbody.analysis.profile module.
The ~pynbody.analysis.profile.Profile class is meant to be a general-purpose class to satisfy (almost) all simulation profiling needs. Profiles are the most elementary ways to begin analyzing a simulation, so the ~pynbody.analysis.profile.Profile class is designed to be an extension of the syntax implemented in the ~pynbody.snapshot.SimSnap class and its derivatives.
Creating a ~pynbody.analysis.profile.Profile instance simply defines the bins etc. Importantly, it also stores lists of particle indices corresponding to each bin, so you can easily identify where the particles belong.
Let's quickly load a snapshot:
In [1]: import pynbody; from pynbody.analysis import profile; import matplotlib.pylab as plt
In [1]: s = pynbody.load('testdata/g15784.lr.01024'); s.physical_units()
In [2]: h = s.halos()
In [3]: pynbody.analysis.angmom.faceon(h[1])
We're centered on the main halo (as in the cookbook example above) and we make the ~pynbody.analysis.profile.Profile instance:
In [3]: p = profile.Profile(h[1].s, rmin='.01 kpc', rmax='50 kpc')
With the default parameters, the profile is made in the xy-plane. To make a spherically-symmetric 3D profile, specify ndim=3 when creating the profile.
In [3]: pdm_sph = profile.Profile(s.d, rmin = '.01 kpc', rmax = '250 kpc')
Even though we use s.d here (i.e. the full snapshot, not just halo 1), the whole snapshot is still centered on halo 1.
Note
You can pass unit strings to rmin and rmax and the conversion will be done automatically into whatever the current units of the snapshot are so you don't have to explicitly do any unit conversions.
Many profiling functions are already implemented -- see the ~pynbody.analysis.profile.Profile documentation for a full list with brief descriptions. You can also check the available profiles in your session using ~pynbody.analysis.profile.Profile.derivable_keys just like you would for a ~pynbody.snapshot.SimSnap:
In [3]: p.derivable_keys()
Additionally, any existing array can be 'profiled'. For example, if the metallicity [Fe/H] is a derived field stored under 'feh', then plotting a metallicity profile is as simple as:
In [4]: plt.plot(p['rbins'].in_units('kpc'),p['feh'],'k')
@savefig profile_fig1.png width=5in In [5]: plt.xlabel('$R$ [kpc]'); plt.ylabel('[Fe/H]')
If the array doesn't exist but is deriveable (check with s.derivable_keys()), it is automatically calculated.
In addition, you can define your own profiling functions in your code by using the Profile.profile_property decorator:
@profile.Profile.profile_property
def random(self):
import numpy as np
return np.random.rand(self.nbins)Now this will be automatically derivable for any newly-created profile as 'random'.
You can calculate derivatives of profiles automatically. For instance, you might be interested in d phi / dr if you're looking at a disk. This is as easy as attaching a d_ to the profile name. For example:
In [6]: p_all = profile.Profile(s, rmin='.01 kpc', rmax='250 kpc')
In [6]: p_all['pot'][0:10] # returns the potential profile
In [7]: p_all['d_pot'][0:10] # returns d phi / dr from p["phi"]
Similarly straightforward is the calculation of dispersions and root-mean-square values. You simply need to attach a _disp or _rms as a suffix to the profile name. To get the stellar velocity dispersion:
In [7]: plt.clf(); plt.plot(p['rbins'].in_units('kpc'),p['vr_disp'].in_units('km s^-1'),'k')
@savefig profile_fig2.png width=5in In [6]: plt.xlabel('$R$ [kpc]'); plt.ylabel('$sigma{r}$')
In addition to doing this by hand, you can make a ~pynbody.analysis.profile.QuantileProfile that can return any desired quantile range. By default, this is the mean +/- 1-sigma:
In [5]: p_quant = profile.QuantileProfile(h[1].s, rmin = '0.1 kpc', rmax = '50 kpc')
In [6]: plt.clf(); plt.plot(p_quant['rbins'], p_quant['feh'][:,1], 'k')
In [6]: plt.fill_between(p_quant['rbins'], p_quant['feh'][:,0], p_quant['feh'][:,2], color = 'Grey', alpha=0.5)
@savefig profile_quant.png width=5in In [6]: plt.xlabel('$R$ [kpc]'); plt.ylabel('[Fe/H]')
Radial profiles are nice, but sometimes we want a "profile" using a different quantity. We might want to know, for example, how the mean metallicity varies as a function of age in the stars. ~pynbody.analysis.profile.Profile calls the function ~pynbody.analysis.profile.Profile._calculate_x by default and this simply returns the 3D or xy-plane radial distance, depending on the value of ndim. We can specify a different function using the calc_x keyword. Often these are simple so a lambda function can be used (e.g. if we just want to return an array) or can also be more complicated functions. For example, to make the profile of stars in halo 1 according to their age:
In [6]: s.s['age'].convert_units('Gyr')
In [5]: p_age = profile.Profile(h[1].s, calc_x = lambda x: x.s['age'], rmax = '10 Gyr')
In [6]: plt.clf(); plt.plot(p_age['rbins'], p_age['feh'], 'k', label = 'mean [Fe/H]')
In [6]: plt.plot(p_age['rbins'], p_age['feh_disp'], 'k--', label = 'dispersion')
In [6]: plt.xlabel('Age [Gyr]'); plt.ylabel('[Fe/H]')
@savefig profile_fig4.png width=5in In [6]: plt.legend()
For analyzing disk structure, it is frequently useful to have a profile in the z-direction. This is done with the ~pynbody.analysis.profile.VerticalProfile which behaves in the same way as the ~pynbody.analysis.profile.Profile. Unlike in the basic class, you must specify the radial range and maximum z to be used:
In [5]: p_vert = profile.VerticalProfile(h[1].s, '3 kpc', '5 kpc', '5 kpc')
In [5]: plt.clf(); plt.plot(p_vert['rbins'].in_units('pc'), p_vert['density'].in_units('Msol pc^-3'),'k')
@savefig profile_fig5.png width=5in In [5]: plt.xlabel('$z$ [pc]'); plt.ylabel(r'$rho{star}$ [M$_{odot}$ pc$^{-3}$]')
Similarly, one can make inclined profiles using the ~pynbody.analysis.profile.VerticalProfile, but the snapshot needs to be rotated first:
In [5]: s.rotate_x(60) # rotate the snapshot by 60-degrees
In [5]: p_inc = profile.InclinedProfile(h[1].s, 60, rmin = '0.1 kpc', rmax = '50 kpc')