diff --git a/c_source/prep/prepObservation.cpp b/c_source/prep/prepObservation.cpp index ec72a0f..073898f 100644 --- a/c_source/prep/prepObservation.cpp +++ b/c_source/prep/prepObservation.cpp @@ -113,7 +113,7 @@ void loadData(const string& fname, NYU_VData & data) // Maubert - avoid double counting of last particle in array if EOF is after a blank line if (!f) { - continue; + continue; } // End Maubert d.uniqueID = d.index; @@ -273,7 +273,7 @@ void generateSurfaceMask(prepObservation_info& args , output_data.Lmax = Rmax; -// PMS - write a small text file with galaxy position +// PMS - write a small text file with galaxy position (for diagnostic purposes) FILE *fp; fp = fopen("galaxies.txt", "w"); for (int i = 0; i < data.size(); i++) { @@ -289,6 +289,9 @@ void generateSurfaceMask(prepObservation_info& args , cout << format("Rmax is %g, surface volume is %g") % (Rmax/100) % (volume/(4*M_PI)) << endl; volume *= Rmax*Rmax*Rmax/3/1e6; numToInsert = (int)floor(volume*args.density_fake_arg); + // TEST NOT USING MOCK PARTICLES + numToIntsert = 0 + // END TEST cout << format("3d volume to fill: %g (Mpc/h)^3") % volume << endl; cout << format("Will insert %d particles") % numToInsert << endl; @@ -325,7 +328,7 @@ void generateSurfaceMask(prepObservation_info& args , } while (!stop_here); -// PMS : write mock galaxies to a small file for plotting +// PMS : write mock galaxies to a small file for diagnostic purposes fprintf(fp, "%e %e %e\n", (p.xyz[0]), (p.xyz[1]), @@ -344,7 +347,7 @@ void generateSurfaceMask(prepObservation_info& args , fclose(fp); // PMS - // TEST - insert mock galaxies along box edge + // TEST - insert mock galaxies along box edge - this is for tesselation safety fp = fopen("mock_boundary.txt", "w"); double dx[3]; dx[0] = output_data.box[0][1] - output_data.box[0][0]; @@ -410,8 +413,9 @@ void generateSurfaceMask(prepObservation_info& args , rmax = args.zMax_arg * LIGHT_SPEED; } - //for (int q = 0; q < 0; q++) { - for (int q = 0; q < full_mask_list.size(); q++) { + // TEST NOT USING BOUNDARY PARTICLES + for (int q = 0; q < 0; q++) { + //for (int q = 0; q < full_mask_list.size(); q++) { vec3 v = mask.pix2vec(full_mask_list[q]); Position p; diff --git a/examples/example_observation/example_observation.py b/examples/example_observation/example_observation.py index 219cb1a..2c7da8d 100644 --- a/examples/example_observation/example_observation.py +++ b/examples/example_observation/example_observation.py @@ -83,7 +83,7 @@ newSample = Sample( volumeLimited = True, # HEALpix mask file - set to None to auto-compute - maskFile = None, + maskFile = "", #maskFile = inputDataDir+"/example_observation_mask.fits", # if maskFile blank, desired resolution for HEALpix diff --git a/python_source/backend/launchers.py b/python_source/backend/launchers.py index 7ddcf29..a171c4b 100644 --- a/python_source/backend/launchers.py +++ b/python_source/backend/launchers.py @@ -47,25 +47,36 @@ LIGHT_SPEED = 299792.458 # ----------------------------------------------------------------------------- def launchPrep(sample, binPath, workDir=None, inputDataDir=None, - zobovDir=None, figDir=None, logFile=None, useComoving=False, + outputDir=None, figDir=None, logFile=None, useComoving=False, continueRun=None, regenerate=False): if sample.dataType == "observation": sampleName = sample.fullName if regenerate: - inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par" - outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName + inputParameterFlag = "inputParameter " + outputDir + + "/zobov_slice_"+sampleName+".par" + outputFile = outputDir + "/regenerated_zobov_slice_" + sampleName else: inputParameterFlag = "" - outputFile = zobovDir+"/zobov_slice_" + sampleName + outputFile = outputDir + "/zobov_slice_" + sampleName if sample.dataFile == "": datafile = inputDataDir+"/"+sampleName else: datafile = inputDataDir+"/"+sample.dataFile - maskFile = sample.maskFile + if sample.maskFile == "": + maskFile = outputDir + "/constructed_mask.fits" + figureOutMask(dataFile, sample.nsideForMask, maskFile) + else: + maskFile = sample.maskFile + + edgeGalFile = outputDir + "/galaxy_edge_flags.txt" + edgeMaskFile = outputDir + "/mask_edge_map.fits" + findEdgeGalaxies(dataFile, maskFile, edgeGalFile, edgeMaskFile, + sample.zBoundary[0], sample.zBoundary[1], sample.omegaM, + useComoving) if useComoving: useComovingFlag = "useComoving" @@ -84,7 +95,7 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, %s omegaM %g """ % (datafile, maskFile, outputFile, - zobovDir+"/zobov_slice_"+sampleName+".par", + outputDir+"/zobov_slice_"+sampleName+".par", sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity, useComovingFlag, inputParameterFlag, sample.omegaM) @@ -108,21 +119,21 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, if os.access(parmFile, os.F_OK): os.unlink(parmFile) if os.access("contour_map.fits", os.F_OK): - os.system("mv %s %s" % ("contour_map.fits", zobovDir)) - os.system("mv %s %s" % ("mask_map.fits", zobovDir)) + os.system("mv %s %s" % ("contour_map.fits", outputDir)) + os.system("mv %s %s" % ("mask_map.fits", outputDir)) if os.access("comoving_distance.txt", os.F_OK): - os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) + os.system("mv %s %s" % ("comoving_distance.txt", outputDir)) if os.access("mask_index.txt", os.F_OK): - os.system("mv %s %s" % ("mask_index.txt", zobovDir)) - os.system("mv %s %s" % ("total_particles.txt", zobovDir)) + os.system("mv %s %s" % ("mask_index.txt", outputDir)) + os.system("mv %s %s" % ("total_particles.txt", outputDir)) if os.access("galaxies.txt", os.F_OK): - os.system("mv %s %s" % ("galaxies.txt", zobovDir)) - os.system("mv %s %s" % ("mock_galaxies.txt", zobovDir)) - os.system("mv %s %s" % ("mock_boundary.txt", zobovDir)) - os.system("mv %s %s" % ("mock_sphere.txt", zobovDir)) + os.system("mv %s %s" % ("galaxies.txt", outputDir)) + os.system("mv %s %s" % ("mock_galaxies.txt", outputDir)) + os.system("mv %s %s" % ("mock_boundary.txt", outputDir)) + os.system("mv %s %s" % ("mock_sphere.txt", outputDir)) else: # simulation sampleName = sample.fullName @@ -142,16 +153,16 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, if prevSubSample == -1: inputParameterFlag = "" - outputFile = zobovDir+"/zobov_slice_" + sampleName + "_ss" + \ + outputFile = outputDir+"/zobov_slice_" + sampleName + "_ss" + \ thisSubSample keepFraction = float(thisSubSample) subSampleLine = "subsample %g" % keepFraction resubSampleLine = "" firstSubSample = keepFraction else: - inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+\ + inputParameterFlag = "inputParameter " + outputDir+"/zobov_slice_"+\ sampleName+"_ss"+prevSubSample+".par" - outputFile = zobovDir+"/zobov_slice_" + sampleName + "_ss" + \ + outputFile = outputDir+"/zobov_slice_" + sampleName + "_ss" + \ thisSubSample keepFraction = float(thisSubSample)/float(prevSubSample) subSampleLine = "subsample %s" % firstSubSample @@ -243,8 +254,8 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, # remove intermediate files if (prevSubSample != -1): - os.unlink(zobovDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample+".par") - os.unlink(zobovDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample) + os.unlink(outputDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample+".par") + os.unlink(outputDir+"/zobov_slice_"+sampleName+"_ss"+prevSubSample) doneLine = "Done! %5.2e\n" % keepFraction if not jobSuccessful(logFile, doneLine): @@ -256,27 +267,27 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, if jobSuccessful(logFile, doneLine): print("done") # place the final subsample - os.system("mv %s %s" % (zobovDir+"/zobov_slice_"+sampleName+"_ss"+\ - prevSubSample, zobovDir+"/zobov_slice_"+sampleName)) - os.system("mv %s %s" % (zobovDir+"/zobov_slice_"+sampleName+"_ss"+\ - prevSubSample+".par", zobovDir+"/zobov_slice_"+sampleName+".par")) + os.system("mv %s %s" % (outputDir+"/zobov_slice_"+sampleName+"_ss"+\ + prevSubSample, outputDir+"/zobov_slice_"+sampleName)) + os.system("mv %s %s" % (outputDir+"/zobov_slice_"+sampleName+"_ss"+\ + prevSubSample+".par", outputDir+"/zobov_slice_"+sampleName+".par")) if os.access("comoving_distance.txt", os.F_OK): - os.system("mv %s %s" % ("comoving_distance.txt", zobovDir)) + os.system("mv %s %s" % ("comoving_distance.txt", outputDir)) if os.access(parmFile, os.F_OK): os.unlink(parmFile) if os.access("mask_index.txt", os.F_OK): - os.system("mv %s %s" % ("mask_index.txt", zobovDir)) - os.system("mv %s %s" % ("total_particles.txt", zobovDir)) - #os.system("mv %s %s" % ("sample_info.txt", zobovDir)) + os.system("mv %s %s" % ("mask_index.txt", outputDir)) + os.system("mv %s %s" % ("total_particles.txt", outputDir)) + #os.system("mv %s %s" % ("sample_info.txt", outputDir)) # add to sample info file if sample.dataType == "observation": (boxVol, nbar) = getSurveyProps(sample.maskFile, sample.zRange[0], sample.zRange[1], sample.zRange[0], sample.zRange[1], "all", - useComoving=useComoving) + sample.omegaM, useComoving=useComoving) else: iX = float(sample.mySubvolume[0]) iY = float(sample.mySubvolume[1]) @@ -290,16 +301,16 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, boxVol = (xMax-xMin)*(yMax-yMin)*(zMax-zMin) nbar = 1.0 - numTracers = int(open(zobovDir+"/mask_index.txt", "r").read()) - numTotal = int(open(zobovDir+"/total_particles.txt", "r").read()) + numTracers = int(open(outputDir+"/mask_index.txt", "r").read()) + numTotal = int(open(outputDir+"/total_particles.txt", "r").read()) meanSep = (1.*numTracers/boxVol/nbar)**(-1/3.) # save this sample's information - with open(zobovDir+"/sample_info.dat", mode='wb') as output: + with open(outputDir+"/sample_info.dat", mode='wb') as output: pickle.dump(sample, output, pickle.HIGHEST_PROTOCOL) - fp = open(zobovDir+"/sample_info.txt", 'w') + fp = open(outputDir+"/sample_info.txt", 'w') fp.write("Sample name: %s\n" % sample.fullName) fp.write("Sample nickname: %s\n" % sample.nickName) fp.write("Data type: %s\n" % sample.dataType) @@ -323,13 +334,13 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None, fp.close() # ----------------------------------------------------------------------------- -def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, +def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None, numZobovDivisions=None, numZobovThreads=None, mergingThreshold=0.2): sampleName = sample.fullName - datafile = zobovDir+"zobov_slice_"+sampleName + datafile = outputDir+"zobov_slice_"+sampleName logFile = logDir+"/zobov_"+sampleName+".out" @@ -339,8 +350,8 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, os.unlink(vozScript) if sample.dataType == "observation": - maskIndex = open(zobovDir+"/mask_index.txt", "r").read() - totalPart = open(zobovDir+"/total_particles.txt", "r").read() + maskIndex = open(outputDir+"/mask_index.txt", "r").read() + totalPart = open(outputDir+"/total_particles.txt", "r").read() maxDen = mergingThreshold*float(maskIndex)/float(totalPart) else: maskIndex = -1 @@ -349,18 +360,18 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, print(" WARNING! You are using a single ZOBOV division with a simulation. Periodic boundaries will not be respected!") if not (continueRun and jobSuccessful(logFile, "Done!\n")): - for fileName in glob.glob(zobovDir+"/part._"+sampleName+".*"): + for fileName in glob.glob(outputDir+"/part._"+sampleName+".*"): os.unlink(fileName) - if os.access(zobovDir+"/adj_"+sampleName+".dat", os.F_OK): - os.unlink(zobovDir+"adj_"+sampleName+".dat") + if os.access(outputDir+"/adj_"+sampleName+".dat", os.F_OK): + os.unlink(outputDir+"adj_"+sampleName+".dat") - if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK): - os.unlink(zobovDir+"/voidDesc_"+sampleName+".out") + if os.access(outputDir+"/voidDesc_"+sampleName+".out", os.F_OK): + os.unlink(outputDir+"/voidDesc_"+sampleName+".out") cmd = [binPath+"/vozinit", datafile, "0.1", "1.0", str(numZobovDivisions), \ "_"+sampleName, str(numZobovThreads), \ - binPath, zobovDir, str(maskIndex)] + binPath, outputDir, str(maskIndex)] log = open(logFile, 'w') subprocess.call(cmd, stdout=log, stderr=log) log.close() @@ -375,13 +386,13 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, sample.selFunFile != None: # load volumes - volFile = zobovDir+"/vol_"+sampleName+".dat" + volFile = outputDir+"/vol_"+sampleName+".dat" with open(volFile, mode="rb") as File: numPartTot = np.fromfile(File, dtype=np.int32,count=1) vols = np.fromfile(File, dtype=np.float32,count=numPartTot) # load redshifts - partFile = zobovDir+"/zobov_slice_"+sample.fullName + partFile = outputDir+"/zobov_slice_"+sample.fullName with open(partFile, mode="rb") as File: chk = np.fromfile(File, dtype=np.int32,count=1) Np = np.fromfile(File, dtype=np.int32,count=1) @@ -412,7 +423,7 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, redshifts = np.fromfile(File, dtype=np.float32,count=Np) chk = np.fromfile(File, dtype=np.int32,count=1) - # z + # redshift chk = np.fromfile(File, dtype=np.int32,count=1) redshifts = np.fromfile(File, dtype=np.float32,count=Np) chk = np.fromfile(File, dtype=np.int32,count=1) @@ -428,29 +439,29 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, for i in range(len(vols)): vols[i] *= selfunc(redshifts[i]) - volFile = zobovDir+"/vol_weighted_"+sampleName+".dat" + volFile = outputDir+"/vol_weighted_"+sampleName+".dat" with open(volFile, mode='wb') as File: numPartTot.astype(np.int32).tofile(File) vols.astype(np.float32).tofile(File) - volFileToUse = zobovDir+"/vol_weighted_"+sampleName+".dat" + volFileToUse = outputDir+"/vol_weighted_"+sampleName+".dat" else: - volFileToUse = zobovDir+"/vol_"+sampleName+".dat" + volFileToUse = outputDir+"/vol_"+sampleName+".dat" cmd = [binPath+"/jozov2", \ - zobovDir+"/adj_"+sampleName+".dat", \ + outputDir+"/adj_"+sampleName+".dat", \ volFileToUse, \ - zobovDir+"/voidPart_"+sampleName+".dat", \ - zobovDir+"/voidZone_"+sampleName+".dat", \ - zobovDir+"/voidDesc_"+sampleName+".out", \ + outputDir+"/voidPart_"+sampleName+".dat", \ + outputDir+"/voidZone_"+sampleName+".dat", \ + outputDir+"/voidDesc_"+sampleName+".out", \ str(maxDen), str(maskIndex)] log = open(logFile, 'a') subprocess.call(cmd, stdout=log, stderr=log) log.close() # don't need the subbox files - for fileName in glob.glob(zobovDir+"/part._"+sampleName+".*"): + for fileName in glob.glob(outputDir+"/part._"+sampleName+".*"): os.unlink(fileName) if jobSuccessful(logFile, "Done!\n"): @@ -467,19 +478,19 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, # ----------------------------------------------------------------------------- def launchPrune(sample, binPath, - summaryFile=None, logFile=None, zobovDir=None, + summaryFile=None, logFile=None, outputDir=None, continueRun=None, useComoving=False, mergingThreshold=0.2, boundaryTolerance=1.0): sampleName = sample.fullName numVoids = sum(1 for line in \ - open(zobovDir+"/voidDesc_"+sampleName+".out")) + open(outputDir+"/voidDesc_"+sampleName+".out")) numVoids -= 2 if sample.dataType == "observation": - mockIndex = open(zobovDir+"/mask_index.txt", "r").read() - totalPart = open(zobovDir+"/total_particles.txt", "r").read() + mockIndex = open(outputDir+"/mask_index.txt", "r").read() + totalPart = open(outputDir+"/total_particles.txt", "r").read() maxDen = mergingThreshold*float(mockIndex)/float(totalPart) observationLine = " --isObservation" #periodicLine = " --periodic=''" @@ -497,7 +508,7 @@ def launchPrune(sample, binPath, if sample.minVoidRadius == -1: minRadius = -1 - for line in open(zobovDir+"/sample_info.txt"): + for line in open(outputDir+"/sample_info.txt"): if "Estimated mean tracer separation" in line: minRadius = float(line.split()[5]) break @@ -510,13 +521,13 @@ def launchPrune(sample, binPath, if not (continueRun and (jobSuccessful(logFile, "NetCDF: Not a valid ID\n") \ or jobSuccessful(logFile, "Done!\n"))): cmd = binPath - cmd += " --partFile=" + zobovDir+"/zobov_slice_"+str(sampleName) - cmd += " --voidDesc=" + zobovDir+"/voidDesc_"+str(sampleName)+".out" - cmd += " --void2Zone="+zobovDir+"/voidZone_"+str(sampleName)+".dat" - cmd += " --zone2Part=" + zobovDir+"/voidPart_"+str(sampleName)+".dat" - cmd += " --partVol=" + zobovDir+"/vol_"+str(sampleName)+".dat" - cmd += " --partAdj=" + zobovDir+"/adj_"+str(sampleName)+".dat" - cmd += " --extraInfo=" + zobovDir+"/zobov_slice_"+str(sampleName)+\ + cmd += " --partFile=" + outputDir+"/zobov_slice_"+str(sampleName) + cmd += " --voidDesc=" + outputDir+"/voidDesc_"+str(sampleName)+".out" + cmd += " --void2Zone="+outputDir+"/voidZone_"+str(sampleName)+".dat" + cmd += " --zone2Part=" + outputDir+"/voidPart_"+str(sampleName)+".dat" + cmd += " --partVol=" + outputDir+"/vol_"+str(sampleName)+".dat" + cmd += " --partAdj=" + outputDir+"/adj_"+str(sampleName)+".dat" + cmd += " --extraInfo=" + outputDir+"/zobov_slice_"+str(sampleName)+\ ".par" cmd += " --tolerance=" + str(boundaryTolerance) cmd += " --mockIndex=" + str(mockIndex) @@ -529,7 +540,7 @@ def launchPrune(sample, binPath, cmd += periodicLine cmd += useComovingFlag cmd += " --omegaM=" + str(sample.omegaM) - cmd += " --outputDir=" + zobovDir + cmd += " --outputDir=" + outputDir cmd += " --sampleName=" + str(sampleName) log = open(logFile, 'w') #log.write(f"Command is {cmd}\n") @@ -635,7 +646,7 @@ def launchVelocityStack(sample, stack, binPath, velField_file, thisDataPortion=None, logDir=None, voidDir=None, runSuffix=None, - zobovDir=None, + outputDir=None, summaryFile=None, continueRun=None, dataType=None, prefixRun=""): diff --git a/python_source/backend/surveyTools.py b/python_source/backend/surveyTools.py index 90a5ab0..51b6099 100644 --- a/python_source/backend/surveyTools.py +++ b/python_source/backend/surveyTools.py @@ -29,24 +29,20 @@ __all__=['getSurveyProps', 'getNside', 'figureOutMask', 'findEdgeGalaxies'] # ----------------------------------------------------------------------------- # returns the volume and galaxy density for a given redshit slice -def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectionFuncFile=None, useComoving=False): +def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, + omegaM, selectionFuncFile=None, useComoving=False): - LIGHT_SPEED = 299792.458 + #LIGHT_SPEED = 299792.458 mask = healpy.read_map(maskFile) area = (1.*np.size(np.where(mask > 0)) / np.size(mask)) * 4.*np.pi if useComoving: - zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=0.27) - zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=0.27) - selFunMin = LIGHT_SPEED/100.*angularDiameter(selFunMin, Om=0.27) - selFunMax = LIGHT_SPEED/100.*angularDiameter(selFunMax, Om=0.27) + zmin = LIGHT_SPEED/100.*angularDiameter(zmin, Om=omegaM) + zmax = LIGHT_SPEED/100.*angularDiameter(zmax, Om=omegaM) else: - zmin = zmin * 3000 - zmax = zmax * 3000 - selFunMin *= 3000 - selFunMax *= 3000 - volume = area * (zmax**3 - zmin**3) / 3 + zmin = zmin * LIGHT_SPEED/100. + zmax = zmax * LIGHT_SPEED/100. if selectionFuncFile == None: nbar = 1.0 @@ -90,40 +86,91 @@ def getSurveyProps(maskFile, zmin, zmax, selFunMin, selFunMax, portion, selectio # ----------------------------------------------------------------------------- # returns the nside resolution from the given maskfile def getNside(maskFile): - nside = 1.0 - return nside + mask = healpy.read_map(maskFile) + return healpy.get_nside(mask) # ----------------------------------------------------------------------------- -# computes the mask from a given datafile and writes it to a file +# 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. + + return phi, theta + +# ----------------------------------------------------------------------------- +# computes the mask from a given galaxy datafile and writes it to a file def figureOutMask(galFile, nside, outMaskFile): npix = healpy.nside2npix(nside) mask = np.zeros((npix)) -for line in open(galFile): - line = line.split() - RA = np.float(line[3]) - Dec = np.float(line[4]) - z = np.float(line[5]) + for line in open(galFile): + line = line.split() + RA = np.float(line[3]) + Dec = np.float(line[4]) + z = np.float(line[5]) - phi = np.pi/180.*RA - theta = Dec*np.pi/180. - theta = np.pi/2. - Dec*np.pi/180. - pos = np.zeros((3)) + phi, theta = convertAngle(RA, Dec) - pix = healpy.ang2pix(nside, theta, phi) - mask[pix] = 1. + pix = healpy.ang2pix(nside, theta, phi) + mask[pix] = 1. - healpy.write_map(outMaskFile, mask) + healpy.write_map(outMaskFile, mask) return mask # ----------------------------------------------------------------------------- # figures out which galaxies live on a mask edge, and also writes the edge # map to an auxillary file -def findEdgeGalaxies(galFile, maskFile): +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) + else: + zmin *= LIGHT_SPEED/100. + zmax *= LIGHT_SPEED/100. + + + mask = healpy.read_map(maskFile) + nside = healpy.get_nside(mask) + npix = healpy.nside2npix(nside) + edgeMask = np.zeros((npix)) + + edgeFile = open(edgeGalFile, "w") + + for line in open(galFile): + line = line.split() + RA = np.float(line[3]) + Dec = np.float(line[4]) + z = np.float(line[5]) + + if useComoving: + z = LIGHT_SPEED/100.*angularDiameter(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): + + # check the redshift boundaries + + if isOnEdge: + edgeFile.write("1") + else: + edgeFile.write("0") + + edgeFile.close() + healpy.write_map(edgeMaskFile, edgeMask) return