rename helper function from angularDiameter to comovingDistance since that's what it actually is

python_source/backend/cosmologyTools.py

@@ -25,7 +25,7 @@ import scipy.integrate as integrate
 import os
 from backend import *
 
-__all__=['expansion', 'angularDiameter', 'aveExpansion']
+__all__=['expansion', 'comovingDistance', 'aveExpansion']
 
 # returns 1/E(z) for the given cosmology
 def expansion(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
@@ -38,8 +38,8 @@ def expansion(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
 def eosDE(z, w0 = -1.0, wa = 0.0):
   return w0 + wa*z/(1+z)
   
-# returns D_A(z) for the given cosmology
-def angularDiameter(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
+# returns the comoving distance for the given cosmology
+def comovingDistance(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0):
   da = integrate.quad(expansion, 0.0, z, args=(Om, Ot, w0, wa))[0]
   return da
   

python_source/backend/surveyTools.py

@@ -1,5 +1,5 @@
 #+
-#   VIDE -- Void IDentification and Examination -- ./python_tools/vide/apTools/chi2/cosmologyTools.py
+#   VIDE -- Void IDentification and Examination
 #   Copyright (C) 2010-2014 Guilhem Lavaux
 #   Copyright (C) 2011-2014 P. M. Sutter
 #
@@ -38,8 +38,8 @@ def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion,
   area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi
 
   if useComoving:
-    zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=omegaM)
-    zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=omegaM)
+    zmin = LIGHT_SPEED/100.*comovingDistance(zmin, Om=omegaM)
+    zmax = LIGHT_SPEED/100.*comovingDistance(zmax, Om=omegaM)
   else:
     zmin = zmin * LIGHT_SPEED/100.
     zmax = zmax * LIGHT_SPEED/100. 
@@ -95,11 +95,11 @@ def getNside(maskFile):
 # helper function to convert RA,dec to phi,theta
 def convertAngle(RA, Dec):
 
-    phi   = np.pi/180.*RA
-    theta = Dec*np.pi/180.
-    theta = np.pi/2. - Dec*np.pi/180.
+  phi   = np.pi/180.*RA
+  theta = Dec*np.pi/180.
+  theta = np.pi/2. - Dec*np.pi/180.
 
-  return phi, theta
+  return (phi, theta)
 
 # -----------------------------------------------------------------------------
 # computes the mask from a given galaxy datafile and writes it to a file
@@ -110,16 +110,16 @@ def figureOutMask(galFile, nside, outMaskFile):
 
   for line in open(galFile):
     line = line.split()
-    RA  = np.float(line[3])
-    Dec = np.float(line[4])
-    z   = np.float(line[5])
+    RA  = float(line[3])
+    Dec = float(line[4])
+    z   = float(line[5])
 
     phi, theta = convertAngle(RA, Dec)
 
     pix = healpy.ang2pix(nside, theta, phi)
     mask[pix] = 1.
  
-    healpy.write_map(outMaskFile, mask)
+  healpy.write_map(outMaskFile, mask, overwrite=True)
 
   return mask
 
@@ -130,8 +130,8 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, edgeMaskFile,
                      zmin, zmax, omegaM, useComoving=False):
 
   if useComoving:
-    zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=omegaM)
-    zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=omegaM)
+    zmin = LIGHT_SPEED/100.*comovingDistance(zmin, Om=omegaM)
+    zmax = LIGHT_SPEED/100.*comovingDistance(zmax, Om=omegaM)
   else:
     zmin *= LIGHT_SPEED/100.
     zmax *= LIGHT_SPEED/100.
@@ -146,31 +146,37 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, edgeMaskFile,
 
   for line in open(galFile):
     line = line.split()
-    RA  = np.float(line[3])
-    Dec = np.float(line[4])
-    z   = np.float(line[5])
+    RA  = float(line[3])
+    Dec = float(line[4])
+    z   = float(line[5])
 
     if useComoving:
-      z = LIGHT_SPEED/100.*angularDiameter(z, Om=omegaM)
+      z = LIGHT_SPEED/100.*comovingDistance(z, Om=omegaM)
     else:
       z *= LIGHT_SPEED/100.
 
     phi, theta = convertAngle(RA, Dec)
 
-    isOnEdge = False
-
     # check the mask edges
     neighbors = healpy.get_all_neighbors(nside, theta, phi)
-    isOnEdge = any(p == 0 for p in neighbors):
+    isOnMaskEdge = any(p == 0 for p in neighbors)
   
     # check the redshift boundaries
+    tol = 0.01 # TODO: made this user-adjustable
+    zBuffer = (zmax-zmin)*tol
+    isOnHighZEdge = (z >= zmax-buffer)
+    isOnLowZEdge  = (z <= zmin+buffer)
 
-    if isOnEdge:
+    if isOnMaskEdge:
       edgeFile.write("1")
+    elif isOnHighZEdge:
+      edgeFile.write("2")
+    elif isOnLowZEdge:
+      edgeFile.write("3")
     else:
       edgeFile.write("0")
 
   edgeFile.close()
-  healpy.write_map(edgeMaskFile, edgeMask)
+  healpy.write_map(edgeMaskFile, edgeMask, overwrite=True)
 
   return

python_source/vide_pipeline/vide_prepare_simulation.py

@@ -226,16 +226,16 @@ dataSampleList.append(newSample)
 
     zBox = float(redshift)
     Om = float(omegaM)
-    zBoxMpc = LIGHT_SPEED/100.*vp.angularDiameter(zBox, Om=Om)
+    zBoxMpc = LIGHT_SPEED/100.*vp.comovingDistance(zBox, Om=Om)
     boxMaxMpc = zBoxMpc + lbox
 
     # converter from redshift to comoving distance   
     zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000)
     zVsDX = np.zeros(len(zVsDY))
     for i in range(len(zVsDY)):
-      zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om)
+      zVsDX[i] = vp.comovingDistance(zVsDY[i], Om=Om)
 
-    boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \
+    boxWidthZ = np.interp(vp.comovingDistance(zBox,Om=Om)+100. / \
                 LIGHT_SPEED*lbox, zVsDX, zVsDY)-zBox
 
     for iSlice in range(numSlices):