[Yt-dev] Radius Calculation

[Matthew Turk]

Oops, just as a note, there's an error in the paste -- the second plot
of "NewRadius" should be of "Radius".

Thanks guys...

-Matt

On Fri, Apr 17, 2009 at 7:16 PM, Matthew Turk <matthewturk at gmail.com> wrote:
> Hi guys,
>
> Can somebody review this?  I chatted with Jeff, and I believe it is
> correct for periodic boxes, but I need somebody to look at it before I
> commit it, as it is a rather crucial aspect of yt.  This changes the
> way radius is calculated and should now work with periodic boundary
> conditions.
>
> Here's the new field definition:
>
> def _NewRadius(field, data):
>    center = data.get_field_parameter("center")
>    DW = data.pf["DomainRightEdge"] - data.pf["DomainLeftEdge"]
>    radius = na.zeros(data["x"].shape, dtype='float64')
>    for i, ax in enumerate('xyz'):
>        r = na.abs(data[ax] - center[i])
>        radius += na.minimum(r, na.abs(DW[i]-r))**2.0
>    na.sqrt(radius, radius)
>    return radius
>
> the final argument in na.sqrt is the 'out' argument, to avoid temporary arrays.
>
> I've committed a sample script here:
>
> http://paste.enzotools.org/show/101/
>
> which loads up a data dump, calculates the radius from [0.1, 0.1, 0.1]
> at all points in the root grid, then plots the output.  It gives
> correct min & max values for the radius.
>
> If I can get somebody else to sign off on this, I will commit the
> necessary changes to _Radius and _ParticleRadius, as well as grid &
> point selection for spheres.  We'll then have a fully-periodic
> profiler, which will close #165.
>
> -Matt
>
> On Thu, Apr 16, 2009 at 11:39 AM, Matthew Turk <matthewturk at gmail.com> wrote:
>> Hi guys,
>>
>> I've been exploring the idea of changing radius to be correct for
>> periodic boxes.  Right now it is incorrect; each component (x,y,z)
>> should not have a distance greater than 0.5 * domain_size.  The
>> easiest way to do this would be:
>>
>> rx = min( abs( x - center_x ), abs( x - center_x - domain_x) )
>>
>> I think the best way to do this is to set up a NumPy ufunc
>>
>> http://docs.scipy.org/doc/numpy/user/c-info.beyond-basics.html#creating-a-new-universal-function
>>
>> that accepts three arrays, along with the center and the domain width,
>> and then returns the radius.  What do you all think?  Alternatively,
>> I'm thinking maybe just a regular function that gets the arrays and
>> returns one would be better; the ufunc machinery is a bit complicated
>> and I might get confused.  Once I come up with it, will somebody be
>> able to look over my work?
>>
>> -Matt
>>
>