Profiles

The Profile class is a general-purpose class that can be used to bin information in a simulation in various ways. The Profile class is designed to be an extension of the syntax implemented in the SimSnap class.

Creating a Profile instance defines the bins, after which specific calculations can be carried out.

To illustrate, let’s load a snapshot:

In [1]: import pynbody;
   ...: from pynbody.analysis import profile;
   ...: import matplotlib.pylab as plt
   ...: 

In [2]: s = pynbody.load('testdata/gasoline_ahf/g15784.lr.01024')

In [3]: s.physical_units()

In [4]: h = s.halos()

In [5]: pynbody.analysis.faceon(h[0])
Out[5]: <Transformation translate, offset_velocity, rotate>

The final command here centres the origin on the main halo and puts the disk in the xy-plane (see Aligning the Snapshot). Now we can make the Profile instance:

In [6]: p = profile.Profile( h[0].star, rmin = '.05 kpc', rmax = '50 kpc')

Note

You can pass either floating point values or unit strings to rmin and rmax.

With the default parameters, a 2D profile is created, in the xy-plane. Our use of the faceon() function ensures that the stellar disk is ready for analysis in this way. We can now plot the density profile:

In [7]: plt.plot(p['rbins'].in_units('kpc'),p['density'].in_units('Msol kpc^-2'),'k')
Out[7]: [<matplotlib.lines.Line2D at 0x7f19684bfb50>]

In [8]: plt.xlabel('$R$ [kpc]'); plt.ylabel(r'$\Sigma_{\star}$ [M$_{\odot}$ kpc$^{-2}$]'); plt.semilogy()
Out[8]: []

The binning is linear by default, but we can see here that the profile is not well-sampled at large radii. We can change the binning to logarithmic by specifying type='log':

In [9]: p = profile.Profile( h[0].star, rmin = '.05 kpc', rmax = '50 kpc', type='log')

In [10]: plt.plot(p['rbins'].in_units('kpc'),p['density'].in_units('Msol kpc^-2'),'k')
Out[10]: [<matplotlib.lines.Line2D at 0x7f19681791d0>]

In [11]: plt.xlabel('$R$ [kpc]'); \
   ....: plt.ylabel(r'$\Sigma_{\star}$ [M$_{\odot}$ kpc$^{-2}$]'); \
   ....: plt.semilogy()
   ....: 
Out[11]: []

To make a spherically-symmetric 3D profile, specify ndim=3 when creating the profile.

In [12]: pdm_3d = profile.Profile(s.dm, rmin = '.01 kpc', rmax = '500 kpc', ndim = 3)

Even though we use s.dm here (i.e. dark matter from the full snapshot, not just halo 0), the whole snapshot is still centered on halo 0 following our earlier call to faceon(). This allows us to explore that far outer reaches of the halo around the galaxy. Let’s now plot the dark matter density profile:

In [13]: plt.plot(pdm_3d['rbins'].in_units('kpc'),pdm_3d['density'].in_units('Msol kpc^-3'),'k')
Out[13]: [<matplotlib.lines.Line2D at 0x7f19684e7310>]

In [14]: plt.xlabel('$r$ [kpc]'); plt.ylabel(r'$\rho_{\rm DM}$ [M$_{\odot}$ kpc$^{-3}$]'); plt.loglog()
Out[14]: []

Mass-weighted average quantities

The above examples illustrate the most basic use of profiling, to generate binned density estimates. One may also generate mass-weighted averages of any quantity that is either stored in the snapshot or derivable from it. For example, the sample snapshot being used above has metallicity information from which an Fe/H estimate can be derived by pynbody.

In [15]: plt.plot(p['rbins'].in_units('kpc'),p['feh'],'k')
Out[15]: [<matplotlib.lines.Line2D at 0x7f196860ebd0>]

In [16]: plt.xlabel('$R$ [kpc]'); plt.ylabel('[Fe/H]')
Out[16]: Text(0, 0.5, '[Fe/H]')

Special quantities

As well as straight-forward densities and mass-weighted averages, there are a number of special profiling functions implemented. To see a full list, use the pynbody.analysis.profile.Profile.derivable_keys() method or consult the list of functions in pynbody.analysis.profile.

For example, the mass enclosed within a given radius is given by mass_enc:

In [17]: plt.plot(p['rbins'].in_units('kpc'), p['mass_enc'], 'k')
Out[17]: [<matplotlib.lines.Line2D at 0x7f196b541bd0>]

In [18]: plt.xlabel('$R$ [kpc]'); plt.ylabel(r'$M_{\star}(<R)$')
Out[18]: Text(0, 0.5, '$M_{\\star}(<R)$')

See the Profile documentation for a full list with brief descriptions. You can also check the available profiles in your session using derivable_keys().

Note

You can also define your own profiling functions in your code by using the Profile.profile_property decorator; these become available in just the same way as the built-in profiling functions. If you wish to do this, the best place to start is by studying the implementation of the existing profile properties in the profile module.

Surface brightnesses

Some of the derivable quantities take parameters. For example, surface brightness profiles are given by sb and on consulting the docstring, this turns out to take the band as an input. Parameters are passed in to the string using commas. For example, to get the Johnson U-band surface brightness profile, we ask for sb,u, or for R-band sb,r:

In [19]: plt.plot(p['rbins'].in_units('kpc'), p['sb,u'], 'b', label="U band");
   ....: plt.plot(p['rbins'].in_units('kpc'), p['sb,r'], 'r', label="R band");
   ....: 

In [20]: plt.xlabel('$R$ [kpc]'); plt.ylabel(r'SB/mag/arcsec$^2$');
   ....: plt.legend()
   ....: 

Note

Surface brightnesses are calculated using SSP tables described further in the luminosity module.

Rotation curves

Another useful special quantity is the rotation curve, which can be calculated using the v_circ key:

In [21]: p_dm = pynbody.analysis.profile.Profile(h[0].dm, min=.05, max=50, type = 'log')

In [22]: p_gas = pynbody.analysis.profile.Profile(h[0].gas, min=.05, max=50, type = 'log')

In [23]: p_all = pynbody.analysis.profile.Profile(h[0], min=.05, max=50, type = 'log')

In [24]: for prof, name in zip([p_all, p_dm, p, p_gas],['total', 'dm', 'stars', 'gas']):
   ....:     plt.plot(prof['rbins'], prof['v_circ'], label=name)
   ....: 

In [25]: plt.xlabel('$R$ [kpc]');

In [26]: plt.ylabel('$v_{circ}$ [km/s]');

In [27]: plt.legend()
Out[27]: <matplotlib.legend.Legend at 0x7f1968190e50>

As the above example makes clear, the circular velocity is estimated from the gravitational force generated by particles known to the profile object, rather than the entire snapshot.

Calculating Derivatives and Dispersions

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 [28]: p_all = profile.Profile(s, rmin='.01 kpc', rmax='250 kpc')

In [29]: p_all['pot'][0:10] # returns the potential profile
Out[29]: 
SimArray([-1883673.73832153, -1774864.42980055, -1722130.45724663,
          -1690546.57353947, -1668762.42991661, -1652787.25382017,
          -1640343.20373383, -1630153.24386831, -1621515.5868703 ,
          -1614131.6419176 ], 'km**2 s**-2')

In [30]: p_all['d_pot'][0:10] # returns d phi / dr from p["phi"]
Out[30]: 
array([Unit("4.35e+04 kpc**-1"), Unit("3.23e+04 kpc**-1"),
       Unit("1.69e+04 kpc**-1"), Unit("1.07e+04 kpc**-1"),
       Unit("7.55e+03 kpc**-1"), Unit("5.68e+03 kpc**-1"),
       Unit("4.53e+03 kpc**-1"), Unit("3.77e+03 kpc**-1"),
       Unit("3.20e+03 kpc**-1"), Unit("2.77e+03 kpc**-1")], dtype=object)

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 [31]: plt.plot(p['rbins'].in_units('kpc'), p['vr_disp'].in_units('km s^-1'), 'k')
Out[31]: [<matplotlib.lines.Line2D at 0x7f196b6798d0>]

In [32]: plt.xlabel('$R$ [kpc]'); \
   ....: plt.ylabel('$\sigma_{r}$')
   ....: 
Out[32]: Text(0, 0.5, '$\\sigma_{r}$')

In addition to doing this by hand, you can make a QuantileProfile that can return any desired quantile range. By default, this is the mean +/- 1-sigma:

In [33]: p_quant = profile.QuantileProfile( h[0].s, rmin = '0.1 kpc', rmax = '50 kpc')

In [34]: plt.clf(); plt.plot(p_quant['rbins'], p_quant['feh'][:,1], 'k')
Out[34]: [<matplotlib.lines.Line2D at 0x7f19681bd810>]

In [35]: plt.fill_between(p_quant['rbins'], p_quant['feh'][:,0], p_quant['feh'][:,2], color = 'Grey', alpha=0.5)
Out[35]: <matplotlib.collections.PolyCollection at 0x7f196b6555d0>

In [36]: plt.xlabel('$R$ [kpc]'); plt.ylabel('[Fe/H]')
Out[36]: Text(0, 0.5, '[Fe/H]')

Vertical Profiles

For analyzing disk structure, it is frequently useful to have a profile in the z-direction. This is done with the VerticalProfile which behaves in the same way as the Profile. Unlike in the basic class, you must specify the radial range and maximum z to be used:

In [37]: p_vert = profile.VerticalProfile( h[0].s, '3 kpc', '5 kpc', '5 kpc')

In [38]: plt.clf(); plt.plot(p_vert['rbins'].in_units('pc'), p_vert['density'].in_units('Msol pc^-3'),'k')
Out[38]: [<matplotlib.lines.Line2D at 0x7f196b607b90>]

In [39]: plt.xlabel('$z$ [pc]'); plt.ylabel(r'$\rho_{\star}$ [M$_{\odot}$ pc$^{-3}$]')
Out[39]: Text(0, 0.5, '$\\rho_{\\star}$ [M$_{\\odot}$ pc$^{-3}$]')

Profiles with arbitrary x-axes

Radial profiles are nice, but sometimes we want a profile using a different quantity on the x-axis. We might want to know, for example, how the mean metallicity varies as a function of age in the stars. Profile by default uses either the 3D or xy-plane radial distance, depending on the value of ndim. But 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 0 according to their age:

In [40]: s.s['age'].convert_units('Gyr')

In [41]: p_age = profile.Profile( h[0].s,
   ....:                          calc_x = lambda x: x.s['age'],
   ....:                          rmax = '10 Gyr' )
   ....: 

In [42]: plt.clf(); plt.plot(p_age['rbins'], p_age['feh'], 'k', label = 'mean [Fe/H]')
Out[42]: [<matplotlib.lines.Line2D at 0x7f196b78d190>]

In [43]: plt.plot(p_age['rbins'], p_age['feh_disp'], 'k--', label = 'dispersion')
Out[43]: [<matplotlib.lines.Line2D at 0x7f196b665310>]

In [44]: plt.xlabel('Age [Gyr]'); plt.ylabel('[Fe/H]')
Out[44]: Text(0, 0.5, '[Fe/H]')

In [45]: plt.legend()
Out[45]: <matplotlib.legend.Legend at 0x7f196b7e66d0>