[Yt-svn] yt-commit r833 - in trunk/yt/lagos: . lightcone
britton at wrangler.dreamhost.com
britton at wrangler.dreamhost.com
Fri Oct 24 19:08:18 PDT 2008
Author: britton
Date: Fri Oct 24 19:08:17 2008
New Revision: 833
URL: http://yt.spacepope.org/changeset/833
Log:
Adding light cone generator to yt. A few updates coming soon, so
hang tight for a few days.
24 October, 2008
Britton Smith
Added:
trunk/yt/lagos/lightcone/
trunk/yt/lagos/lightcone/LightCone.py
trunk/yt/lagos/lightcone/LightConeProjection.py
trunk/yt/lagos/lightcone/__init__.py
Modified:
trunk/yt/lagos/__init__.py
trunk/yt/lagos/setup.py
Modified: trunk/yt/lagos/__init__.py
==============================================================================
--- trunk/yt/lagos/__init__.py (original)
+++ trunk/yt/lagos/__init__.py Fri Oct 24 19:08:17 2008
@@ -56,6 +56,9 @@
import time
+from Cosmology import *
+from EnzoCosmology import *
+
if ytcfg.getboolean("lagos","useswig"):
try:
from yt.enki import EnzoInterface
Added: trunk/yt/lagos/lightcone/LightCone.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/lightcone/LightCone.py Fri Oct 24 19:08:17 2008
@@ -0,0 +1,388 @@
+from yt.lagos.lightcone import *
+import numpy as na
+import random as rand
+import tables as h5
+
+class LightCone(object):
+ def __init__(self,EnzoParameterFile,LightConeParameterFile):
+ self.EnzoParameterFile = EnzoParameterFile
+ self.LightConeParameterFile = LightConeParameterFile
+ self.enzoParameters = {}
+ self.lightConeParameters = {}
+ self.redshiftOutputs = []
+ self.timeOutputs = []
+ self.allOutputs = []
+ self.lightConeSolution = []
+ self.projectionStack = []
+ self.projectionWeightFieldStack = []
+
+ # Set some parameter defaults.
+ self._SetParameterDefaults()
+
+ # Read parameters.
+ self._ReadLightConeParameterFile()
+ self._ReadEnzoParameterFile()
+
+ # Calculate redshifts for dt data dumps.
+ if (self.enzoParameters.has_key('dtDataDump')):
+ self._CalculateTimeDumpRedshifts()
+
+ # Combine all data dumps.
+ self._CombineDataOutputs()
+
+ # Calculate maximum delta z for each data dump.
+ self._CalculateDeltaZMax()
+
+ def CalculateLightConeSolution(self):
+ "Create list of projections to be added together to make the light cone."
+
+ # Make light cone using minimum number of projections.
+ if (self.lightConeParameters['UseMinimumNumberOfProjections']):
+
+ z_Tolerance = 1e-4
+ z = self.lightConeParameters['InitialRedshift']
+
+ # fill redshift space with projections
+ while ((z > self.lightConeParameters['FinalRedshift']) and
+ (na.fabs(z - self.lightConeParameters['FinalRedshift']) > z_Tolerance)):
+ # Sort data outputs by proximity to current redsfhit.
+ self.allOutputs.sort(key=lambda obj:na.fabs(z - obj['redshift']))
+ # For first data dump, choose closest to desired redshift.
+ if (len(self.lightConeSolution) == 0):
+ self.lightConeSolution.append(self.allOutputs[0])
+ # Start with data dump closest to desired redshift and move backward
+ # until one is within delta z of last output in solution list.
+ else:
+ output = self.allOutputs[0]
+ while (z > output['redshift']):
+ output = output['previous']
+ if (output is None):
+ print "CalculateLightConeSolution: search for data output went off the end the stack."
+ print "Could not calculate light cone solution."
+ return
+ if (output['redshift'] == self.lightConeSolution[-1]['redshift']):
+ print "CalculateLightConeSolution: No data dump between z = %f and %f." % \
+ ((self.lightConeSolution[-1]['redshift'] - self.lightConeSolution[-1]['deltaz']),
+ self.lightConeSolution[-1]['redshift'])
+ print "Could not calculate light cone solution."
+ return
+ self.lightConeSolution.append(output)
+ z = self.lightConeSolution[-1]['redshift'] - self.lightConeSolution[-1]['deltaz']
+
+ # Make light cone using maximum number of projections (minimum spacing).
+ else:
+ deltazMin = 0.01
+
+ # Sort data outputs by proximity to current redsfhit.
+ self.allOutputs.sort(key=lambda obj:na.fabs(self.lightConeParameters['InitialRedshift'] - obj['redshift']))
+ # For first data dump, choose closest to desired redshift.
+ self.lightConeSolution.append(self.allOutputs[0])
+
+ nextOutput = self.lightConeSolution[-1]['next']
+ while (nextOutput is not None):
+ if (nextOutput['redshift'] <= self.lightConeParameters['FinalRedshift']):
+ break
+ if ((self.lightConeSolution[-1]['redshift'] - nextOutput['redshift']) > deltazMin):
+ self.lightConeSolution.append(nextOutput)
+ nextOutput = nextOutput['next']
+
+ print "CalculateLightConeSolution: Used %d data dumps to get from z = %f to %f." % (len(self.lightConeSolution),
+ self.lightConeParameters['InitialRedshift'],
+ self.lightConeParameters['FinalRedshift'])
+
+ # Calculate projection sizes, and get random projection axes and centers.
+ rand.seed(self.lightConeParameters['RandomSeed'])
+ co = lagos.Cosmology(HubbleConstantNow = (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
+ OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
+ OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
+
+ for q in range(len(self.lightConeSolution)):
+ del self.lightConeSolution[q]['previous']
+ del self.lightConeSolution[q]['next']
+ if (q == len(self.lightConeSolution) - 1):
+ z_next = self.lightConeParameters['FinalRedshift']
+ else:
+ z_next = self.lightConeSolution[q+1]['redshift']
+ # Calculate fraction of box required for a depth of delta z
+ self.lightConeSolution[q]['DepthBoxFraction'] = co.ComovingRadialDistance(z_next,self.lightConeSolution[q]['redshift']) * \
+ self.enzoParameters['CosmologyHubbleConstantNow'] / self.enzoParameters['CosmologyComovingBoxSize']
+ # Simple error check to make sure more than 100% of box depth is never required.
+ if (self.lightConeSolution[q]['DepthBoxFraction'] > 1.0):
+ print "Warning: box fraction required to go from z = %f to %f is %f" % (self.lightConeSolution[q]['redshift'],z_next,
+ self.lightConeSolution[q]['DepthBoxFraction'])
+ print "Full box delta z is %f, but it is %f to the next data dump." % (self.lightConeSolution[q]['deltaz'],
+ self.lightConeSolution[q]['redshift']-z_next)
+
+ # Calculate fraction of box required for width corresponding to requested image size.
+ scale = co.AngularScale_1arcsec_kpc(self.lightConeParameters['FinalRedshift'],self.lightConeSolution[q]['redshift'])
+ size = self.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 * scale / 1000.0
+ boxSizeProper = self.enzoParameters['CosmologyComovingBoxSize'] / (self.enzoParameters['CosmologyHubbleConstantNow'] *
+ (1.0 + self.lightConeSolution[q]['redshift']))
+ self.lightConeSolution[q]['WidthBoxFraction'] = size / boxSizeProper
+
+ # Get random projection axis and center.
+ self.lightConeSolution[q]['ProjectionAxis'] = rand.randint(0,2)
+ self.lightConeSolution[q]['ProjectionCenter'] = [rand.random(),rand.random(),rand.random()]
+
+ # Clear out some stuff.
+ del self.allOutputs
+ del self.redshiftOutputs
+ del self.timeOutputs
+ del co
+
+ def ProjectLightCone(self,field,weight_field=None):
+ "Create projections for light cone, then add them together."
+
+ if not(self.lightConeParameters['OutputDir'].endswith("/")):
+ self.lightConeParameters['OutputDir'] += "/"
+
+ pixels = int(self.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 / \
+ self.lightConeParameters['ImageResolutionInArcSeconds'])
+
+ for q,output in enumerate(self.lightConeSolution):
+ name = "%s%s_%04d_%04d" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'],
+ q,len(self.lightConeSolution))
+ output['object'] = EnzoStaticOutput(output['filename'])
+ frb = LightConeProjection(output,field,pixels,weight_field=weight_field,
+ save_image=self.lightConeParameters['SaveLightConeSlices'],name=name)
+ if (weight_field is not None):
+ # Data come back normalized by the weight field.
+ # Undo that so it can be added up for the light cone.
+ self.projectionStack.append(frb[field]*frb['weight_field'])
+ self.projectionWeightFieldStack.append(frb['weight_field'])
+ else:
+ self.projectionStack.append(frb[field])
+
+ # Unless this is the last slice, delete the dataset object.
+ # The last one will be saved to make the plot collection.
+ if (q < len(self.lightConeSolution) - 1):
+ del output['object']
+
+ # Add up slices to make light cone projection.
+ if (weight_field is None):
+ lightConeProjection = sum(self.projectionStack)
+ else:
+ lightConeProjection = sum(self.projectionStack) / sum(self.projectionWeightFieldStack)
+
+ filename = "%s%s" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
+
+ # Save the last fixed resolution buffer for the plot collection,
+ # but replace the data with the full light cone projection data.
+ frb.data[field] = lightConeProjection
+
+ # Make a plot collection for the light cone projection.
+ pc = raven.PlotCollection(self.lightConeSolution[-1]['object'],center=[0.5,0.5,0.5])
+ pc.add_fixed_resolution_plot(frb,field)
+ pc.save(filename)
+
+ # Return the plot collection so the user can remake the plot if they want.
+ return pc
+
+ def SaveLightConeSolution(self,file="light_cone.dat"):
+ "Write out a text file with information on light cone solution."
+
+ print "Saving light cone solution to %s." % file
+
+ f = open(file,'w')
+ f.write("EnzoParameterFile = %s\n" % self.EnzoParameterFile)
+ f.write("LightConeParameterFile = %s\n" % self.LightConeParameterFile)
+ f.write("\n")
+ for parameter in self.lightConeParameters:
+ f.write("%s = %s\n" % (parameter,str(self.lightConeParameters[parameter])))
+ f.write("\n")
+ for q,output in enumerate(self.lightConeSolution):
+ f.write("Proj %04d, %s, z = %f, depth/box = %f, width/box = %f, axis = %d, center = %f, %f, %f\n" %
+ (q,output['filename'],output['redshift'],output['DepthBoxFraction'],output['WidthBoxFraction'],
+ output['ProjectionAxis'],output['ProjectionCenter'][0],output['ProjectionCenter'][1],output['ProjectionCenter'][2]))
+ f.close()
+
+ def SaveLightConeStack(self,filename=None):
+ "Save the light cone projection stack as a 3d array in and hdf5 file."
+ if (filename is None):
+ filename = "%s%s_data" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
+ filename += ".h5"
+
+ if (len(self.projectionStack) == 0):
+ print "SaveLightConeStack: no projection data loaded."
+ return
+
+ print "Writing light cone data to %s." % filename
+
+ output = h5.openFile(filename, "a")
+
+ self.projectionStack = na.array(self.projectionStack)
+ output.createArray("/", "data",self.projectionStack)
+
+ if (len(self.projectionWeightFieldStack) > 0):
+ self.projectionWeightFieldStack = na.array(self.projectionWeightFieldStack)
+ output.createArray("/","weight",self.projectionWeightFieldStack)
+
+ output.close()
+
+ def _CalculateDeltaZMax(self):
+ "Calculate delta z that corresponds to full box length going from z to (z - delta z)."
+ co = lagos.Cosmology(HubbleConstantNow = (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
+ OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
+ OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
+
+ d_Tolerance = 1e-4
+ max_Iterations = 100
+
+ for output in self.allOutputs:
+ z = output['redshift']
+
+ # Calculate delta z that corresponds to the length of the box at a given redshift.
+ # Use Newton's method to calculate solution.
+ z1 = z
+ z2 = z1 - 0.1 # just an initial guess
+ distance1 = 0.0
+ iteration = 1
+
+ # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
+ distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+
+ while ((na.fabs(distance2-self.enzoParameters['CosmologyComovingBoxSize'])/distance2) > d_Tolerance):
+ m = (distance2 - distance1) / (z2 - z1)
+ z1 = z2
+ distance1 = distance2
+ z2 = ((self.enzoParameters['CosmologyComovingBoxSize'] - distance2) / m) + z2
+ distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+ iteration += 1
+ if (iteration > max_Iterations):
+ print "CalculateDeltaZMax: Warning - max iterations exceeded for z = %f (delta z = %f)." % (z,na.fabs(z2-z))
+ break
+ output['deltaz'] = na.fabs(z2-z)
+
+ del co
+
+ def _CalculateTimeDumpRedshifts(self):
+ "Calculates redshift values for the dt data dumps."
+ ec = lagos.EnzoCosmology(HubbleConstantNow = \
+ (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
+ OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
+ OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'],
+ InitialRedshift = self.enzoParameters['CosmologyInitialRedshift'])
+
+ z_Tolerance = 1e-3
+ z = self.enzoParameters['CosmologyInitialRedshift']
+ index = 0
+ while ((z > self.enzoParameters['CosmologyFinalRedshift']) and
+ (na.fabs(z-self.enzoParameters['CosmologyFinalRedshift']) > z_Tolerance)):
+ t = ec.InitialTime + (self.enzoParameters['dtDataDump'] * ec.TimeUnits * index)
+ z = ec.ComputeRedshiftFromTime(t)
+ filename = "%s/%s%04d/%s%04d" % (self.enzoParameters['GlobalDir'],
+ self.enzoParameters['DataDumpDir'],
+ index,
+ self.enzoParameters['DataDumpDir'],
+ index)
+ self.timeOutputs.append({'index':index,'redshift':z,'filename':filename})
+ index += 1
+
+ def _CombineDataOutputs(self):
+ "Combines redshift and time data into one sorted list."
+ self.allOutputs = self.redshiftOutputs + self.timeOutputs
+ self.allOutputs.sort(reverse=True,key=lambda obj:obj['redshift'])
+ for q in range(len(self.allOutputs)):
+ del self.allOutputs[q]['index']
+ if (q == 0):
+ self.allOutputs[q]['previous'] = None
+ self.allOutputs[q]['next'] = self.allOutputs[q+1]
+ elif (q == len(self.allOutputs) - 1):
+ self.allOutputs[q]['previous'] = self.allOutputs[q-1]
+ self.allOutputs[q]['next'] = None
+ else:
+ self.allOutputs[q]['previous'] = self.allOutputs[q-1]
+ self.allOutputs[q]['next'] = self.allOutputs[q+1]
+
+ def _ReadEnzoParameterFile(self):
+ "Reads an Enzo parameter file looking for cosmology and output parameters."
+ lines = open(self.EnzoParameterFile).readlines()
+ for line in lines:
+ if line.find("#") >= 0: # Keep the commented lines
+ line=line[:line.find("#")]
+ line=line.strip()
+ if len(line) < 2:
+ continue
+ try:
+ param, vals = map(str,line.split("="))
+ param = param.strip()
+ vals = vals.strip()
+ except ValueError:
+ print "Skipping line: %s" % line
+ continue
+ if EnzoParameterDict.has_key(param):
+ t = map(EnzoParameterDict[param], vals.split())
+ if len(t) == 1:
+ self.enzoParameters[param] = t[0]
+ else:
+ self.enzoParameters[param] = t
+ elif param.startswith("CosmologyOutputRedshift["):
+ index = param[param.find("[")+1:param.find("]")]
+ self.redshiftOutputs.append({'index':int(index),
+ 'redshift':float(vals)})
+
+ # Add filenames to redshift outputs.
+ for output in self.redshiftOutputs:
+ output["filename"] = "%s/%s%04d/%s%04d" % (self.enzoParameters['GlobalDir'],
+ self.enzoParameters['RedshiftDumpDir'],
+ output['index'],
+ self.enzoParameters['RedshiftDumpDir'],
+ output['index'])
+
+ def _ReadLightConeParameterFile(self):
+ "Reads a light cone parameter file."
+ lines = open(self.LightConeParameterFile).readlines()
+ for line in lines:
+ if line.find("#") >= 0: # Keep the commented lines
+ line=line[:line.find("#")]
+ line=line.strip()
+ if len(line) < 2:
+ continue
+ try:
+ param, vals = map(str,line.split("="))
+ param = param.strip()
+ vals = vals.strip()
+ except ValueError:
+ print "Skipping line: %s" % line
+ continue
+ if LightConeParameterDict.has_key(param):
+ t = map(LightConeParameterDict[param], vals.split())
+ if len(t) == 1:
+ self.lightConeParameters[param] = t[0]
+ else:
+ self.lightConeParameters[param] = t
+
+ def _SetParameterDefaults(self):
+ "Set some default parameters to avoid problems if they are not in the parameter file."
+ self.enzoParameters['GlobalDir'] = ""
+ self.enzoParameters['RedshiftDumpName'] = "RD"
+ self.enzoParameters['RedshiftDumpDir'] = "RD"
+ self.enzoParameters['DataDumpName'] = "DD"
+ self.enzoParameters['DataDumpDir'] = "DD"
+ self.lightConeParameters['UseMinimumNumberOfProjections'] = 1
+ self.lightConeParameters['OutputDir'] = "./"
+ self.lightConeParameters['OutputPrefix'] = "LightCone"
+
+EnzoParameterDict = {"CosmologyCurrentRedshift": float,
+ "CosmologyComovingBoxSize": float,
+ "CosmologyOmegaMatterNow": float,
+ "CosmologyOmegaLambdaNow": float,
+ "CosmologyHubbleConstantNow": float,
+ "CosmologyInitialRedshift": float,
+ "CosmologyFinalRedshift": float,
+ "dtDataDump": float,
+ "RedshiftDumpName": str,
+ "RedshiftDumpDir": str,
+ "DataDumpName": str,
+ "DataDumpDir": str,
+ "GlobalDir" : str}
+
+LightConeParameterDict = {"InitialRedshift": float,
+ "FinalRedshift": float,
+ "FieldOfViewInArcMinutes": float,
+ "ImageResolutionInArcSeconds": float,
+ "RandomSeed": int,
+ "UseMinimumNumberOfProjections": int,
+ "SaveLightConeSlices": int,
+ "OutputDir": str,
+ "OutputPrefix": str}
Added: trunk/yt/lagos/lightcone/LightConeProjection.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/lightcone/LightConeProjection.py Fri Oct 24 19:08:17 2008
@@ -0,0 +1,159 @@
+from yt.lagos.lightcone import *
+import random as rand
+
+def LightConeProjection(lightConeSlice,field,pixels,weight_field=None,save_image=False,name=""):
+ "Create a single projection to be added into the light cone stack."
+
+ # Use some projection parameters to seed random number generator to make unique node name.
+ rand.seed(lightConeSlice['ProjectionCenter'][0] + lightConeSlice['ProjectionCenter'][1] + lightConeSlice['ProjectionCenter'][2] +
+ lightConeSlice['ProjectionAxis'] + lightConeSlice['DepthBoxFraction'] + lightConeSlice['WidthBoxFraction'])
+ node_name = "%s_%s_%09d" % (field,weight_field,rand.randint(1,999999999))
+
+ print "Making projection at z = %f from %s." % (lightConeSlice['redshift'],lightConeSlice['filename'])
+
+ # Make an Enzo data object.
+ dataset_object = lightConeSlice['object']
+
+ # 1. The Depth Problem
+ # Excise a region from the box with the specified depth and make projection with that.
+ region_center = [0.5,0.5,0.5]
+
+ region_left = [0.0,0.0,0.0]
+ region_left[lightConeSlice['ProjectionAxis']] = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] - \
+ 0.5 * lightConeSlice['DepthBoxFraction']
+
+ region_right = [1.0,1.0,1.0]
+ region_right[lightConeSlice['ProjectionAxis']] = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] + \
+ 0.5 * lightConeSlice['DepthBoxFraction']
+
+ periodic_region = dataset_object.h.periodic_region(region_center,region_left,region_right)
+
+ # Make plot collection.
+ pc = raven.PlotCollection(dataset_object,center=region_center)
+
+ # Make projection.
+ pc.add_projection(field,lightConeSlice['ProjectionAxis'],weight_field=weight_field,source=periodic_region,use_colorbar=True,
+ node_name=node_name)
+
+ # Serialize projection data.
+ pc.plots[0].data._serialize(node_name,force=True)
+
+ # Delete region, maybe save some ram.
+ del periodic_region
+
+ # 2. The Tile Problem
+ # Tile projection to specified width.
+
+ # Original projection data.
+ original_px = pc.plots[0].data['px']
+ original_py = pc.plots[0].data['py']
+ original_pdx = pc.plots[0].data['pdx']
+ original_pdy = pc.plots[0].data['pdy']
+ original_field = pc.plots[0].data[field]
+ original_weight_field = pc.plots[0].data['weight_field']
+
+ # Copy original into offset positions to make tiles.
+ for x in range(int(na.ceil(lightConeSlice['WidthBoxFraction']))):
+ for y in range(int(na.ceil(lightConeSlice['WidthBoxFraction']))):
+ if ((x + y) > 0):
+ pc.plots[0].data['px'] = na.concatenate([pc.plots[0].data['px'],original_px+x])
+ pc.plots[0].data['py'] = na.concatenate([pc.plots[0].data['py'],original_py+y])
+ pc.plots[0].data['pdx'] = na.concatenate([pc.plots[0].data['pdx'],original_pdx])
+ pc.plots[0].data['pdy'] = na.concatenate([pc.plots[0].data['pdy'],original_pdy])
+ pc.plots[0].data[field] = na.concatenate([pc.plots[0].data[field],original_field])
+ pc.plots[0].data['weight_field'] = na.concatenate([pc.plots[0].data['weight_field'],original_weight_field])
+
+ # Delete originals.
+ del original_px
+ del original_py
+ del original_pdx
+ del original_pdy
+ del original_field
+ del original_weight_field
+
+ # 3. The Shift Problem
+ # Shift projection by random x and y offsets.
+
+ offset = lightConeSlice['ProjectionCenter']
+ # Delete depth coordinate.
+ del offset[lightConeSlice['ProjectionAxis']]
+
+ # Shift x and y positions.
+ pc.plots[0]['px'] -= offset[0]
+ pc.plots[0]['py'] -= offset[1]
+
+ # Wrap off-edge cells back around to other side (periodic boundary conditions).
+ pc.plots[0]['px'][pc.plots[0]['px'] < 0] += na.ceil(lightConeSlice['WidthBoxFraction'])
+ pc.plots[0]['py'][pc.plots[0]['py'] < 0] += na.ceil(lightConeSlice['WidthBoxFraction'])
+
+ # After shifting, some cells have fractional coverage on both sides of the box.
+ # Find those cells and make copies to be placed on the other side.
+
+ # Cells hanging off the right edge.
+ add_x_px = pc.plots[0]['px'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_x_px -= na.ceil(lightConeSlice['WidthBoxFraction'])
+ add_x_py = pc.plots[0]['py'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_x_pdx = pc.plots[0]['pdx'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_x_pdy = pc.plots[0]['pdy'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_x_field = pc.plots[0][field][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_x_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+
+ # Cells hanging off the left edge.
+ add2_x_px = pc.plots[0]['px'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+ add2_x_px += na.ceil(lightConeSlice['WidthBoxFraction'])
+ add2_x_py = pc.plots[0]['py'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+ add2_x_pdx = pc.plots[0]['pdx'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+ add2_x_pdy = pc.plots[0]['pdy'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+ add2_x_field = pc.plots[0][field][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+ add2_x_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+
+ # Cells hanging off the top edge.
+ add_y_px = pc.plots[0]['px'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_y_py = pc.plots[0]['py'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_y_py -= na.ceil(lightConeSlice['WidthBoxFraction'])
+ add_y_pdx = pc.plots[0]['pdx'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_y_pdy = pc.plots[0]['pdy'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_y_field = pc.plots[0][field][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+ add_y_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+
+ # Cells hanging off the bottom edge.
+ add2_y_px = pc.plots[0]['px'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+ add2_y_py = pc.plots[0]['py'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+ add2_y_py += na.ceil(lightConeSlice['WidthBoxFraction'])
+ add2_y_pdx = pc.plots[0]['pdx'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+ add2_y_pdy = pc.plots[0]['pdy'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+ add2_y_field = pc.plots[0][field][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+ add2_y_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+
+ # Add the hanging cells back to the projection data.
+ pc.plots[0].data['px'] = na.concatenate([pc.plots[0]['px'],add_x_px,add_y_px,add2_x_px,add2_y_px])
+ pc.plots[0].data['py'] = na.concatenate([pc.plots[0]['py'],add_x_py,add_y_py,add2_x_py,add2_y_py])
+ pc.plots[0].data['pdx'] = na.concatenate([pc.plots[0]['pdx'],add_x_pdx,add_y_pdx,add2_x_pdx,add2_y_pdx])
+ pc.plots[0].data['pdy'] = na.concatenate([pc.plots[0]['pdy'],add_x_pdy,add_y_pdy,add2_x_pdy,add2_y_pdy])
+ pc.plots[0].data[field] = na.concatenate([pc.plots[0][field],add_x_field,add_y_field,add2_x_field,add2_y_field])
+ pc.plots[0].data['weight_field'] = na.concatenate([pc.plots[0]['weight_field'],
+ add_x_weight_field,add_y_weight_field,add2_x_weight_field,add2_y_weight_field])
+
+ # Delete original copies of hanging cells.
+ del add_x_px,add_y_px,add2_x_px,add2_y_px
+ del add_x_py,add_y_py,add2_x_py,add2_y_py
+ del add_x_pdx,add_y_pdx,add2_x_pdx,add2_y_pdx
+ del add_x_pdy,add_y_pdy,add2_x_pdy,add2_y_pdy
+ del add_x_field,add_y_field,add2_x_field,add2_y_field
+ del add_x_weight_field,add_y_weight_field,add2_x_weight_field,add2_y_weight_field
+
+ # Tiles were made rounding up the width to the nearest integer.
+ # Cut off the edges to get the specified width.
+ pc.set_xlim(0,lightConeSlice['WidthBoxFraction'])
+ pc.set_ylim(0,lightConeSlice['WidthBoxFraction'])
+
+ # Save an image if requested.
+ if (save_image):
+ pc.save(name)
+
+ # Create fixed resolution buffer to return back to the light cone object.
+ # These buffers will be stacked together to make the light cone.
+ frb = raven.FixedResolutionBuffer(pc.plots[0].data,(0,lightConeSlice['WidthBoxFraction'],0,lightConeSlice['WidthBoxFraction']),
+ (pixels,pixels),antialias=True)
+
+ return frb
Added: trunk/yt/lagos/lightcone/__init__.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/lightcone/__init__.py Fri Oct 24 19:08:17 2008
@@ -0,0 +1,6 @@
+from yt.mods import *
+
+from LightConeProjection import *
+from LightCone import *
+
+
Modified: trunk/yt/lagos/setup.py
==============================================================================
--- trunk/yt/lagos/setup.py (original)
+++ trunk/yt/lagos/setup.py Fri Oct 24 19:08:17 2008
@@ -23,6 +23,7 @@
config.add_extension("Interpolators", "yt/lagos/Interpolators.c")
#config.add_extension("DepthFirstOctree", "yt/lagos/DepthFirstOctree.c")
config.add_subpackage("hop")
+ config.add_subpackage("lightcone")
H5dir = check_for_hdf5()
if H5dir is not None:
include_dirs=[os.path.join(H5dir,"include")]
More information about the yt-svn
mailing list