From a9e040540b6f4afa1480b0c5c16c2860bcddf1b4 Mon Sep 17 00:00:00 2001
From: "P.M. Sutter" <sutter@iap.fr>
Date: Sun, 9 Feb 2014 06:11:06 -0600
Subject: [PATCH] fixed bug in handling subhalos

---
 .../pipeline_source/prepareCatalogs.in.py     | 44 +++++++++++++------
 .../void_python_tools/backend/launchers.py    | 14 +++---
 2 files changed, 40 insertions(+), 18 deletions(-)

diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py
index ba7cd42..1a805cd 100644
--- a/python_tools/pipeline_source/prepareCatalogs.in.py
+++ b/python_tools/pipeline_source/prepareCatalogs.in.py
@@ -264,6 +264,19 @@ newSample.addStack(0.0, 5.0, 85, 90, False, False, rescaleMode="rv")
 newSample.addStack(0.0, 5.0, 90, 95, False, False, rescaleMode="rv")
 newSample.addStack(0.0, 5.0, 95, 100, False, False, rescaleMode="rv")
   """
+  elif stackMode == "sample_thick":
+    stackInfo = """
+# {zMin}, {zMax}, {minRadius}
+newSample.addStack({zMin}, {zMax}, 10, 20, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 20, 30, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 30, 40, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 40, 50, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 50, 60, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 60, 70, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 70, 80, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 80, 90, True, False, rescaleMode="rmax")
+newSample.addStack({zMin}, {zMax}, 90, 100, True, False, rescaleMode="rmax")
+  """
 
   elif stackMode == "log":
     stackInfo = ""
@@ -279,15 +292,15 @@ newSample.addStack(0.0, 5.0, 95, 100, False, False, rescaleMode="rv")
       
       rStart = rEnd
 
-#  elif stackMode == "auto":
-#    stackInfo = """
-#newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+2, True, False, rescaleMode="rv")
-#newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+4, True, False, rescaleMode="rv")
-#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False, rescaleMode="rv")
-#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False, rescaleMode="rv")
-#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False, rescaleMode="rv")
-#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False, rescaleMode="rv")
-#               """
+  elif stackMode == "auto":
+    stackInfo = """
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+2, True, False, rescaleMode="rv")
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}  , 2*{minRadius}+4, True, False, rescaleMode="rv")
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False, rescaleMode="rv")
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False, rescaleMode="rv")
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False, rescaleMode="rv")
+newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False, rescaleMode="rv")
+               """
   else:
     stackInfo = """
 # {zMin}, {zMax}, {minRadius}
@@ -840,20 +853,25 @@ if (args.halos or args.all) and haloFileBase != "":
 
         # do the check again since now we've filtered out subhalos
         numPart = 0
-        for line in enumerate(tempFile):
+        fileTemp = open(tempFile, 'r')
+        for line in fileTemp:
           numPart += 1
- 
+        fileTemp.close()
+
         actualDen = 1.*numPart / lbox**3
         keepFraction = haloDen / actualDen
         if numPart < numPartExpect:
           print " ERROR: not enough galaxies to support that density! Maximum is %g" % (1.*numPart / lbox**3)
           exit(-1)
-        
+
         numKept = 0
+        inFile = open(dataFile, 'r')
+        outFile = open(outFileName, 'a')
         for (iLine,line) in enumerate(inFile):
           if np.random.uniform() > keepFraction: continue
           outFile.write(line)
-          numKept += 1 
+          numKept += 1
+        inFile.close()
 
         outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n")
         outFile.close()
diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py
index 2fa66fc..9d9b393 100644
--- a/python_tools/void_python_tools/backend/launchers.py
+++ b/python_tools/void_python_tools/backend/launchers.py
@@ -43,6 +43,8 @@ import void_python_tools.apTools as vp
 NetCDFFile = Dataset
 ncFloat = 'f8' # Double precision
 
+LIGHT_SPEED = 299792.458
+
 # -----------------------------------------------------------------------------
 def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, 
                    zobovDir=None, figDir=None, logFile=None, useComoving=False,
@@ -561,6 +563,7 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir,
 def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
                 voidDir=None, freshStack=True, runSuffix=None,
                 zobovDir=None, useLightCone=None, useComoving=None,
+                omegaM=None,
                 INCOHERENT=False, ranSeed=None, summaryFile=None, 
                 continueRun=None, dataType=None, prefixRun="",
                 idList=None, rescaleOverride=None):
@@ -589,12 +592,13 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
   zMin = max(sample.zRange[0], stack.zMin)
   zMax = min(sample.zRange[1], stack.zMax)
   if useComoving or not useLightCone:
-    zMin = 3000.*vp.angularDiameter(zMin, Om=0.27)
-    zMax = 3000.*vp.angularDiameter(zMax, Om=0.27)
-    print min(sample.zRange[1],stack.zMax)*3000, zMax
+    zMin = LIGHT_SPEED/100.*vp.angularDiameter(zMin, Om=omegaM)
+    zMax = LIGHT_SPEED/100.*vp.angularDiameter(zMax, Om=omegaM)
+    print min(sample.zRange[1],stack.zMax)*LIGHT_SPEED/100., zMax
+    print max(sample.zRange[0],stack.zMin)*LIGHT_SPEED/100., zMin
   else:
-    zMin *= 3000
-    zMax *= 3000
+    zMin *= LIGHT_SPEED/100.
+    zMax *= LIGHT_SPEED/100.
 
   if dataType == "observation":
     obsFlag = "observation"