[yt-users] simple calculation question

Stephen Skory stephenskory at yahoo.com
Thu Aug 12 13:24:51 PDT 2010


Hi Robert,


>How would I go about doing this in yt given the halo finder outputs.  Thanks for 
>
>the assistance!


I was about to step through the whole process for you, but I realized that an 
example script would be far more useful. I can't quite parse the units of your 
box, so distances are in code units below and you can change that if you need 
to. Remember that the halo velocity output of the HopAnalysis.out file is in 
cgs, not code units.

This double loop is not very scalable. If you are doing large lists of halos, 
you might consider using a kD-tree which will speed things up significantly.

---------

# Essential modules
from yt.mods import *
from yt.math_utils import *

# The halo finder output file
infile = "HopAnalysis.out"
# Halo separation distance in code units.
dmin = 0.01

# Define some lists to store our data. Depending on application,
# it is usally best to convert these into Numpy arrays, but it's not
# required for this.
Xctr, Yctr, Zctr, vx, vy, vz = [], [], [], [], [], []

# Now we read it in.
data = file(infile, "r")
for line in data:
    # Skip the header lines
    if "#" in line: continue
    # Splits the line up by spaces
    line = lines.split()
    # Store the data using the conventional columns for yt-halofinder output
    Xctr.append(float(line[7]))
    Yctr.append(float(line[8]))
    Zctr.append(float(line[9]))
    xv.append(float(line[10]))
    yv.append(float(line[11]))
    zv.append(float(line[12]))
# Close the text file.
data.close()

# Now the double loop.
for i in range(0,len(Xctr)):
    for j in range(i+1, len(Xctr)):
        # This function periodic_dist comes from math_utils,
        # and d is in code units because every input is in code units.
        d = periodic_dist([Xctr[i], Yctr[i], Zctr[i]],
            [Xctr[j], Yctr[j], Zctr[j]], 1.0)
        if d < dmin:
            dvx = vx[i] - vx[j]
            dvy = vy[i] - vy[j]
            dvz = vz[i] - vz[j]
            
            vrel = math.sqrt(dvx*dvx + dvy*dvy + dvz*dvz)



 _______________________________________________________
sskory at physics.ucsd.edu           o__  Stephen Skory
http://physics.ucsd.edu/~sskory/ _.>/ _Graduate Student
________________________________(_)_\(_)_______________




More information about the yt-users mailing list