[Yt-svn] yt-commit r472 - trunk/yt/lagos

[mturk at wrangler.dreamhost.com]

Author: mturk
Date: Wed May 14 14:06:39 2008
New Revision: 472
URL: http://yt.spacepope.org/changeset/472

Log:
Added a brief, non-pretty progress bar to IsBound's C function call.  (Couldn't
get it to work by calling Python code.)  Made C code the default.



Modified:
   trunk/yt/lagos/DerivedQuantities.py
   trunk/yt/lagos/PointCombine.c

Modified: trunk/yt/lagos/DerivedQuantities.py
==============================================================================
--- trunk/yt/lagos/DerivedQuantities.py	(original)
+++ trunk/yt/lagos/DerivedQuantities.py	Wed May 14 14:06:39 2008
@@ -155,7 +155,7 @@
 add_quantity("BaryonSpinParameter", function=_BaryonSpinParameter,
              combine_function=_combBaryonSpinParameter, n_ret=4)
 
-def _IsBound(data, truncate = True, include_thermal_energy = False, use_c = False):
+def _IsBound(data, truncate = True, include_thermal_energy = False):
     # Kinetic energy
     bv_x,bv_y,bv_z = data.quantities["BulkVelocity"]()
     kinetic = 0.5 * (data["CellMass"] * (
@@ -170,35 +170,10 @@
     # We only divide once here because we have velocity in cgs, but radius is
     # in code.
     G = 6.67e-8 / data.convert("cm") # cm^3 g^-1 s^-2
-    # the slow way (someone get clever!)
-    total_cells = len(data['x'])
-    cells_done = 0
-    potential = 0.0
-    pb = get_pbar("Calculating gravitational potential ", (0.5*(total_cells**2 - total_cells)))
-    import time
-    t1 = time.time()
-    if use_c:
-        pot = PointCombine.FindBindingEnergy(data["CellMass"],
-                                       data["x"], data["y"], data["z"],
-                                       truncate, kinetic/(2*G))
-        potential += 2*G*pot
-    else:
-        t1 = time.time()
-        for q in xrange(total_cells-1):
-            pot = data['CellMass'][(q+1):total_cells] / \
-                (((data['x'][(q+1):total_cells]-data['x'][q])**2 + \
-                  (data['y'][(q+1):total_cells]-data['y'][q])**2 + \
-                  (data['z'][(q+1):total_cells]-data['z'][q])**2)**(0.5))
-            potential += 2 * G * data['CellMass'][q] * pot.sum()
-            cells_done += (total_cells - q - 1)
-            pb.update(cells_done)
-            if truncate and (potential > kinetic):
-                pb.finish()
-                break
-    pb.finish()
-    print "Took %0.9e" % (time.time() - t1)
-    #if truncate: return [potential > kinetic]
-    return [(potential / kinetic)]
+    pot = 2*G*PointCombine.FindBindingEnergy(data["CellMass"],
+                                  data['x'],data['y'],data['z'],
+                                  truncate, kinetic/(2*G))
+    return [(pot / kinetic)]
 def _combIsBound(data,bound):
     return bound
 add_quantity("IsBound",function=_IsBound,combine_function=_combIsBound,n_ret=1,

Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c	(original)
+++ trunk/yt/lagos/PointCombine.c	Wed May 14 14:06:39 2008
@@ -844,19 +844,21 @@
     }
 
     /* Do the work here. */
-    int q_outer, q_inner;
+    int q_outer, q_inner, n_q = PyArray_SIZE(mass);
     double this_potential, total_potential;
     total_potential = 0;
     npy_float64 mass_o, x_o, y_o, z_o;
     npy_float64 mass_i, x_i, y_i, z_i;
 
-    for (q_outer = 0; q_outer < PyArray_SIZE(mass) ; q_outer++) {
+    for (q_outer = 0; q_outer < n_q ; q_outer++) {
         this_potential = 0;
         mass_o = *(npy_float64*) PyArray_GETPTR1(mass, q_outer);
         x_o = *(npy_float64*) PyArray_GETPTR1(x, q_outer);
         y_o = *(npy_float64*) PyArray_GETPTR1(y, q_outer);
         z_o = *(npy_float64*) PyArray_GETPTR1(z, q_outer);
-        for (q_inner = q_outer+1; q_inner < PyArray_SIZE(mass); q_inner++) {
+        fprintf(stderr,"Calculating Potential: % 12i / % 12i\r", q_outer, n_q);
+        fflush(stdout); fflush(stderr);
+        for (q_inner = q_outer+1; q_inner < n_q; q_inner++) {
             mass_i = *(npy_float64*) PyArray_GETPTR1(mass, q_inner);
             x_i = *(npy_float64*) PyArray_GETPTR1(x, q_inner);
             y_i = *(npy_float64*) PyArray_GETPTR1(y, q_inner);
@@ -868,9 +870,12 @@
         }
         total_potential += this_potential;
         if ((truncate == 1) && (total_potential > kinetic_energy)){
+            fprintf(stderr, "Truncating!\r");
             break;
         }
     }
+    fprintf(stderr,"\n");
+    fflush(stdout); fflush(stderr);
 
     Py_DECREF(mass);
     Py_DECREF(x);