;+
; NAME:
;	CenterOfMass
; PURPOSE:
;	Calculate center of mass of a distribution of masses in n dimensions
; CATEGORY:
;	gen/idl/toolbox
; CALLING SEQUENCE:
	FUNCTION CenterOfMass, mass_array, loc_array
; INPUTS:
;	mass_array		array of any structure; type: integer or float
; OPTIONAL INPUT PARAMETERS:
;	loc_array		array[n, size(mass_array,/dim)]
;						i.e. same structure as mass_array with one extra
;						leading dimension
; OUTPUTS:
;	P				scalar, or array[n]; double precision
;						coordinates of center of mass. In the one-dimensional
;						case (n=1) a scalar is returned.
; INCLUDE:
	@compile_opt.pro				; On error, return to caller
; CALLS:
;	SubArray, SyncDims
; SIDE EFFECTS:
;	The center of mass calculation is done in double precision.
; PROCEDURE:
;	Two separate cases are calculated:
; >	1. Only the mass_array is specified as an n-dimensional array.
;		The center of mass is returned as an n-element array (scalar if n=1).
;		The array index in each dimension is used as the distance scale:
;			mass_array = [[1,2],[3,4]]
;			P = centerofmass(mass_array)
;		results in
;			P = [0.6,0.7]
; >	2. If both mass_array and loc_array are specified then loc_array, with
;			leading dimension of n elements, specifies the n-dimensional locations
;			of the masses in mass_array, e.g. in the 3-dimensional case, if mass_array
;			is a 10x10 array then loc_array would be a 3x10x10 array, with
;			loc_array[0,*,*] the x-coordinate, loc_array[1,*,*] the y-coordinate and
;			loc_array[2,*,*] the z-ccordinate.
;		The above example again:
;			mass_array = [1,2,3,4]
;			loc_array  = [[0,0],[1,0],[0,1],[1,1]]
;			P = centerofmass(mass_array)
;		results in
;			P = [0.6,0.7]
; MODIFICATION HISTORY:
;	NOV-1999, Paul Hick (UCSD/CASS)
;	FEB-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		internally the procedure converts the mass array to a normalized mass by
;		dividing by the total mass. This results in a divide by zero if the sum is
;		zero (presumably because all masses are zero). To avoid this one unit is
;		added to each mass.
;-

mass = mass_array

IF total(mass, /nan) EQ 0 THEN mass += 1	; Presumably all zero masses

mass = mass/total(mass, /nan, /double)			; Mass normalized to total mass

IF n_elements(loc_array) EQ 0 THEN BEGIN

	s = size(mass)
	n = s[0]
	center = dblarr(n, /nozero)

	; For each coordinate i the mass array is collapsed by summing
	; over all dimensions except i. This combines all masses with
	; the same i-coordinate. The resulting array[s[1+i]] is multiplied
	; by the appropriate i-coordinate values to find the
	; position of the mean value in dimension i.

	FOR i=0,n-1 DO BEGIN
		submass = mass
		FOR j=n-1,0,-1 DO IF j NE i THEN submass = total(submass, 1+j, /nan, /double)
		center[i] = total(submass*dindgen(s[1+i]), /nan, /double)
	ENDFOR

ENDIF ELSE BEGIN

	s = size(loc_array)
	m = n_elements(mass)						; # of masses
	n = s[s[0]+2]/m								; # dimensions

	center = dblarr(n, /nozero)					; Center of mass coordinates

	; The use of SubArray to calculate the center of mass position
	; requires that the 1st dimension of loc_array is equal to n
	; (the # dimensions). This may not be the case for 'degenerate'
	; cases: both mass_array and loc_array might be scalars (s[0]=0)
	; or, in the one-dimensional case, mass_array and loc_array
	; probably would have the same array structure (with no leading
	; dimension of 1 on loc_array).

	degenerate = s[0] EQ 0 OR s[1] NE n

	IF degenerate THEN loc_array = reform([loc_array],n,m, /overwrite)

	FOR i=0,n-1 DO center[i] = total(mass*SubArray(loc_array,element=i), /nan, /double)

	IF degenerate THEN SyncDims, loc_array, size=s ; Reestablish original structure

ENDELSE

IF n EQ 1 THEN center = center[0]				; Return scalar for 1-dim case

RETURN, center  &  END