Added code snippet to show how to extract particles from a given halo

doc/RTD/source/VELOCIraptorInterface/output.rst

@@ -58,6 +58,64 @@ Besides the ``.catalog_particles`` file, there is also a
 but only for the unbound particles, a particle can only be present in one of
 these two lists. 
 
+Extracting the particles in a given halo
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The ``.catalog_particles`` file returns particle IDs that need to be matched
+with those in your snapshot to find the particles in the file that you
+wish to extract. The python snippet below should give you an idea of how to
+go about doing this for the bound particles.
+
+First, we need to extract the offset from the ``.catalog_group`` file, and
+work out how many _bound_ partices are in our halo. We can do this by
+looking at the next offset. Then, we can ID match those with the snapshot
+file, and get the mask for the _positions_ in the file that correspond
+to our bound particles. (Note this requires ``numpy > 1.15.0``).
+
+.. code-block:: python
+   :linenos:
+
+   import numpy as np
+   import h5py
+
+   snapshot_file = h5py.File("swift_snapshot.hdf5", "r")
+   group_file = h5py.File("velociraptor_output.catalog_group", "r")
+   particles_file = h5py.File("velociraptor_output.catalog_particles", "r")
+
+   halo = 100
+   # Grab the start position in the particles file to read from
+   halo_start_position = group_file["Offset"][halo]
+   halo_end_position = group_file["Offset"][halo + 1]
+   # We're done with that file now, best to close earlier rather than later
+   group_file.close()
+
+   # Get the relevant particle IDs for that halo; this includes particles
+   # of _all_ types.
+   particle_ids_in_halo = particles_file["Particle_IDs"][
+       halo_start_position:halo_end_position
+   ]
+   # Again, we're done with that file.
+   particles_file.close()
+
+   # Now, the tricky bit. We need to create the correspondance between the
+   # positions in the snapshot file, and the ids.
+
+   # Let's look for the dark matter particles in that halo.
+   particle_ids_from_snapshot = snapshot_file["PartType1/ParticleIDs"][...]
+
+   _, indices_v, indices_p = np.intersect1d(
+       particle_ids_in_halo,
+       particle_ids_from_snapshot,
+       assume_unique=True,
+       return_indices=True,
+   )
+
+   # indices_p gives the positions in the particle file where we will find
+   # the co-ordinates that we're looking for! To get the positions of all of
+   # those particles,
+   particle_positions_in_halo = snapshot_file["PartType1/Coordinates"][indices_p]
+
+
 Catalog_parttypes file
 ----------------------