From 7ba3b2a98d1e26ce5a402cb8ea12c8c90eb169e4 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 23 Aug 2013 13:13:55 -0500 Subject: [PATCH] prunevoids now can handle coordinate changes to cosmoving space in observations --- c_tools/stacking/pruneVoids.cpp | 75 +++++++++++++++++-- pipeline/generateCatalog.py | 3 +- .../void_python_tools/backend/launchers.py | 8 +- 3 files changed, 79 insertions(+), 7 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index c857496..6f8f60f 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -43,6 +43,8 @@ #include "assert.h" #include "voidTree.hpp" #include "loadZobov.hpp" +#include +#include #define LIGHT_SPEED 299792.458 #define MPC2Z 100./LIGHT_SPEED @@ -82,6 +84,32 @@ typedef struct voidStruct { float ellip; } VOID; +// this defines the expansion function that we will integrate +// Laveaux & Wandelt (2012) Eq. 24 +struct my_expan_params { double Om; double w0; double wa; }; +double expanFun (double z, void * p) { + struct my_expan_params * params = (struct my_expan_params *)p; + double Om = (params->Om); + double w0 = (params->w0); + double wa = (params->wa); + + //const double h0 = 1.0; + const double h0 = 0.71; + double ez; + + double wz = w0 + wa*z/(1+z); + + ez = Om*pow(1+z,3) + (1.-Om); + //ez = Om*pow(1+z,3) + pow(h0,2) * (1.-Om)*pow(1+z,3+3*wz); + + ez = sqrt(ez); + //ez = sqrt(ez)/h0; + + ez = 1./ez; + + return ez; +} + void openFiles(string outputDir, string sampleName, string prefix, string dataPortion, int mockIndex, int numKept, @@ -127,7 +155,34 @@ int main(int argc, char **argv) { return 1; } - int i, p, p2, numPartTot, numZonesTot, dummy, iVoid, iZ; + // initialize cosmology integrator and interpolator + gsl_function expanF; + expanF.function = &expanFun; + struct my_expan_params expanParams; + expanParams.Om = 0.27; + expanParams.w0 = -1.0; + expanParams.wa = 0.0; + expanF.params = &expanParams; + double result, error; + size_t nEval; + + int iZ, numZ = 4000; + double maxZ = 2.0, z, *dL, *redshifts; + dL = (double *) malloc(numZ * sizeof(double)); + redshifts = (double *) malloc(numZ * sizeof(double)); + for (iZ = 0; iZ < numZ; iZ++) { + z = iZ * maxZ/numZ; + gsl_integration_qng(&expanF, 0.0, z, 1.e-6, 1.e-6, &result, &error, &nEval); + dL[iZ] = result*LIGHT_SPEED/100.; + //printf("HERE %e %e\n", z, dL[iZ]); + redshifts[iZ] = z; + } + gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, 4000); + gsl_interp_init(interp, dL, redshifts, numZ); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + + + int i, p, p2, numPartTot, numZonesTot, dummy, iVoid; int numVoids, mockIndex; double tolerance; FILE *fp; @@ -541,11 +596,21 @@ int main(int argc, char **argv) { pow(voids[iVoid].barycenter[1] - boxLen[1]/2.,2) + pow(voids[iVoid].barycenter[2] - boxLen[2]/2.,2)); voids[iVoid].redshiftInMpc = voids[iVoid].redshiftInMpc; - redshift = voids[iVoid].redshiftInMpc; - nearestEdge = fabs(redshift-args.zMax_arg*LIGHT_SPEED/100.); + + + if (args.useLCDM_flag) { + redshift = gsl_interp_eval(interp, dL, redshifts, + voids[iVoid].redshiftInMpc, acc); + //printf("HELLO %e %e\n", redshift, args.zMax_arg); + nearestEdge = fabs((redshift-args.zMax_arg)*LIGHT_SPEED/100.); + voids[iVoid].redshift = redshift; + } else { + redshift = voids[iVoid].redshiftInMpc; + nearestEdge = fabs(redshift-args.zMax_arg*LIGHT_SPEED/100.); + voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.; + } //nearestEdge = fmin(fabs(redshift-args.zMin_arg*LIGHT_SPEED/100.), // fabs(redshift-args.zMax_arg*LIGHT_SPEED/100.)); - voids[iVoid].redshift = voids[iVoid].redshiftInMpc/LIGHT_SPEED*100.; } else { @@ -844,7 +909,7 @@ void openFiles(string outputDir, string sampleName, fprintf(*fpSkyPositions, "# RA, dec, redshift, radius (Mpc/h), void ID\n"); *fpShapes = fopen((outputDir+"/"+prefix+"shapes_"+dataPortion+"_"+sampleName).c_str(), "w"); - fprintf(*fpShapes, "# void ID, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n"); + fprintf(*fpShapes, "# void ID, ellip, eig(1), eig(2), eig(3), eigv(1)-x, eiv(1)-y, eigv(1)-z, eigv(2)-x, eigv(2)-y, eigv(2)-z, eigv(3)-x, eigv(3)-y, eigv(3)-z\n"); } // end openFiles diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 87d50a5..741a709 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -137,7 +137,8 @@ for sample in dataSampleList: PRUNE_PATH = CTOOLS_PATH+"/stacking/pruneVoids" launchPrune(sample, PRUNE_PATH, - logFile=logFile, zobovDir=zobovDir, continueRun=continueRun) + logFile=logFile, zobovDir=zobovDir, + useLCDM=useLCDM, continueRun=continueRun) # ------------------------------------------------------------------------- if (startCatalogStage <= 4) and (endCatalogStage >= 4): diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index cc72294..de7e6ee 100644 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -301,7 +301,7 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, # ----------------------------------------------------------------------------- def launchPrune(sample, binPath, summaryFile=None, logFile=None, zobovDir=None, - continueRun=None): + continueRun=None, useLCDM=False): sampleName = sample.fullName @@ -329,6 +329,11 @@ def launchPrune(sample, binPath, periodicLine = " --periodic='" + getPeriodic(sample) + "'" + if useLCDM: + useLCDMFlag = "--useLCDM" + else: + useLCDMFlag = "" + if not (continueRun and (jobSuccessful(logFile, "NetCDF: Not a valid ID\n") \ or jobSuccessful(logFile, "Done!\n"))): cmd = binPath @@ -348,6 +353,7 @@ def launchPrune(sample, binPath, cmd += " --numVoids=" + str(numVoids) cmd += observationLine cmd += periodicLine + cmd += useLCDMFlag cmd += " --outputDir=" + zobovDir cmd += " --sampleName=" + str(sampleName) cmd += " &> " + logFile