!  N-body calculation for general multi-step integrators to 6-th order.
!  Uses O(N^2) processors to calculate accelerations.  Integration is
!  started by one-step method as for nb1step.fcm, with SS steps of starter
!  per step of multi-step integrator.  Parameters and basic common data
!  taken from nb1step.h.  Additional common data for multistep method is
!  defined in nbmstep.h.  Acceleration function is the same as for nb1step.fcm.
!
!  File:  nbmstep.f90		Modules:  nbmsmod.f90, nbrmod.f90
!
!  Multi-step integrator:
!  Perform ns steps with time-step h.  X0, V0 are taken to be current
!  state, and are updated with final state.  Velocities and accelerations
!  from earlier times are stored in VH, AH (see below).
!  X1, V1 may be used in the process.
!  X's, V's are calculated redundantly.
!
!     INTERFACE
!     subroutine Multistep(h, ns)
!     real h
!     integer ns
!     END INTERFACE
!

!  Data arrays (position, velocity, acceleration, mass) are (N, 3),
!  as for nbring.f90.  History of velocities, accelerations is stored
!  in VH, AH, each (NH, N, 3).
!
!  The module nbrmod.f90 is used as in nbring.f90.
!
!  The module nbmsmod.f90 defines the parameter
!
!	    NH		- max number of "historical" points
!
!  the real arrays
!
!	VH, AH	- Velocities, accelerations from earlier steps
!
!  and the integer array
!
!	    HI(1:NH)	- subscripts for V, H for earlier steps.
!
!  System wide conventions regarding data:
!
!	VH(HI(u),:,:) has V for u steps earlier, similarly AH
!


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

!  Initialize coordinates, logical flags.

      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


!------- Startup --------------------------------

!  Initialize NN steps of h, using SS steps of the one-step integrator
!  for each step of dt.  To make it easy to use standard integrators,
!  we calculate some accelerations more than once.  This is very small
!  overhead for high-order integrators with SS > 1.

      subroutine Startup(h, NN, SS)
      USE nbrmod
      USE nbmsmod
      real h
      integer NN, SS

      INTERFACE
      subroutine Onestep(h, S)
      USE nbrmod
      IMPLICIT NONE
      real h
      integer S
      END INTERFACE

      real  hs
      integer u, uh

      hs = h/SS
      do, u = 1, NN
	uh     = NN+1-u
	HI(uh) = uh
	VH(uh,:,:) = V0
        call Grav(X0);	AH(uh,:,:) = Ac	! Calculated again in Onestep
	call Onestep(hs, SS)		! Advances h
!        print 23, u*h                   !Diagnostic
! 23     FORMAT (/ " simulated time = ", F9.3)
!        print 24, Ac(1,:)                    ! For testing
! 24     FORMAT ("Ac(1,:):  ", 3(1X, F10.6))
!        call print_state(4)             ! Diagnostic
      end do
      end subroutine Startup


!------- nbmstep  -------------------------------

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

      program nb1step
      USE nbmsmod
      USE nbrmod

      integer ns, np, nv, k, PX, NN, SS
      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
      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)
!   Scale masses and distribute, set actv
      Gm = G*M
!   Get starting information
      print '(" NN SS > "$)'
      read *, NN, SS
      if ( NN > NH ) then
	print *, "Too many start-up steps for this program."
	go to 99					! Exit
      end if ! NN > NH
!   Start multistep and run first group.
	call Startup(dt, NN, SS)
	call Multistep(dt, ns - NN)
	print 23, ns*dt
!        print 24, Ac(1,:)                    ! For testing
	call print_state(PX)
!   Finish np runs of ns steps, timestep dt
      do, k = 2, np
        call Multistep(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
 99   continue						! For exit
      end program nb1step

!*******************************************************
!  Data files:  See nbring.f90
!*******************************************************
