[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

*Subject*: Re: Center of mass???*From*: "J.D. Smith" <jdsmith(at)astro.cornell.edu>*Date*: Tue, 09 Nov 1999 15:51:06 -0500*Newsgroups*: comp.lang.idl-pvwave*Organization*: Cornell University*References*: <3826E4D0.6E84F53B@mrc.ks.se> <MPG.12909a27cf5a46b198994e@news.frii.com> <38279336.C13C76F8@ucsd.edu> <MPG.12914882e623d683989951@news.frii.com>*Sender*: verified_for_usenet(at)cornell.edu (jts11 on vodka.tn.cornell.edu)*Xref*: news.doit.wisc.edu comp.lang.idl-pvwave:17194

David Fanning wrote: > > 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 > > Well, it seems to work in the simple-minded cases I > tried it on. :-) > > I'll write it up in an article if no one else has > serious objections. > Why stop at 3 dimensions? How about any dimension: function com, arr,DOUBLE=dbl s=size(arr,/DIMENSIONS) n=n_elements(s) tot=total(arr,DOUBLE=dbl) if keyword_set(dbl) then ret=dblarr(n,/NOZERO) else ret=fltarr(n,/NOZERO) for i=0,n-1 do begin tmp=arr for j=0,i-1 do tmp=total(tmp,1,DOUBLE=dbl) for j=i+1,n-1 do tmp=total(tmp,2,DOUBLE=dbl) ret[i]=total(findgen(s[i])*tmp,DOUBLE=dbl)/tot endfor return,ret end I would have posted this yesterday, but I couldn't help but feel there must be a better way to do it. Here's why: In the 3D example above, Paul calculates total(array,3) twice, which means once too many. It could have been saved between steps. The number of times total() is called on the array in this straightforward, brute-force method is n*(n-1) (for each of n dimensions total the array over the other n-1), not counting the final index total. And many of those are repeats to the exact same total() call with the exact same array! Since we needed to do all directional sub-totals, there must be a more efficient way, analogous to the 3D case of saving total(array,3). After a bit of scratching on paper, I found that I could do it with only (n-1)(n+2)/2 calls to total(), by working simultaneously from the last dimension backward and the first dimension forward, saving sub-totals which can be shared by subsequent calls in both cases. For 10 dimensions that becomes 54 vs. 90 calls to total() (though only 5 vs. 6 for 3 dimens -- 1 saved call as described above). The problem is how to code this model. Two reciprocating mutually recursive functions should do the trick, but I haven't had time to explore. Any takers? Anybody think they can save more calls to total()? This may be folly, since few will use it for more than 3D, but it's an interesting problem. JD P.S. Another wrinkle: The total()'s you'd prefer to save are the ones that do the most work... those which total the "largest" dimensions. So sorting by dimension size first of all might speed things up even more. -- J.D. Smith |*| WORK: (607) 255-5842 Cornell University Dept. of Astronomy |*| (607) 255-6263 304 Space Sciences Bldg. |*| FAX: (607) 255-5875 Ithaca, NY 14853 |*|

**References**:**Re: Center of mass???***From:*David Fanning

- Prev by Date:
**Re: Object Widgets** - Next by Date:
**openr,/xdr & readu,structure problems** - Prev by thread:
**Re: Center of mass???** - Next by thread:
**Re: where() and IEEE NaN and other stuff...** - Index(es):