!  N-body calculation using the modified Euler method,
!  which is globally 2nd-order.
!  Uses a ring organization to calculate accelerations.
!
!  Data arrays (position, velocity, acceleration, mass)
!  are organized:
!
!          1    Particle number (i) [axis 1] -->      N
!	  +--------------------------------------------+
!   |  x  |					       |
!  (c) y  |					       | [axis 2]
!   v  z  |					       |
!	  +--------------------------------------------+
!
!  The processor array, parameter N, arrays
!
!	X0, X1	- Current, next positions
! 	V0, V1	- velocities
!	 M, GM	- mass, G*mass (redundant)
!	 c,  i	- array coordinates
!  and
!	    NB	- integer <= N, actual no. of bodies
!	     G	- real, the gravitational constant
!
!  are defined in the file nbrmod.f90
!
!  Source file:  nbring.f90	Module with global data:  nbrmod.f90
!  Acceleration:  grav.f90, etc.


!------- setup() -----------------------------------------------------

      subroutine setup()
      USE nbrmod

!  Initialize coordinates
      c = SPREAD( [1:3],  1, N)
      i = SPREAD([1:N], 2, 3)
      end subroutine setup


!------- init() ------------------------------------------------------

!  Load initial positions, velocities, also masses,
!  gravitational const.  All are read from logical
!  unit 10.  Replicates masses over x,y,z.

      subroutine init()
      USE nbrmod

      integer k

      X0 = 0.0		! For unused particles
      V0 = 0.0
      do k = 1, NB
        read (10, *) X0(k, 1), X0(k, 2), X0(k, 3)
      end do
      do k = 1, NB
        read (10, *) V0(k, 1), V0(k, 2), V0(k, 3)
      end do
      M = 0.0		! For unused particles
      do k = 1, NB
        read (10, *) M(k)
      end do
      read (10, *) G
      end subroutine init


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

!  Insert nx asteroids at radius ar with total mass am.
!  This assumes that an orbit of radius 1 about origin
!  corresponds to speed 2*pi (e.g., au-earth-yr units).

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

      real, parameter:: TWOPI = 6.2831853
      real, array(1:N) :: sx, cx, arg
      real spd
      integer nn, ii, cc

      print 21, nx, ar, am
 21   FORMAT (" Making ", I4, " asteroids, radius ", &
              F10.5, ", mass ", F8.5)
      nn = NB + nx
      spd = TWOPI/sqrt(ar)
      do ii = NB+1, nn
	arg(ii)  = (TWOPI*(ii - NB))/nx
        sx(ii)   = sin(arg(ii));  cx(ii)   = cos(arg(ii))
	X0(ii,1) =  cx(ii)*ar;    V0(ii,1) = -sx(ii)*spd
	X0(ii,2) =  sx(ii)*ar;    V0(ii,2) =  cx(ii)*spd
	X0(ii,3) = 0.0;           V0(ii,3) = 0.0
      end do
      M(NB+1:nn) = am/nx
      NB = nn
      end subroutine mkast


!------- print_state(PX) ---------------------------------------------

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

      subroutine print_state(PX)
      USE nbrmod
      integer PX	! Print control

      integer k, np

      np = NB
      if ( (PX >= 0) .and. (PX < NB) ) np = PX
      if ( np > 0 ) print 21
 21   FORMAT (/ "         x         y         z     ", &
                "    vx        vy        vz        m")
      do k = 1, np
        print 22, k, X0(k,1), X0(k,2), X0(k,3), &
                     V0(k,1), V0(k,2), V0(k,3), M(k)
      end do
 22   FORMAT (1X, I3, 6F10.5, E10.2)
      end subroutine print_state


!------- print_x(ax, na, nf) ------------------------------------------

!  Print ax(na:nf, :), nf-na <= 9 (for testing)

      subroutine print_x(ax, na, nf)
      USE nbrmod
      real, array(1:N,1:3) :: ax
      integer na, nf

      integer u, lm

      lm = min(nf, na+9)
      do u = 1, 3
        print 21, ax(na:lm, u)
 21     FORMAT ( 10(1X : E9.2 :) )
      end do
      end subroutine print_x

!------- print_c(Xf, Xc, Af, Ac, na, nb) -------------------------------

!  Print conveyor status in positions na:nb (for testing)

      subroutine print_c(Xf, Xc, Af, Ag, na, nf)
      USE nbrmod
      real, array(1:N,1:3) :: Xf, Xc, Af, Ag
      integer na, nf

      integer kk
      print 22, (kk, kk = na, nf)
 22   FORMAT ( 10(1X : I9 :) )
      print *, "  - Xf ----------"
      call print_x(Xf, na, nf)
      print *, "  - Xc ----------"
      call print_x(Xc, na, nf)
      print *, "  - Af ----------"
      call print_x(Af, na, nf)
      print *, "  - Ac ----------"
      call print_x(Ag, na, nf)
      end subroutine print_c


!------- modEuler(h, ns)  -----------------------------------------


!  Perform ns steps of the modified Euler method with time-step h.
!  X0, V0 are taken to be current state, and are updated
!  with final state.  X1, V1 are used as working arrays.

      subroutine modEuler(h, ns)
      USE nbrmod
      IMPLICIT NONE
      real h, tol
      integer ns

      INTERFACE
	subroutine Grav(X)
        USE nbrmod
        real, dimension(1:N, 1:3) :: X
	end subroutine Grav
      END INTERFACE

      real h2
      integer k

      h2 = 0.5*h	! h/2, for convenience
      do k = 1, ns	! ns steps
	call Grav(X0)   	! current acceleration 
        V1 = V0 + h2*Ac         ! half-way velocity
	X1 = X0 + h2*V0         ! half-way position
        call Grav(X1)           ! half-way acceleration
	V0 = V0 + h*Ac		! New velocities
	X0 = X0 + h*V1     	! New positions
      end do    	! ns steps
      end subroutine modEuler


!------- nbControl  ----------------------------------------------------------

!  Main program:  reads problem and control information
!  from data file named by user, performs calculation.
!  Print control:
!	  0 - just performance data
!	< 0 - state for all particles
!	> 0 - state for number of particles specified

      program nbControl
      USE nbrmod
      IMPLICIT NONE

      INTERFACE
      subroutine modEuler(h, ns)
      real h, tol
      integer ns
      end subroutine modEuler
      END INTERFACE

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

      print '(" Data file> "$)'    ! Request data file
      read *, fnm
      print '(" Print ctl> "$)'	! Request print control
      read *, PX
!   Begin initialization
      open(10, file=fnm, status = "OLD")
      read (10, *) NB
      call setup()     ! Setup geometric structure
      call init()      ! Read data for individual bodies
!   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
!   Make extra asteroids as needed to fill arrays.
      if ( NB < N ) call mkast(N - NB, 3.4, 0.05)
      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.0"
      call print_state(PX)
!   Scale masses
      GM = G*m
!   np runs of ns steps, timestep dt
      do k = 1, np
        call modEuler(dt, ns)
        print 23, k*ns*dt
 23     FORMAT (/ " simulated time = ", F9.3)
        print 24, Ac(1,:)                    ! For testing
 24     FORMAT ("Ac(1,:):  ", 3(1X, F10.6))
	call print_state(PX)
      end do
      end program nbControl

!*******************************************************
!  Data files:
!
!   N			No. of explicit bodies
!
!   x y z		Position of body 1
!     :                      :
!   x y z		Position of body N
!
!   vx vy va		Velocity of body 1
!      :                     :
!   vx vy va		Velocity of body N
!
!   m			Mass of body 1
!   :			      :
!   m			Mass of body N
!
!   G			Gravitational constant
!
!   nx			Number,
!   ar am		radius, mass of asteroid group
!      :                        :     0 or more groups
!   nx			Number,
!   ar am		radius, mass of asteroid group
!   -1			Stopper
!
!   ns np dt		Number of time steps per phase,
!   			  number of phases,
!   			  delta-t (per time step)
!
!  Further data in the file will be ignored.
!*******************************************************
