diff --git a/python_tools/void_python_tools/partUtil/partUtil.py b/python_tools/void_python_tools/partUtil/partUtil.py index 46cca76..0f8a1f7 100644 --- a/python_tools/void_python_tools/partUtil/partUtil.py +++ b/python_tools/void_python_tools/partUtil/partUtil.py @@ -156,3 +156,146 @@ def shiftPart(inPart, newCenter, periodicLine, boxLen): np.copysign(boxLen[2], part[shiftUs,2]) return part + +# ----------------------------------------------------------------------------- + +class Bunch: + def __init__(self, **kwds): + self.__dict__.update(kwds) + +class Catalog: + numVoids = 0 + numPartTot = 0 + numZonesTot = 0 + boxLen = np.zeros((3)) + part = None + zones2Parts = None + void2Zones = None + voids = None + +# ----------------------------------------------------------------------------- +def loadVoidCatalog(sampleDir): + #print " Loading particle data..." + sys.stdout.flush() + + catalog = Catalog() + with open(sampleDir+"/sample_info.dat", 'rb') as input: + sample = pickle.load(input) + + print "Loading info..." + infoFile = sampleDir+"/zobov_slice_"+sample.fullName+".par" + File = NetCDFFile(infoFile, 'r') + ranges = np.zeros((3,2)) + ranges[0][0] = getattr(File, 'range_x_min') + ranges[0][1] = getattr(File, 'range_x_max') + ranges[1][0] = getattr(File, 'range_y_min') + ranges[1][1] = getattr(File, 'range_y_max') + ranges[2][0] = getattr(File, 'range_z_min') + ranges[2][1] = getattr(File, 'range_z_max') + catalog.boxLen[0] = ranges[0][1] - ranges[0][0] + catalog.boxLen[1] = ranges[1][1] - ranges[1][0] + catalog.boxLen[2] = ranges[2][1] - ranges[2][0] + File.close() + + print "Loading all particles..." + partData, boxLen, volNorm, isObservationData = loadPart(sampleDir) + numPartTot = len(partData) + catalog.numPartTot = numPartTot + catalog.part = [] + for i in xrange(len(partData)): + catalog.part.append(Bunch(x = partData[i][0], + y = partData[i][1], + z = partData[i][2], + volume = 0, + ra = 0, + dec = 0, + redshift = 0, + uniqueID = 0)) + + + print "Loading volumes..." + volFile = sampleDir+"/vol_"+sample.fullName+".dat" + File = file(volFile) + chk = np.fromfile(File, dtype=np.int32,count=1) + vols = np.fromfile(File, dtype=np.float32,count=numPartTot) + for ivol in xrange(len(vols)): + catalog.part[ivol].volume = vols[ivol] / volNorm + + print "Loading voids..." + fileName = sampleDir+"/voidDesc_central_"+sample.fullName+".out" + catData = np.loadtxt(fileName, comments="#", skiprows=2) + catalog.voids = [] + for line in catData: + catalog.voids.append(Bunch(iVoid = int(line[0]), + voidID = int(line[1]), + coreParticle = line[2], + coreDens = line[3], + zoneVol = line[4], + zoneNumPart = line[5], + numZones = int(line[6]), + voidVol = line[7], + numPart = int(line[8]), + densCon = line[9], + voidProp = line[10], + radius = pow(line[7]/volNorm*3./4./np.pi, 1./3.), + barycenter = np.zeros((3)))) + + print "Read %d voids" % len(catalog.voids) + + print "Loading barycenters..." + iLine = 0 + for line in open(sampleDir+"/barycenters_central_"+sample.fullName+".out"): + line = line.split() + catalog.voids[iLine].barycenter[0] = float(line[1]) + catalog.voids[iLine].barycenter[1] = float(line[2]) + catalog.voids[iLine].barycenter[2] = float(line[3]) + iLine += 1 + + + print "Loading zone-void membership info..." + zoneFile = sampleDir+"/voidZone_"+sample.fullName+".dat" + catalog.void2Zones = [] + File = file(zoneFile) + numZonesTot = np.fromfile(File, dtype=np.int32,count=1) + catalog.numZonesTot = numZonesTot + for iZ in xrange(numZonesTot): + numZones = np.fromfile(File, dtype=np.int32,count=1) + catalog.void2Zones.append(Bunch(numZones = numZones, + zoneIDs = [])) + + for p in xrange(numZones): + zoneID = np.fromfile(File, dtype=np.int32,count=1) + catalog.void2Zones[iZ].zoneIDs.append(zoneID) + + + print "Loading particle-zone membership info..." + zonePartFile = sampleDir+"/voidPart_"+sample.fullName+".dat" + catalog.zones2Parts = [] + File = file(zonePartFile) + chk = np.fromfile(File, dtype=np.int32,count=1) + numZonesTot = np.fromfile(File, dtype=np.int32,count=1) + for iZ in xrange(numZonesTot): + numPart = np.fromfile(File, dtype=np.int32,count=1) + catalog.zones2Parts.append(Bunch(numPart = numPart, + partIDs = [])) + + for p in xrange(numPart): + partID = np.fromfile(File, dtype=np.int32,count=1) + catalog.zones2Parts[iZ].partIDs.append(partID) + + return catalog + + + +# ----------------------------------------------------------------------------- +def getVoidPart(catalog, voidID): + + partOut = [] + for iZ in xrange(catalog.void2Zones[voidID].numZones): + zoneID = catalog.void2Zones[voidID].zoneIDs[iZ] + for p in xrange(catalog.zones2Parts[zoneID].numPart): + partID = catalog.zones2Parts[zoneID].partIDs[p] + partOut.append(catalog.part[partID]) + + return partOut +