# HG changeset patch # User Sam Skillman # Date 1247101228 25200 # Branch trunk # Node ID 356c0e2a81b234b2d99794f680e7380b0e1987ad # Parent 7caae7bbcb13c6b8485c0db265f7f69dbe1db4e3 [svn r1369] Added support for photon_field=True/False(default). If true, lightcone frb's are reduced by 4*pi*dL^2 where dL is the luminosity distance determined through lagos.Cosmology. Also added support for non-constant dtDataDumps in case any other poor soul decided it was a good idea to get better time resolution at the end of the simulation. Might not be used much, but it doesn't hurt. diff -r 7caae7bbcb13c6b8485c0db265f7f69dbe1db4e3 -r 356c0e2a81b234b2d99794f680e7380b0e1987ad yt/extensions/lightcone/LightCone.py --- a/yt/extensions/lightcone/LightCone.py Wed Jul 08 16:50:27 2009 -0700 +++ b/yt/extensions/lightcone/LightCone.py Wed Jul 08 18:00:28 2009 -0700 @@ -33,6 +33,7 @@ import os import numpy as na import random as rand +import yt.lagos.Cosmology as cosmo class LightCone(object): def __init__(self,EnzoParameterFile,LightConeParameterFile,verbose=True): @@ -79,8 +80,10 @@ # Calculate redshifts for dt data dumps. if (self.enzoParameters.has_key('dtDataDump')): - self._CalculateTimeDumpRedshifts() - + if self.lightConeParameters['dtSplit'] is None: + self._CalculateTimeDumpRedshifts() + else: + self._CalculateSplitTimeDumpRedshifts(self.lightConeParameters['dtSplit'],self.lightConeParameters['dt_new']) # Combine all data dumps. self._CombineDataOutputs() @@ -250,7 +253,7 @@ del haloMaskCube def ProjectLightCone(self,field,weight_field=None,apply_halo_mask=False,node=None, - save_stack=True,save_slice_images=False,flatten_stack=False,**kwargs): + save_stack=True,save_slice_images=False,flatten_stack=False,photon_field=False,**kwargs): "Create projections for light cone, then add them together." # Clear projection stack. @@ -273,6 +276,19 @@ frb = LightConeProjection(output,field,self.pixels,weight_field=weight_field, save_image=save_slice_images, name=name,node=node,**kwargs) + if ytcfg.getint("yt","__parallel_rank") == 0: + if photonfield: + #need to decrement the flux by the luminosity distance. Assume field in frb is in erg/s/cm^2/Hz + co = lagos.Cosmology(HubbleConstantNow = (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']), + OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'], + OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow']) + dL = co.LuminosityDistance(self.lightConeParameters['ObserverRedshift'],output['redshift']) #in Mpc + boxSizeProper = self.enzoParameters['CosmologyComovingBoxSize'] / (self.enzoParameters['CosmologyHubbleConstantNow'] * + (1.0 + output['redshift'])) + pixelarea = (boxSizeProper/self.pixels)**2 #in proper cm^2 + factor = pixelarea/(4.0*na.pi*dL**2) + mylog.info("Distance to slice = %e" % dL) + frb[field] *= factor #in erg/s/cm^2/Hz on observer's image plane. if ytcfg.getint("yt","__parallel_rank") == 0: if (weight_field is not None): @@ -582,6 +598,35 @@ self.timeOutputs.append({'index':index,'redshift':z,'filename':filename}) index += 1 + def _CalculateSplitTimeDumpRedshifts(self,break_index,new_dt): + "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)): + if index <= break_index: + t = ec.InitialTime + (self.enzoParameters['dtDataDump'] * ec.TimeUnits * index) + z = ec.ComputeRedshiftFromTime(t) + else: + t = ec.InitialTime + (self.enzoParameters['dtDataDump'] * ec.TimeUnits * (break_index)) + t += (new_dt * ec.TimeUnits * (index-break_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 @@ -727,7 +772,7 @@ def _SetParameterDefaults(self): "Set some default parameters to avoid problems if they are not in the parameter file." - self.enzoParameters['GlobalDir'] = "" + self.enzoParameters['GlobalDir'] = "./" self.enzoParameters['RedshiftDumpName'] = "RD" self.enzoParameters['RedshiftDumpDir'] = "RD" self.enzoParameters['DataDumpName'] = "DD" @@ -738,6 +783,8 @@ self.lightConeParameters['ObserverRedshift'] = 0.0 self.lightConeParameters['OutputDir'] = "./" self.lightConeParameters['OutputPrefix'] = "LightCone" + self.lightConeParameters['dtSplit'] = None + self.lightConeParameters['dt_new'] = 0.0 EnzoParameterDict = {"CosmologyCurrentRedshift": float, "CosmologyComovingBoxSize": float, @@ -763,4 +810,6 @@ "MinimumCoherentBoxFraction": float, "MinimumSliceDeltaz": float, "OutputDir": str, - "OutputPrefix": str} + "OutputPrefix": str, + "dtSplit": float, + "dt_new": float}