!  Function to compute weighted average of 3-vectors, as for center of
!    mass of a system of particles.
!
!  File:  center.f90

!  Calculate center of points X(n, 3) with masses W(n).  Result is a
!  3-vector, the center.

      function center(n, X, W)
      integer n			! Array dimension
      real    X(n, 3)		! The points
      real    W(n)              ! The masses
      real    center(3) 	! The result

      center = SUM( SPREAD(W, 2, 3)*X, 1 )/SUM(W)
      end function center

!  ============================================================================

!  Simple program to test `center' using artificial data.  These data are
!  viewed as positions and velocities for a collection of particles.

!  Test data:
!
!     Positions		Velocities	Masses
!     1    -1    2	 2  4  6	  2
!     2    -2    4	 2  4  6	  2
!     .  .  .		  . . .		. . .
!   512  -512 1024	 2  4  6	  2
!   513  -513 1026	 2  4  6	  1
!     .  .  .		  . . .		. . .
!  1024 -1024 2048	 2  4  6	  1

      program ctest

      INTERFACE
      function center(n, X, W)
      integer n			! Array dimension
      real    X(n, 3)		! The points
      real    W(n)              ! The masses
      real    center(3)   	! The result
      end function center
      END INTERFACE

      real    X(1024, 3), V(1024, 3)	! Positions, velocities
      real    M(1024)		        ! Masses
      real    cX(3), cV(3)        	! Center pos, vel
      integer I(1024, 3), J(1024, 3)	! Coordinates

!   Set up test data

      I = SPREAD( (/ (k, k=1,1024) /), 2, 3)	! Row
      J = SPREAD( (/ (k, k=1,3) /), 1, 1024)	! Column
      X = I*SPREAD( (/ 1, -1, 2 /), 1, 1024)	! Simple position data
      V = 2*J					! All have same velocity
      M = (/ (2, k=1,512), (1, k=1,512) /)	! Masses
      print *, X(1,:)				! Print representative
      print *, X(2,:)				!   entries
      print *, V(1,:)
      print *, M(512), M(513)

!   Calculate centers and report
      print *, "Two center calculations on 1024 particles:"
      cX = center(1024, X, M)
      cV = center(1024, V, M)
      print *, "cX = ", cX
      print *, "cV = ", cV

      end program ctest

