import yt
	import glob
	from yt.analysis_modules.particle_trajectories.api import ParticleTrajectories

	# If your particles are separate from your grid data, you have to do this:

	f_plts = glob.glob("/work/02418/rkashyap/flash4/super3D/run_89_small/super3d_hdf5_plt_cnt_0[0-2][0-9][0-9]")
	f_plts.sort()

	f_parts = glob.glob("/work/02418/rkashyap/flash4/super3D/run_89_small/super3d_hdf5_part_0[0-2][0-9][0-9]")
	f_parts.sort()

	fields = ["particle_dens","particle_temp","particle_pres","particle_ni56","particle_velocity_x", "particle_velocity_y", "particle_velocity_z"]

	#fields = ["particle_dens","particle_velocity_x", "particle_velocity_y", "particle_velocity_z"]

	indices = [1,2,3,67537,35164,64875]

	#indices = [64875]

	pfs = []

	for f_plt, f_part in zip(f_plts, f_parts):
	    pfs.append(yt.load(f_plt, particle_filename=f_part))

	ts = yt.DatasetSeries(pfs)

	trajs = ParticleTrajectories(ts, indices, fields=fields)


	# If your particles are in the same HDF5 file as your grid data,
	# you can call ParticleTrajectories like this:
	# trajs = ParticleTrajectories(filenames, indices, fields=fields)
	# without having to load each one up as above--it will do it for you.

	trajs.write_out_h5("super3D_89_trajectories.h5") # This will write all trajectories to a single file