Re: Center of mass???

["Dick Jackson" <dickjackson(at)dfanning.com>]

Paul Hick (pphick@ucsd.edu) writes:

> Reasoning by analogy to the 2D case, this should work, I think:
>
> xcm = Total( Total(Total(array,3),2) * Indgen(s[0])) / totalMass
> ycm = Total( Total(Total(array,3),1) * Indgen(s[1])) / totalMass
> zcm = Total( Total(Total(array,2),1) * Indgen(s[2])) / totalMass

Right, but as a wise man once told me (Dr. Coyote, or something like that),
the fastest dimension to run across with things like Total is usually the
second dimension. Empirical testing confirms this, so I propose the
following, now extended to do 2D or 3D. I also pulled the "/ totalMass"
inside a bit to keep the numbers closer to 1, lessen the possibility of
overflow and perhaps maintain more precision.

FUNCTION CenterOfMass, array

s = Size(array, /Dimensions)
totalMass = Total(array)

CASE Size(array, /N_Dimensions) OF
   2: BEGIN
      xcm = Total(Total(array,2) / totalMass * Indgen(s[0]))
      ycm = Total(Total(array,1) / totalMass * Indgen(s[1]))
      Return, [xcm, ycm]
   END
   3: BEGIN
      totalAcross2 = Total(array, 2)    ; 2 is fastest dim to total across
      xcm = Total(Total(totalAcross2,   2) / totalMass * Indgen(s[0]))
      ycm = Total(Total(Total(array,1), 2) / totalMass * Indgen(s[1]))
      zcm = Total(Total(totalAcross2,   1) / totalMass * Indgen(s[2]))
      Return, [xcm, ycm, zcm]
   END
ENDCASE

END


; Time testing was done as follows:

array = findgen(200, 200, 200)
t0 = systime(1)
print, CenterOfMass(array)
print, systime(1)-t0

My timings went from 1.2 seconds (with the previous approach) down to 0.8.

I'd love to see the CASE statement disappear. Who will dare to generalize
this to N dimensions, while ensuring that we total over dimension 2 wherever
possible?

--
Cheers,
-Dick

Dick Jackson             Fanning Software Consulting, Canadian Office
djackson@dfanning.com    Calgary, Alberta   Voice/Fax: (403) 242-7398
    Coyote's Guide to IDL Programming: http://www.dfanning.com/