[yt-users] halos

[gso at physics.ucsd.edu]

Hi, Matt and Stephen

Yeah our little discussion got side tracked after other stuff came along. 
I've put in some more thought on this over the past 2 weeks, and I think
I'm ready to start working on this.

"First let me correct one of your misconceptions. There is no halo data
>> stored in the .yt file."

Yeah I know, I was saying "If" we have access to it in the original email:-)

"you need to decide
>> how you want to deal with particles. Halos are generally not
>> ellipsoidal."

I think it will be sufficient for the analysis I want to do, which is
comparing the stars and dm particle mass to the baryon mass using the
radius coming from the halo profiler and center of mass and draw a sphere
(I think the way it was ACTUALLY done in the past was to calculating the
quantities inside a box centered at center of mass with sides 2*r instead
of a sphere, so it's even worse, an ellipsoid would be an improvement
already).  I'm open to suggestions on drawing other container types around
haloes, but I've only been thinking of ellipsoid so far.

I've been mostly a YT user so I haven't thought of doing things in YT for
the benefit for others (kindda selfish I know...), but I think for this
analysis it's much simpler once there's a defined data container, and this
container may be useful in other areas as well.  This is my first attempt
to answering these questions, if someone can point me to how they're done
for the other types of data containers would probably enhance my response
in more of a YT framework.

1) How to tell if a grid intersects an ellipse
2) How to tell if a grid is completely enclosed in an ellipse
3) How to tell if a given x,y,z value is enclosed in an ellipse

Right now I'm envisioning solving them this way...(WARNING: a bit long)

FIRST we need the equation of the ellipsoid:
[(x-x0)/a]^2+[(y-y0)/b]^2+[(z-z0)/c]^2=1

step 1:
find the longest axis-- I envision two ways
1) Find the distance between the two furthest points and draw a line and
call it 2*a ~O(N^2) calculation where N is the number of particles.
2) Find the furthest distance of a particle away from the center of mass
and draw a line and call it 1*a ~O(N), center of mass already found by
halo profiler.

of course method 2 would be less accurate and potentially include a lot of
area empty of particles on the opposite side of the furthest particle, but
it's still way better than drawing a box.

step 2:
find the point furthest away from the longest axis call it p2
with the axis of longest distance a or 2*a defined as center of mass to a
pt or distance between two points
(skipping math that's ~O(N))
with this point's (x,y) solve the equation of an ellipse (2D)
[(x-x0)/a]^2+[(y-y0)/b]^2 = 1 for b
(this is non-trivial as Matt and I discussed, requires rotation of
coordinate axes so that the pts lie in the same plane)

step 3:
find the pt furthest away from the plane as defined by
1) containing the line of distance 2*a and the point p2
2) containing the line from center of mass to furthest point distance a
away and the point p2
(both operations should be only ~O(N), with again rotations)

Once we found the point p3 furthest from that plane
using p3's (x,y,z) solve the eqn [(x-x0)/a]^2+[(y-y0)/b]^2+[(z-z0)/c]^2=1
for c

SECOND we check for cell centers that are within this ellipsoid.
[(x-x0)/a]^2+[(y-y0)/b]^2+[(z-z0)/c]^2  <= 1
This is certainly a short cut and gross exaggeration for small ellipsoids,
so maybe set a threshold like the how the halo profiler does not include
haloes of less than certain pts, we don't consider ellipsoids of less than
some number of cells.

we don't need to check every cells, just the ones within a radius of "a"
around the center located at
1) mid pt between the line length 2*a or
2) center of mass of halo

Instead of rotating and translating every cell center (x,y,z) to the
coordinate of the ellipsoid, we can probably rotate and translate the
equation of the ellipse that we're checking each cell center by to save
computation complexity.

If any of this seems illogical in terms of YT or long winded, I apologize
in advance, I don't know too much about the inner workings of YT.  I can
move the discussion of the math I omitted to the developer's mailing list
or continue offline with Matt and Stephen whichever one you guys prefer
:-)

From
G.S.

> Hi Stephen and Geoffrey,
>
> On Fri, Mar 4, 2011 at 4:12 PM, Stephen Skory <stephenskory at yahoo.com>
> wrote:
>> Geoffrey,
>>
>>> Stephen, I would like to be able to calculate the a, b, c in the eqn of
>>> an
>>> ellipsoid for the stars/dm particles in a Halo, and see how much mass
>>> is
>>> in particles, and how much is in baryons.  If I have the prescription
>>> to
>>> do the calculation, can this be easily achieved if I have access to the
>>> halo data if it was stored as binary in the .yt file?  Or would this be
>>> better done while it is still  in memory with the halo profiler
>>> running?
>>
>>
>> This is actually not so easy in yt presently.
>>
>>
>> First let me correct one of your misconceptions. There is no halo data
>> stored in the .yt file. Halo data is only ever stored in the
>> HopAnalysis.out text file, and the hdf5 .h5 particle data files (called
>> from halos.write_particle_lists()).
>>
>> What makes your task difficult is that yt has no ellipsoid data
>> container. It has a sphere, cylinder, inclined box, and rectilinear
>> region. You are more than welcome to write one for yourself! It's a bit
>> convoluted where you make the additions, but the trail beings in
>> yt/data_objects/data_containers.py. For example, you can follow how the
>> AMRSphereBase works. If this is something you want to try, I can write a
>> more detailed explanation of how to do this and the files involved. You
>> should also get yourself a bitbucket account and make your own fork of
>> yt if you are up to this challenge.
>>
>
> I did write up a short description of how to do this:
>
> http://yt.enzotools.org/doc/advanced/creating_datatypes.html
>
> I think I have the start of an ellipsoidal container from last fall,
> when Geoffrey was originally thinking about doing this.  It wasn't
> very much, however, as the project got sidelined by some other things.
>
>> As I see it you have two tasks. First, is the ability to define, or at
>> least discriminate, cells that are inside your ellipsoid. If you can do
>> this, this will get your baryonic content. Second, you need to decide
>> how you want to deal with particles. Halos are generally not
>> ellipsoidal. Do you want to find the halo, find the ellipsoid for that
>> halo, and then grab the particles inside that region, which will not be
>> completely identical to what the halo finder found? Or do you want to
>> sick with the identified particles?
>>
>> To answer your specific question, as it stands right now this would be
>> easiest to do as you find the halos. I am currently thinking about how
>> to import previously found halos back into yt. But I think you have
>> bigger challenges to figure out first: namely how to define an
>> ellipsoidal region in yt/python for your baryons.
>
> I agree; this is going to be the rate-limiting step, even if you do
> have the three vectors, three axis values and center corresponding to
> the description of the ellipse.
>
> It's probably best if you think about it in two stages:
>
> 1) How to tell if a grid intersects an ellipse
> 2) How to tell if a grid is completely enclosed in an ellipse
> 3) How to tell if a given x,y,z value is enclosed in an ellipse
>
> Once you have these written out, like Stephen says you should fork yt
> on BitBucket, add the code to the data_containers.py file, and we'll
> review it an pull it back in.
>
> Good luck!
>
> -Matt
>
>>
>>
>> Stephen Skory
>> stephenskory at yahoo.com
>> http://stephenskory.com/
>> 510.621.3687 (google voice)
>>
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>>
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>