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

[britton at wrangler.dreamhost.com]

Author: britton
Date: Fri May  9 13:58:46 2008
New Revision: 453
URL: http://yt.spacepope.org/changeset/453

Log:
Added IsBound derived quantity which returns TRUE if group of cells is 
gravitationally bound.

9 May, 2008

Britton Smith



Modified:
   trunk/yt/lagos/DerivedQuantities.py

Modified: trunk/yt/lagos/DerivedQuantities.py
==============================================================================
--- trunk/yt/lagos/DerivedQuantities.py	(original)
+++ trunk/yt/lagos/DerivedQuantities.py	Fri May  9 13:58:46 2008
@@ -28,6 +28,7 @@
 """
 
 from yt.lagos import *
+from yt.funcs import get_pbar
 
 quantity_info = {}
 
@@ -152,3 +153,34 @@
     return spin
 add_quantity("BaryonSpinParameter", function=_BaryonSpinParameter,
              combine_function=_combBaryonSpinParameter, n_ret=4)
+
+def _IsBound(data):
+    # Kinetic energy
+    bv_x,bv_y,bv_z = data.quantities["BulkVelocity"]()
+    v_x2 = (data["x-velocity"] - bv_x)**2
+    v_y2 = (data["y-velocity"] - bv_y)**2
+    v_z2 = (data["z-velocity"] - bv_z)**2
+    kinetic = 0.5 * (data["CellMass"] * (v_x2 + v_y2 + v_z2)).sum()
+
+    # Gravitational potential energy
+    G = 6.67e-8 # 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 ", (na.array(range(total_cells)).sum())) # Progress bar
+    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 * data['CellMass'][q] * pot.sum()
+        cells_done += (total_cells - q - 1)
+        pb.update(cells_done)
+    potential *= G
+    pb.finish()
+
+    return (kinetic < potential),1.0
+def _combIsBound(data,bound,weight):
+    return bound
+add_quantity("IsBound",function=_IsBound,combine_function=_combIsBound,n_ret=2)