simple calculation question
Hi yt users, I'm still very new at python syntax so please forgive me =) I've ran a few n-body simulations and I'd like to simply calculate the pairwise velocity between the halo pairs. I simply need to read in the x,y,z positions & vx,vy,vz velocities from the halo finder output, do some calculations, and output the result. I've done this in C but I'd love to learn how to do this in yt. Here's my calculation in C: for(i=0; i<NUM; i++) { for(j=i+1; j<NUM; j++) { dXctr = periodic(Xctr[i] - Xctr[j]) / 1000.; /* in Mpc/h */ dYctr = periodic(Yctr[i] - Yctr[j]) / 1000.; dZctr = periodic(Zctr[i] - Zctr[j]) / 1000.; d = sqrt(dXctr*dXctr + dYctr*dYctr + dZctr*dZctr); /* in Mpc/h */ if(d < 10.) /* less than 10 Mpc/h */ { dvx = vx[i] - vx[j]; dvy = vy[i] - vy[j]; dvz = vz[i] - vz[j]; vrel = sqrt(dvx*dvx + dvy*dvy + dvz*dvz); } } How would I go about doing this in yt given the halo finder outputs. Thanks for the assistance! -Robert Thompson
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@physics.ucsd.edu o__ Stephen Skory http://physics.ucsd.edu/~sskory/ _.>/ _Graduate Student ________________________________(_)_\(_)_______________
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@physics.ucsd.edu o__ Stephen Skory http://physics.ucsd.edu/~sskory/ _.>/ _Graduate Student ________________________________(_)_\(_)_______________
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Sorry for the delayed response! Thank you for your help that example script helped a ton. It'll take some time before I dive into the numpy array magic to speed up large calculations. -Robert On 08/12/2010 01:41 PM, John Wise wrote:
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@physics.ucsd.edu o__ Stephen Skory http://physics.ucsd.edu/~sskory/ _.>/ _Graduate Student ________________________________(_)_\(_)_______________
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
participants (3)
-
John Wise -
Robert Thompson -
Stephen Skory