"""
Using the particles in a late-redshift halo, figure out the enclosing box
for those particles over a series of older redshifts.
"""
from yt.mods import *
import gc
pf_now = load('RD0005/RedshiftOutput0005', data_style="enzo_packed_3d")
hID = 211
# First find the haloes from this run
haloes = parallelHF(pf_now, threshold=80.0, dm_only=True)
# Get the particle IDs that end up in and very close to the halo.
halo_IDs = haloes[hID]['particle_index']
# Now get the data from back in time over all the older dumps and
# write out the extrema to a file.
fp = open("refine_z.txt", "w")
fudge = 0.03
xc = haloes[hID].center_of_mass()[0]
yc = haloes[hID].center_of_mass()[1]
zc = haloes[hID].center_of_mass()[2]
for i in xrange(126):
print i
pf_old = load('DD%04d/data%04d' % (i,i), data_style="enzo_packed_3d")
red = pf_old['CosmologyCurrentRedshift']
all_old = pf_old.h.all_data()
old_IDs = all_old["particle_index"]
# Pick the particles in my halo
select = na.in1d(old_IDs, halo_IDs)
# Now we get the positions of these particles
old_xp = all_old["particle_position_x"][select]
old_yp = all_old["particle_position_y"][select]
old_zp = all_old["particle_position_z"][select]
# Find the extrema in code units, and re-center them by our offsets.
max_x = na.max(old_xp) + 0.5 - xc
max_y = na.max(old_yp) + 0.5 - yc
max_z = na.max(old_zp) + 0.5 - zc
min_x = na.min(old_xp) + 0.5 - xc
min_y = na.min(old_yp) + 0.5 - yc
min_z = na.min(old_zp) + 0.5 - zc
line = "%f %f %f %f %f %f %f\n" % (red, min_x-fudge, min_y-fudge, min_z-fudge,
max_x+fudge, max_y+fudge, max_z+fudge)
fp.write(line)
gc.collect()
# For z=0.
xp = haloes[hID]["particle_position_x"]
yp = haloes[hID]["particle_position_y"]
zp = haloes[hID]["particle_position_z"]
line = "%f %f %f %f %f %f %f\n" % \
(0.0,
na.min(xp) + 0.5 - xc - fudge,
na.min(yp) + 0.5 - yc - fudge,
na.min(zp) + 0.5 - zc - fudge,
na.max(xp) + 0.5 - xc + fudge,
na.max(yp) + 0.5 - yc + fudge,
na.max(zp) + 0.5 - zc + fudge)
fp.write(line)
fp.close()