[yt-users] simple calculation question

[John Wise]

You can also use some numpy "array magic" to speed the double loop up.

http://mail.scipy.org/pipermail/numpy-discussion/2007-April/027203.html

Here's a kD-tree implementation that I've used before when searching for nearest neighbors.

http://pypi.python.org/pypi/scikits.ann/0.2.dev-r803

John

On 12 Aug 2010, at 16:24, Stephen Skory wrote:

> 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
> ________________________________(_)_\(_)_______________
> 
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org