[Yt-svn] yt-commit r1149 - trunk/yt/extensions

[britton at wrangler.dreamhost.com]

Author: britton
Date: Mon Jan 26 09:54:23 2009
New Revision: 1149
URL: http://yt.spacepope.org/changeset/1149

Log:
Added some catches for halos that are not resolved enough for radial profiles.  
I lowered the minimum profile radius from 4*dx_min to 2*dx_min to be slightly 
more generous with small halos.  The variable PROFILE_RADIUS_THRESHOLD sets the 
minimum ratio of r_max to r_min for profiling is currently set to 2.  I've 
tried this and it seems to not fail for values of r_max/r_min near the cutoff 
value.  In case the code gets past this check but still produces invalid 
profiles, I've put another check in place to return zeros for the virial 
quantities if no valid profile data exist.


Modified:
   trunk/yt/extensions/HaloProfiler.py

Modified: trunk/yt/extensions/HaloProfiler.py
==============================================================================
--- trunk/yt/extensions/HaloProfiler.py	(original)
+++ trunk/yt/extensions/HaloProfiler.py	Mon Jan 26 09:54:23 2009
@@ -31,6 +31,8 @@
 import os
 import tables as h5
 
+PROFILE_RADIUS_THRESHOLD = 2
+
 class HaloProfiler(object):
     def __init__(self,dataset,HaloProfilerParameterFile,halos='multiple',radius=0.1):
         self.dataset = dataset
@@ -94,9 +96,14 @@
         for q,halo in enumerate(self.hopHalos):
             filename = "%s/Halo_%04d_profile.dat" % (outputDir,q)
 
+            r_min = 2*self.pf.h.get_smallest_dx() * self.pf['mpc']
+            if (halo['r_max'] / r_min < PROFILE_RADIUS_THRESHOLD):
+                mylog.error("Skipping halo with r_max / r_min = %f." % (halo['r_max']/r_min))
+                continue
+
             sphere = self.pf.h.sphere(halo['center'],halo['r_max']/self.pf.units['mpc'])
             profile = lagos.BinnedProfile1D(sphere,self.haloProfilerParameters['n_bins'],"RadiusMpc",
-                                            (4*self.pf.h.get_smallest_dx() * self.pf['mpc']),halo['r_max'],
+                                            r_min,halo['r_max'],
                                             log_space=True, lazy_reader=True)
             for field in self.profileFields.keys():
                 profile.add_fields(field,weight=self.profileFields[field][0],
@@ -259,12 +266,6 @@
         for field in fields:
             temp_profile[field] = []
 
-        virial = {}
-        if (len(profile[overdensity_field]) < 2):
-            for field in fields:
-                virial[field] = 0.0
-            return virial
-
         for q in range(len(profile[overdensity_field])):
             good = True
             if (profile[overdensity_field][q] != profile[overdensity_field][q]):
@@ -279,6 +280,13 @@
                 for field in fields:
                     temp_profile[field].append(profile[field][q])
 
+        virial = {}
+        if (len(overDensity) < 2):
+            mylog.error("Skipping halo with no valid points in profile.")
+            for field in fields:
+                virial[field] = 0.0
+            return virial
+
         if (overDensity[1] <= self.haloProfilerParameters['VirialOverdensity']):
             index = 0
         elif (overDensity[-1] >= self.haloProfilerParameters['VirialOverdensity']):