Working with Simulations

Sometimes, you need to work with many snapshots across an entire simulation instead of a single snapshot at once. This notebook describes the routines included in snapAnalysis.sim for extracting useful data from entire simualtions. We’ll work with the isolated dark matter halo snapshots included in the test suite. There’s only two snapshots in this test simulation (for space reasons), but these recipes translate directly to simulations including thousands of snapshots.

from snapanalysis import sim
from snapanalysis.snap import snapshot
import astropy.units as u

Obtaining a list of available snapshots

Sometimes, it’s convenient to get a list of all the snapshot files in a particular simulation directory.

sim_dir = "../tests/example_snaps/cdm_snaps"
snaps = sim.get_snaps(sim_dir)
print(snaps)

['../tests/example_snaps/cdm_snaps/snap_000.hdf5'
 '../tests/example_snaps/cdm_snaps/snap_001.hdf5']

Extracting the center-of-mass orbit of a galaxy

To extract the center-of-mass orbit of a galaxy, use the orbit_com function. In a simulation that includes multiple galaxies, you must specify the particle ID range that belongs to the galaxy of interest. You can also supply dictionaries of kwargs to the position and velocity center-finding routines. The columns in the output array are time, x, y, z, vx, vy, vz

sim.orbit_com(sim_dir, part_type=1, select_IDs=None, com_kwargs={"vol_dec": 1.1, "r_start":15.0*u.kpc, "guess": [0,0,0]*u.kpc})

array([[ 0.      ,  0.120152, -0.287273,  0.114837, -2.411327,  0.092466,
        -1.306028],
       [ 0.197182, -0.426582, -0.400833, -0.421794, -0.098579, -0.782643,
        -0.2411  ]])

Extracting rotations

To extract the rotation matrices needed at each time to align a galaxy’s angular momentum to the \(z\)-axis, use the get_alignment_rotations function. Again, you can specify an ID range to select a particular galaxy. By default, get_alignment_rotations will center each snapshot before calculating its angular momentum, but you can also specify a file of precomputed centers from orbit_com if you have one. The resulting array will be the stacked rotation matrices required to align each snapshot’s angular momentum with \(z\), of the shape N_snapshots x 3 x 3.

sim.get_alignment_rotations(sim_dir, part_type=1, r_max=10*u.kpc, select_IDs=None, use_centers=None)

/Users/haydenfoote/Documents/gradschool/research/code/snapAnalysis/src/snapanalysis/snap.py:468: UserWarning: COM: Minimum number of particles reached before COM converged.
  warnings.warn(

array([[[-0.95041252,  0.05663625, -0.30579141],
        [-0.0594857 , -0.99822916,  0.        ],
        [-0.3052499 ,  0.01819022,  0.95209853]],

       [[ 0.53087062,  0.66868122, -0.52061676],
        [-0.78319123,  0.62178091,  0.        ],
        [ 0.32370956,  0.40774248,  0.85379048]]])

Extracting particle orbits

To extract the orbit of a specific particle(s), use get_particle_orbit. Again, by default this function will center and align the snapshots for you, but you may specify precomputed centers and rotations.

import numpy as np

# get 10 random particle ID from the first snapshot
s = snapshot(snaps[0], 1)
s.load_particle_data(["ParticleIDs"])
rand_ids = np.random.choice(s.data_fields["ParticleIDs"], size=10, replace=False)

# and extract their orbits
orbits = sim.get_particle_orbit(sim_dir, part_type=1, ids=rand_ids)

This function returns an orbit object, which keeps track of particle orbits. orbit objects are designed for ease of use with the the Gala package for galactric dynamics. You can access the data as follows.

# The times of each stored point in the orbit are stored in
print(orbits.t)

[0.         0.19718176] Gyr

# positions and velocities of each particle are stored in orbits.pos and orbits.vel.
# these have shape (N_times, 3, N_particles) for consistency with Gala.
print(orbits.pos.shape)
print(orbits.vel.shape)

(2, 3, 10)
(2, 3, 10)

# The IDs attribute is a dictionary whose keys are the particle IDs,
# and whose values are the indices of those particles in the data arrays.
print(orbits.ids)

{'1881': 0, '3490': 1, '4334': 2, '4556': 3, '4881': 4, '6363': 5, '6548': 6, '7777': 7, '8855': 8, '9644': 9}

# Therefore, to access the position of a particular particle at all times,
id = str(int(rand_ids[0]))
print(orbits.pos[:, :, orbits.ids[id]])

[[-14.05476249 -11.31392662 -15.47666254]
 [ 11.23177433   2.7399446  -12.70862759]] kpc