Hi all,
I recently attempted to reproduce fig 3 of
Hallman et al. ApJ 671:27 (2007) using HaloProfiler. I've attached a
png of the plot for those of you unfamiliar with it. It was made on
three of the L7 datasteps, z=0.1, 1.0 and 2.0.
After some trials
and tribulations (of my own creation), I was able to pretty closely
reproduce the z=0.1 results using yt. However, the z=1.0 and 2.0 lines
are quite a bit different. I have a very ugly plot which I could share,
but the punchline is my two lines are quite a bit lower in overall
count to the Hallman fig. I am trying to figure out the source of this
difference.
I have access to Brian's L7 directories on gpfs-wan
where he ran HOP for these datasets. With yt-hop I can create the exact
same halo lists. Well, actually, there is a small deviation in particle
counts for a small subset of the haloes, but they're all well less than
1% of the total number of particles in that halo. Of course, the halo
centers and radii are what's important, and the centers are nearly
identical to Brian's data.
I also found the files that list the
virial quantities for these datasets that made the Hallman figure, and
I confirmed it by over-plotting the data on top of the Hallman figure.
I can't quote specifics because Kraken's scratch is stuck right now,
but comparing the old data to the HaloProfiler data for z=2.0 I saw
that the old data had nearly twice as many virial mass values, and
those virial masses were roughly twice what HaloProfiler gives. The old
virial radius was also roughly twice the HaloProfiler value.
In
order to try to understand this differences, I've inspected the code
that calculated the old virial stats. I've put a quick summary below of
how it was done. It appears to be different than how yt does it, but I
still haven't figured out where a factor of roughly two can come from.
There are those on this list who can better summarize how yt does it.
I
know I haven't fully investigated this, but I'd like to start the
discussion anyway. I need to ask Rick what he used for the input halo
radii below. I ran a quick test using HaloProfiler with double the
input radius for each halo and it didn't change the virial radius or
number of haloes very much, so for HaloProfiler (at least) it's not a
very sensitive variable in that direction.
# read in the enzo data grids, for an input halo define a sphere
-for each grid
-quick check to see if the left and right edges of grid are inside the sphere
-compute some physical constants (like densityconversion)
-for each grid cell compute the distance from the *center* of the cell to the center of the sphere, with periodic boundary conditions considered
-if that distance is less than the full sphere
-for n in bins (equally logarithmically spaced)
-if that distance is less than the radius of the bin
-gas_mass=cell_dens * cellvolume * densityconversion
-profile_gas[n] += gas_mass
-break
-for each particle
-compute distance from each particle to the center of the sphere
-if that distance is less than the size of the sphere
-for n in bins
-if that distance is less than the radius of the bin
-part_mass = particlemass * cellvolume * densityconversion
-profile_dm[n] += part_mass
-break
# now calculate the virial stuff, first some constants
Esquared = OmegaMatterNow * POW(1.0+CurrentRedshift, 3) +
OmegaCurvatureNow * POW(1.0+CurrentRedshift, 2) +
OmegaLambdaNow;
critical_dens = 2.78e11*POW(HubbleConstantNow, 2) * Esquared;
-for n in bins
-annulus[n] = profile_gas[n]+profile_dm[n]
-if n>0
-total[n] += annulus[n-1]
-overdensity[n] = total[n] / (POW((radius of bin n+1)*BoxSize, 3) * 4.0*pi/3.0) / critical_dens
-for n in bins
-if overdensity[n] <= virial_overdensity
-rvir = (radius of bin n-1) + (size of bin n) * (virial_overdensity -
overdensity[n-1]) / (overdensity[n] - overdensity[n-1] + tiny_number)
-break
mvir = POW(rvir*BoxSize, 3)*4.0/3.0*pi*critical_dens*virial_overdensity
_______________________________________________________
sskory(a)physics.ucsd.edu o__ Stephen Skory
http://physics.ucsd.edu/~sskory/ _.>/ _Graduate Student
________________________________(_)_\(_)_______________
