Have you ever wondered how to get the streamlines of the velocity field of a N-body simulation? In practice, streamlines are an extraordinary tool to get insight on what is going on in a simulation by inspecting a single snapshot. They are also extremely useful in many applications that are outside the scope of this post.

Let me just show an example of their power. Figure 1 shows the streamlines of the velocity field of a cosmological N-body simulation. Although the underlying density field is not shown, we can infer its presence by looking at the flux lines. Indeed, we can infer the presence of a “void” at (x,y) = (10,-20) because the flux lines point outward. We can also infer the location of filaments by looking at the intersection of opposing lines.

Doing this kind of diagrams is relatively straightforward once we have the velocity field calculated on a regular grid. However, for N-body simulations (or hydrodynamical SPH simulations) this is not always the case. Here is where Py-SPHViewer can give us a hand.

With Py-SPHViewer we are able to calculate the instantaneous (and continuous) mean velocity field, which can then be passed to matplotlib to construct the streamlines.

In this post, I’ll give the recipe to get incredible nice results in only a few steps.

### Getting the data for this example and drawing the streamlines

In this example I use an output of a low-resolution N-Body cosmological simulation ran with Gadget-2.

The initial conditions were created using Music, and I used my own laptop to run it. You may want to download this file to continue with the example. Figure 2 shows a scatter plot of the data.

In order to get the velocity field of the box we can consider different slices, or just project the full box. I’ll use the latter, and consider the mean (projected) velocity field. In particular, I will consider the density-weighted average. Because of this, the first thing to do is to calculate the density field of the box:

This can be done with QuickView:

import h5py from sphviewer.tools import QuickView #Read the file snap = h5py.File('darkmatter_box.h5py', 'r') pos = snap['PartType1/Coordinates'].value vel = snap['PartType1/Velocities'].value hsml = snap['PartType1/SmoothingLength'].value #Calculate the density field qv = QuickView(pos.T, hsml=hsml, r='infinity', x=32, y=32, z=32,plot=False, extent=[-32,32,-32,32], logscale=False) density_field = qv.get_image() extent = qv.get_extent()

Fig. 3 shows the resulting density field. We now need to do the same for the x and y components of the velocity field:

vfield = [] for i in xrange(2): qv = QuickView(pos.T, vel[:,i], hsml=hsml, r='infinity', x=32, y=32, z=32,plot=False, extent=[-32,32,-32,32], logscale=False) vfield.append(qv.get_image()/density_field)

Where we have smoothed the velocity field (x and y components separately), and divided the resulting field by the density field, which effectively produces a density-weighted average. The resulting fields are shown in Fig. 4

We have all the ingredients to let matplotlib do the magic. Using the streamplot method:

fig = plt.figure(1, figsize=(7,7)) ax = fig.add_subplot(111) X = np.linspace(extent[0], extent[1], 500) Y = np.linspace(extent[2], extent[3], 500) ax.imshow(np.log10(density_field), origin='lower', extent=extent, cmap='bone') v = np.log10(np.sqrt(vfield[0]**2+vfield[1]**2)) color = v/np.max(v) lw = 2*color streams = ax.streamplot(X,Y,vfield[0], vfield[1], color=color, density=1.5, cmap='jet', linewidth=lw, arrowsize=1) ax.set_xlim(extent[0],extent[1]) ax.set_ylim(extent[2],extent[3]) ax.minorticks_on() ax.set_xlabel(r'$\rm X / \ Mpc \ h^{-1}$', size=25) ax.set_ylabel(r'$\rm Y / \ Mpc \ h^{-1}$', size=25) fig.savefig('figura5.png', bbox_inches='tight')

we obtain the output shown in Fig. 5.

## Conclusions

Streamlines are straightforward to calculate for N-body simulations, provided the velocity field is known on a regular grid. The velocity field can be calculated trivially using Py-SPHViewer, and streamlines can be constructed with matplotlib.

Finally, let me just add that streamlines can also be animated to give a really deep insight on what is going on in the simulation. I plan to focus on this in a forthcoming post. Please, do not hesitate to leave your comments, suggestions or questions down below.