C mpfortran -pesize=8 -nofecommon nb2.f
C
C----ROB SALGADO (cps615g)
C    nb2.f  11/09/92 (extended from 11/06/92 PROBLEM 5.1)
C
C
C--NEW AND IMPROVED-OVER-LAST-WEEK'S-HORRIBLE-ATTEMPT-TO-IMPLEMENT-RUNGE-KUTTA
C--INCLUDES CALCULATIONS OF THEORETICALLY-CONSERVED QUANTITIES


C  N-body calculation using the global 4-th order Runge-Kutta
C  Uses O(N^2) processors to calculate accelerations.
C
C  Data arrays (position, velocity, acceleration, mass)
C  are (N, N, 3), organized:
C
C            ________________________________
C     /|  z /_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/|
C   (c)  y /_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ |
C   /   x /_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ /|
C       1 |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   |     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   |     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   P     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   a     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   r     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   t     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   i     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   c     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   l     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   e     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C         |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C  (i)    |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C   |     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|///
C   v     |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//
C       N |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|/
C
C          1   ---  Particle (j)  -->    N
C
C  The parameter N and common block /NBODYPE/
C  with arrays
C
C       X0, X1  - Current, next positions
C       V0, V1  - velocities
C        MI,GJ  - mass(i), G*mass(j)
C        i,j,c  - array coordinates
C         diag  - diagonal entries
C          act  - active in acceleration 
C
C  and common block /NBODYFE/ with
C
C           NB  - integer <= N, actual no. of bodies
C            G  - real, the gravitational constant
C
C  are defined in the file nb2mode.h
C
C  System wide conventions regarding data:
C
C       Entries (i, :, :)  - data for particle i
C  except    MJ (:, j, :) = G*mass(j)
C
C  This is intentionally redundant, and the programs
C  perform many redundant calculations.
C
C  File:  nb2mode.f



C------- setup() ---------------------------------------

C  Initialize coordinates, logical flags.

      subroutine setup()

      INCLUDE 'nb2mode.h'
      i = SPREAD(SPREAD([1:N], 2, N), 3, 3)
      j = SPREAD(SPREAD([1:N], 1, N), 3, 3)
      c = SPREAD(SPREAD([1:3], 1, N), 2, N)
      diag = (i == j)
      actv = .FALSE.
      end subroutine setup



C------- init() ----------------------------------------

C  Load initial positions, velocities, also masses,
C  gravitational const.  All are read from logical
C  unit 10.  Replicates masses over x,y,z.
C  Data for particle i are left in (i, :, :) .

      subroutine init()

      INCLUDE 'nb2mode.h'

      integer k
      real mx
      X0 = 0.0          ! For unused particles
      V0 = 0.0
      do, k = 1, NB
        read (10, *) X0(k, 1, 1), X0(k, 1, 2), X0(k, 1, 3)
      end do
      do, k = 1, NB
        read (10, *) V0(k, 1, 1), V0(k, 1, 2), V0(k, 1, 3)
      end do
      M = 0.0           ! For unused particles
      do, k = 1, NB
        read (10, *) mx
        MI(k, 1, :) = mx
      end do
      read (10, *) G
      X0 = SPREAD(X0(:, 1, :), 2, N)
      V0 = SPREAD(V0(:, 1, :), 2, N)
      MI = SPREAD(MI(:, 1, :), 2, N)
      end subroutine init


C------- mkast(nx, ar, am) -----------------------------

C  Insert nx asteroids at radius ar with total mass am.
C  This assumes that an orbit of radius 1 about origin
C  corresponds to speed 2*pi (e.g., au-earth-yr units).
C  Data for particle i are left in  (i, :, :) .

      subroutine mkast(nx, ar, am)
      integer nx        ! Number of new asteroids
      real ar, am       ! Orbit radius, total mass

      INCLUDE 'nb2mode.h'

      real, parameter:: TWOPI = 6.2831853
      real, array(1:N, 1:N, 1:3) :: sx, cx, arg
      real spd
      integer nn
      print 21, nx, ar, am
 21   FORMAT (" Making ", I4, " asteroids, radius ",
     *        F10.5, ", mass ", F8.5)
      nn = NB + nx
      spd = TWOPI/sqrt(ar)
      where ( (NB < i) .and. (i <= nn) )
        arg = (TWOPI*(i - NB))/nx
        sx = sin(arg)
        cx = cos(arg)
      end where
      where ( (NB < i) .and. (i <= nn) .and. (c == 1) )
        X0 =  cx
        V0 = -sx
      end where
      where ( (NB < i) .and. (i <= nn) .and. (c == 2) )
        X0 =  sx
        V0 =  cx
      end where
      where ( (NB < i) .and. (i <= nn) .and. (c == 3) )
        X0 = 0.0
        V0 = 0.0
      end where
      where ( (NB < i) .and. (i <= nn) )
        X0 = X0*ar
        V0 = V0*spd
        MI = am/nx
      end where
      NB = nn
      end subroutine mkast


C------- print_state(PX,now) -------------------------------MODIFIED BY ROB
C                                                           IN ORDER TO GET
C                                                           CURRENT TIME FOR
C                                                           CALCULATING Cc

C  Print current state of system according to PX.
C  See comments for main program for PX settings.

      subroutine print_state(PX,now)
      integer PX        ! Print control
      real now                                 ! Pass the global clock-reading
      INCLUDE 'nb2mode.h'
      real nrg       ! energy
      real, array(1:N,1:3) :: VT
      real, array(1:N,1:N,1:3) :: Ac,D,R,XT
      integer k, np
      np = NB
      if ( (PX >= 0) .and. (PX < NB) ) np = PX
      print 21
 21   FORMAT (/ "         x         y         z     ",
     *          "    vx        vy        vz        m")
      do, k = 1, np
        print 22, k, X0(k,1,1), X0(k,1,2), X0(k,1,3),
     *               V0(k,1,1), V0(k,1,2), V0(k,1,3), MI(k,1,1)
      end do

C-------------------------------CALCULATE (theoretically) CONSERVED QUANTITIES

C   ---LINEAR MOMENTUM     Cm:=SUM m[i]*V[i]
      VT=SUM(MI*V0,1)
      print 220,  VT(1,1), VT(1,2), VT(1,3),
     * SQRT(VT(1,1)*VT(1,1)+VT(1,2)*VT(1,2)+VT(1,3)*VT(1,3))

C   ---INIT COM-POSITION   Cc:=Cm*t-SUM m[i]*X[i]
      VT=VT*now-SUM(MI*X0,1)
      print 221,  VT(1,1), VT(1,2), VT(1,3),
     * SQRT(VT(1,1)*VT(1,1)+VT(1,2)*VT(1,2)+VT(1,3)*VT(1,3))

C   ---ANGULAR MOMENTUM    Ca:=SUM m[i]*(X[i] x V[i])
      VT=SUM(MI*( CSHIFT(X0*CSHIFT(V0, 1,3), 1,3)
     *              -CSHIFT(X0+CSHIFT(V0,-1,3),-1,3) ),1)
      print 222,  VT(1,1), VT(1,2), VT(1,3),
     * SQRT(VT(1,1)*VT(1,1)+VT(1,2)*VT(1,2)+VT(1,3)*VT(1,3))

C   ---KINETIC ENERGY 
      VT=SUM(MI*V0*V0,1)
      nrg=(VT(1,1)+VT(1,2)+VT(1,3))/2.0
C   ---POTENTIAL ENERGY (variant of Grav(Xf))---
      where (diag)
           XT=X0
      elsewhere
           XT=-1E30
      end where
      XT=SPREAD(MAXVAL(XT,1),1,N)
      Ac=0.0
      where (actv)
           D=XT-X0
      elsewhere
           D=1.0
      end where
      R=sqrt(SPREAD(SUM(D*D,3),3,3))
      where (actv) Ac=GJ/R
      VT=SUM(Ac,2)
      print 223,nrg+SUM(VT)/(6.0)   ! 6=(overcount)2 * (triplecount)3
      print 224

 22   FORMAT (1X, I3, 6F10.5, E10.2)
 220  FORMAT ("LIN MOM   (", 3F12.5, ") ", F12.5)
 221  FORMAT ("INIT COM  (", 3F12.5, ") ", F12.5)
 222  FORMAT ("ANG MOM   (", 3F12.5, ") ",F12.5)
 223  FORMAT (
     * "ENERGY                                           ", F12.5)
 224  FORMAT ("-- -- --")
      end subroutine print_state


C------- print_x(ax, nx, cx) -------------------------------

C  Print ax(1:nx, 1:nx, cx), nx <= 10 (for testing)

      subroutine print_x(ax, nx, cx)

      INCLUDE 'nb2mode.h'

      real, array(1:N, 1:N, 1:3) :: ax
CMPF  ONDPU ax
      integer nx, cx

      integer u, lm
      lm = min(nx, 10)
      do, u = 1, lm
        print 21, ax(u, 1:lm, cx)
 21     FORMAT ( 10(1X : E9.2 :) )
      end do
      end subroutine print_x


C------- Grav(Xf) --------------------------------------

C  Calculate accelerations (Newtonian gravitation)
C  Uses completely parallel calculation, ignoring
C  anti-symmetry of force.  Acceleration for body i is
C  left in  (i, :, :)  of result.
C  Acceleration on body i due to body j is calculated
C  in entries (i, j, :) .

      function Grav(Xf)

      INCLUDE 'nb2mode.h'

      real, array(1:N, 1:N, 1:3) :: Xf, Grav
CMPF  ONDPU Xf

      real, array(1:N, 1:N, 1:3) :: Ac, D, R, Xj
      integer k, lm

C   Spread X0 across j from diagonal:
      where ( diag )
        Xj = Xf
      elsewhere
        Xj = -1E30              ! Enormously negative
      end where
      Xj = SPREAD(MAXVAL(Xj, 1), 1, N)
D     print *, " --- Xj -----------"
D     call print_x(Xj, 7, 1)    ! For testing
C   Now Xj(i,j,:) = X0(j,j,:)

      Ac = 0.0
      where ( actv )
        D  = Xj - Xf            ! Displacement
      elsewhere
        D = 1.0            ! Protects from div by 0
      endwhere
D     print *, " --- D  -----------"
D     call print_x(D, 7, 1)     ! For testing
      R  = sqrt(SPREAD(SUM(D*D, 3), 3, 3))
D     print *, " --- R  -----------"
D     call print_x(R, 7, 1)     ! For testing
D     print *, " MIN(R) = ", MINVAL(R)
      where ( actv )  Ac = GJ*D/R**3
D     print *, " --- Ac  -----------"
D     call print_x(Ac, 7, 1)    ! For testing
      Grav = SPREAD(SUM(Ac, 2), 2, N)
      end function Grav


C------- modeul(h, ns)  --------------------------------

C  Perform ns steps of modified Euler with time-step h.
C  X0, V0 are taken to be current state, and are updated
C  with final state.  X1, V1 are used in the process.
C  Note that X's, V's are calculated redundantly.

      subroutine modeul(h, ns)
      real h
      integer ns

      INCLUDE 'nb2mode.h'

      INTERFACE
        function Grav(Xf)
        INCLUDE 'nbmode.h'
        real, array(1:N, 1:N, 1:3) :: Xf, Grav
CMPF    ONDPU Xf
      END INTERFACE

C----------------------FOR NOW, JUST BETWEEN YOU AND ME,
C                      WE WILL USE THE 4-TH ORDER RUNGE-KUTTA
C                      IF YOU CAN READ THIS MESSAGE, THEN THAT MEANS
C                      I HAVEN'T RENAMED modeul() TO, SAY, rk4()

      real h2
      integer k

      real, array(1:N,1:N,1:3) :: k1,k2,k3,k4

      h2 = 0.5*h        ! h/2, for convenience
      do k = 1, ns      ! ns steps
D     print *, " --- X0  -----------"
D     call print_x(X0, 7, 1)    ! For testing

C       (Modified Euler)
      X1 = X0 + h2*V0         ! Half-way positions
      V1 = V0 + h2*Grav(X0)   ! Half-way velocities

      k1=h*V0                                    ! X'=V
      k2=h*(V0+h2*Grav(X0))
      k3=h*(V0+h2*Grav(X1))
      k4=h*(V0+h*Grav(X0))
      X1=X0+(k1+k4+2.0*(k2+k3))/6.0              ! this will be the new X

      k1=h*Grav(X0)                              ! V'=Grav(X)
      k2=h*(Grav(X0+h2*V0))
      k3=h*(Grav(X0+h2*V1))
      k4=h*(Grav(X0+h*V0))                       ! this is the new 

      V0=V0+(k1+k4+2.0*(k2+k3))/6.0              ! this is the new V
      X0=X1                                      ! this is the new X

C       (Modified Euler)
C       X1 = X0 + h2*V0         ! Half-way positions
C       V1 = V0 + h2*Grav(X0)   ! Half-way velocities
C       X0 = X0 + h*V1          ! Next positions
C       V0 = V0 + h*Grav(X1)    ! Next velocities

      end do
      end subroutine modeul


C------- nb2mode  -------------------------------

C  Main program:  reads problem and control information
C  from data file named by user, performs and times
C  calculation.  Print control:
C         0 - just performance data
C       < 0 - all
C       > 0 - number given

      program nb2mode

      INCLUDE 'nb2mode.h'

      integer ns, np, nv, k, PX, mptm
      real dt, ar, am
      character*50 fnm

      print '("/ Data file> "$)'    ! Request data file
      read *, fnm
      print '("  Print ctl> "$)'        ! Request print control
      read *, PX
C   Begin initialization
      open(10, file=fnm, status = "OLD")
      read (10, *) NB
      call setup()     ! Setup geometric structure
      call init()      ! Read data for individual bodies
C   Read "asteroid" descriptions from file
      read (10, *) ns
      do while ( ns > 0 )
        read (10, *) ar, am
        call mkast(ns, ar, am)
        read (10, *) ns
      end do
      read (10, *) ns, np, dt   ! Integration parameters
      close(10)
      print 21, NB, N
 21   FORMAT (/, X, I4, " bodies, geometry allows ", I4)
      print 22, G, dt
 22   FORMAT ("    G = ", E13.6, "  delta-t = ", E12.5)
      print *, "  simulated time = 0"
      call print_state(PX,0.0)
C   Scale masses and distribute, set actv
      where ( diag )
        GJ = G*MI
      elsewhere
        GJ = -1E30              ! Enormously negative
      end where
      GJ = SPREAD(MAXVAL(GJ, 1), 1, N)
D     print *, " --- GJ  -----------"
D     call print_x(GJ, 7, 1)    ! For testing
      actv = (i <= NB) .and. (j <= NB) .and. .not. diag
D     print *, " ---actv------------"
D     do, k = 1, 7
D       print *, actv(k, 1:7, 1)
D     end do
C   np runs of ns steps, timestep dt
      do, k = 1, np
        call mpTimerStart()
        call modeul(dt, ns)
        mptm = mpTimerElapsed()
        print 23, k*ns*dt, mptm/1000.0
 23     FORMAT (/ " simulated time = ", F9.3,
     *            "  Elapsed time = ", F9.3, " msec.")
        call print_state(PX,k*ns*dt)
      end do
      end program nb2mode

C*******************************************************
C  Data files:  See nbmode.f
C*******************************************************
