Example: Getting the Most Bound Particle

In this example, I show how to get the most bound particle from a halo. This example is also in the example/mbp/ path. The plotting_interactive.py will let matplotlib open an interactive window for the 3D rendering of the mock halo. The plotting.py will save a rendered .png file to the exmaple/mbp/ path so you can see a rendered 3d visualization. Here, I will walk through what each step does.

Note: Before you run this file, be sure that you have generated the test data. To do this, type the following into your terminal:

cd path/to/GROUPy/tests/test_data/
python generate_data.py

when you ls what's in your directory, you should see halos_snapshot_*.hdf5.*.particles files, out_*.list files, and and file called output_snapshot_times.txt. This is the test data that we will be using to run this test.

Getting the Data Ready for Harvest

First, get the particle data files ready to input. Notice how I put the data path minus the *.particle suffix

pd_base = "../../tests/test_data/halos_snapshot_%i.hdf5"
particle_data = [pd_base%i for i in range(10)]

Next, get the halo data files ready to be input. For the actual ROCKSTAR data, you should only use the out_*.list file, as this file has the data on what the ID of the halo is in the next timestep. To reiterate: do not use the .ascii file.

rs_base = "../../tests/test_data/out_%i.list"
halo_data = [rs_base%i for i in range(10)]

Next, get the time information ready to input into the halo. These times need to match up with the a values in each snapshot.

time_file = "../../tests/test_data/output_snapshot_times.txt"
tf = open(time_file, 'r')
times = [float(l) for l in tf]

Creating the Halo object and getting the most bound particle

Let's choose to find the most bound particle for snapshot 4. I chose the file_start to also be 4, but this is not necessary.

t = 4
h = Halo(0, halo_data, times, particle_data, file_start=t)

Now calculate the most bound particle in the halo at time 4. The get_mb_particle function stores the most bound particle, but it also returns it.

MBP = h.get_mb_particle(t)

For fun, let's also just calculate the most bound particle's kinetic energy and potential energy:

K = KE(MBP, t, (h.vx[t],h.vy[t],h.vz[t]))
P = PEg(MBP, h.particle_list, t)

print("Kinetic Energy: %.3e \nPotential Energy: %.3e"%(K,P))

Depending on the data that was generated, your output will look slightly different. However, mine ended up looking like the following:

> Kinetic Energy: 8.652e-33 
> Potential Energy: -1.052e-32