C-----------------------------------------------------------------------------
C  ROB SALGADO (cps615g)
C-----------------------------------------------------------------------------
C
C    *** SORRY FOR THE DELAY....  I HAD TO GO FISH OUT ALL OF MY SOURCE-CODE
C    *** AND FIGURE OUT WHAT IS WHAT AND HOW THEY ARE SUPPOSED TO BE COMPILED
C
C 2.f         PROBLEM 5.2   
C jacobi.f    PROBLEM 5.3  
C nb2.f       PROBLEM 5.1   
C prob2.f     PROBLEM 4.2   (not compilable/executable.. see nb2.f)
C sph.f       PROBLEM 4.1   
C sorshift.f  PROBLEM 6.1
C
C-----------------------------------------------------------------------------
C--------------these are the 'heads' of each included file--------------------
C mpfortran -pesize=8 -d_lines 2.f
C mpfortran -pesize=8 -d_lines jacobi.f
C mpfortran -pesize=8 -nofecommon nb2.f
C echo 'prob2.f is not compilable/executable source code see:  nb2.f'
C mpfortran -pesize=8 sorshift.f
C mpfortran -pesize=8 -d_lines mpfrand.o sph.f 
C-----------------------------------------------------------------------------


C=============================================================================
C=============================================================================
C mpfortran -pesize=8 -d_lines 2.f
C
C----ROB SALGADO (cps615g)
C    2.F 11/12/92 (extended from 11/06/92 PROBLEM 5.2) 
C
C    MORE-SOPHISTICATED INITIAL-VALES
C
      subroutine showme(U)
      integer, parameter :: R=128,C=64
      real, array(R,C) :: U
      integer i,j
      print *,"------------------BOUNDARY POINTS + SAMPLE OF INTERIOR--"
      print 199, [1:C-1:int((C-1)/7)],C
      do, i=1, R-1, int((R-1)/11)
         print 200,i,U(i,1:C-1:int((C-1)/7)),U(i,C)
      end do
      print 200,R,U(R,1:C-1:int((C-1)/7)),U(R,C)
      print *,"-------------------------------------------------------/"
      print *," "
 199  FORMAT("    ",8I8)
 200  FORMAT(I4,8F8.4)
      end subroutine showme

C  Elementary solution of Laplace's equation in 2 dimensions
C  using the red/black SOR iteration.
C  Parameters R, C give array dimensions

C  Perform n cycles of SOR on U with relaxation parameter W
C  Returns updated function, writes error estimate in mxr
C  Error estimate is max abs change over all n iterations.

      function Sor(U, n, W, offset)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: Sor, U
      integer n
      real W, offset
      real maxU1, maxU0
      COMMON maxU1, maxU0

      real, array(R, C) ::  U1, U2, UU
      logical, array(R, C) :: Red
      integer k
      real OnemW                        ! 1 - W

      OnemW = 1.0 - W
C  Test for coord difference i - j being even
      Red = Iand( SPREAD([1:R], 2, C) - SPREAD([1:C], 1, R),
     *            1) == 0
      U1 = U
      U2 = U
      do, k = 1, n
         UU=U
       Where ( Red(2:R-1, 2:C-1) )
        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
     *                           + U1(3:R  , 2:C-1)
     *                           + U1(2:R-1, 1:C-2)
     *                           + U1(2:R-1, 3:C  ) )
        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
     *                    + OnemW*U1(2:R-1, 2:C-1)
       Elsewhere
        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
     *                           + U1(3:R  , 2:C-1)
     *                           + U1(2:R-1, 1:C-2)
     *                           + U1(2:R-1, 3:C  ) )
        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
     *                    + OnemW*U1(2:R-1, 2:C-1)
       end Where
      U1=U1-offset*UU
       maxU0=maxU1
       maxU1 = maxval(abs(U1))  ! maximum absolute value obtained
C  Re-normalize to the maximum value
      U1= U1/maxU1
      end do
C  Re-normalize to the maximum value
      Sor = U1
      print *,"Current:                    ",maxU1
      end function Sor


            
C  Test Redblack SOR

      program TestRedblack

      INTERFACE
        function Sor(U, n, W, offset)
        integer, parameter :: R = 128, C = 64
        real, array(R, C) :: Sor, U
CMPF    ONDPU Sor,U
        integer n
        real W, offset
      END INTERFACE

      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: U0, Z
      integer, array(R, C) :: I, J
      real bsum, eps, err, W, epsx, offset, maxU0, maxU1
      integer it, nits, notsure
      COMMON maxU1, maxU0
C--Ask for number of iterations
C
      print *, "How many NITS? "
      read *,nits

      I = SPREAD( [1:R], 2, C )         !Coordinates
      J = SPREAD( [1:C], 1, R )

C     Z is a solution of the difference equation

      Z = -1.4 + 0.12*I - 0.23*J + 0.001*(I**2 - J**2)
      U0 = Z    ! Install boundary values
C      print  *, "ORIGINAL ANALYTIC SOLUTION... ZERO BCs IMPOSED LATER"
C      call showme(U0)

       bsum = SUM(U0(1,:)) + SUM(U0(R,:)) ! Sum on boundary
     *     + SUM(U0(:,1)) + SUM(U0(:,C))
C
CC     Intialize interior to average on boundary
CC     The "-4" prevents counting corners twice
C
       U0(2:R-1, 2:C-1) = bsum/(2*R + 2*C - 4)
C       print *,"INITIAL GUESS: UNSOPHISTICATED"
C       call showme(U0)
      Z(1,1:C)=0.0
      Z(R,1:C)=0.0
      Z(1:R,1)=0.0
      Z(1:R,C)=0.0
      U0(1,1:C)=0.0
      U0(R,1:C)=0.0
      U0(1:R,1)=0.0
      U0(1:R,C)=0.0

C     Obtain SOR parameter W:
      print *, "SOR param: W >"
      read *, W
      print *, "epsilon: >"
      read *, epsx
      print *, "offset: >"
      read *,offset
      print *, "        ITS, EPS"
      it = 0
      notsure=0
      do while (notsure < 2)
      eps=1.0
       do while ( (abs(eps) > epsx))
        U0 = Sor(U0, nits, W, offset)    !Number of ITerationS=nits
         eps=maxU1-maxU0
        it = it + nits
        print *, it, eps
D        call showme(U0)
       end do
      notsure=notsure+1
      end do

      end program TestRedBlack

C===========================================================================
C===========================================================================
C mpfortran -pesize=8 -d_lines jacobi.f
C
C---ROB SALGADO (cps615g)
C   jacobi.f (3.f) 11/08/92 (extended from 11/06/92 PROBLEM 5.3)
C
C   1  4  1
C   4 -20 4
C   1  4  1
C
 
      subroutine showme(U)
      integer, parameter :: R=128,C=64
      real, array(R,C) :: U
      integer i,j
      print *,"------------------BOUNDARY POINTS + SAMPLE OF INTERIOR--"
      print 199, [1:C-1:int((C-1)/7)],C
      do, i=1, R-1, int((R-1)/11)
         print 200,i,U(i,1:C-1:int((C-1)/7)),U(i,C)
      end do
      print 200,R,U(R,1:C-1:int((C-1)/7)),U(R,C)
      print *,"-------------------------------------------------------/"
      print *," "
 199    FORMAT("    ",8I8)
 200      FORMAT(I4,8F8.4)
      end subroutine showme
 
C  Elementary solution of Laplace's equation in 2 dimensions
C  using the red/black SOR iteration.
C  Parameters R, C give array dimensions
 
 
C---AS USUAL, WE WON'T BOTH CHANGING ALL OF THE NAMES.
C
CC  Perform n cycles of SOR on U with relaxation parameter W
CC  Returns updated function, writes error estimate in mxr
CC  Error estimate is max abs change over all n iterations.
C
      function Sor(U, n, W, mxr)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: Sor, U
      integer n
      real W, mxr
 
      real, array(R, C) ::  U1, U2
      logical, array(R, C) :: Red
      integer k
      real OnemW                        ! 1 - W
 
C
C      OnemW = 1.0 - W
CC  Test for coord difference i - j being even
C      Red = Iand( SPREAD([1:R], 2, C) - SPREAD([1:C], 1, R),
C     *            1) == 0
C      U1 = U
C      U2 = U
C
 
C-----JACOBI FOR  1  4  1 ---------------------------------------------
C                 4 -20 4
C                 1  4  1
C
      U1=U
      do, k = 1, n
 
      U1(2:R-1, 2:C-1)=(4*( U1(2:R-1, 1:C-2)+U1(2:R-1, 3:C  )
     *                     +U1(1:R-2, 2:C-1)+U1(3:R,   2:C-1))
     *                   +( U1(1:R-2, 1:C-2)+U1(1:R-2, 3:C  )
     *                     +U1(3:R  , 1:C-2)+U1(3:R  , 3:C  )))/20.0
 
 
C       Where ( Red(2:R-1, 2:C-1) )
C        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
C     *                           + U1(3:R  , 2:C-1)
C     *                           + U1(2:R-1, 1:C-2)
C     *                           + U1(2:R-1, 3:C  ) )
C        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
C     *                    + OnemW*U1(2:R-1, 2:C-1)
C       Elsewhere
C        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
C     *                           + U1(3:R  , 2:C-1)
C     *                           + U1(2:R-1, 1:C-2)
C     *                           + U1(2:R-1, 3:C  ) )
C        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
C     *                    + OnemW*U1(2:R-1, 2:C-1)
C       end Where
 
      end do
 
      mxr = maxval(abs(U1 - U)) ! error estimate
      Sor = U1
      end function Sor
 
 
C  Test Redblack SOR
 
      program TestRedblack
 
      INTERFACE
        function Sor(U, n, W, mxr)
        integer, parameter :: R = 128, C = 64
        real, array(R, C) :: Sor, U
CMPF    ONDPU Sor,U
        integer n
        real W, mxr
      END INTERFACE
 
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: U0, Z
      integer, array(R, C) :: I, J
      real bsum, eps, err, W
      integer it
 
      I = SPREAD( [1:R], 2, C )         !Coordinates
      J = SPREAD( [1:C], 1, R )
 
C     Z is a solution of the difference equation
 
      Z = -1.4 + 0.12*I - 0.23*J + 0.001*(I**2 - J**2)
 
      U0 = Z    ! Install boundary values
      print *, "U0=Z"
      call showme(U0)
 
      bsum = SUM(U0(1,:)) + SUM(U0(R,:)) ! Sum on boundary
     *     + SUM(U0(:,1)) + SUM(U0(:,C))
 
C     Intialize interior to average on boundary
C     The "-4" prevents counting corners twice
 
      U0(2:R-1, 2:C-1) = bsum/(2*R + 2*C - 4)
      print *, "ORIGINAL, UNSOPHISTICATED"
      call showme(U0)
C     Obtain SOR parameter W:
      print *, " W >"
      read *, W
      print *, " it, eps, err"
      it = 0
      eps = 1.0
      do while ( eps > 0.01 )
        U0 = Sor(U0, 1000, W, eps)
        err = maxval(abs(U0 - Z))       ! Actual error
        it = it + 1000
        print *, it, eps, err
        call showme(U0)
      end do
 
      end program TestRedBlack

C===========================================================================
C===========================================================================
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*******************************************************
C===========================================================================
C===========================================================================
C  echo 'prob2.f is not compilable/executable source code see:  nb2.f'
C
C----ROB SALGADO (cps615g)
C prob2.f (rk.f) 10/12/92 (PROBLEM 4.2 -- SEE PROBLEM 5.1 FOR REVISION)
C
C RUNGE-KUTTA (Prob 2)
C
C ---------------------------------------IGNORABLE:  ORIGINAL COMMENTS
C                              
C
C I compiled the program with the mpf alias
C        mpf -nofecommon nb2.f
C
C This is the result of my attempt to run it.
C I have to admit, I don't know what $ Print ctl> is asking for.
C I'll try this again.
C But for now,  here is the output.
C The subroutine follows.
C 
C
C :/\/\/\/\ 
C
C magpie% a.out
C $ Data file>
C f7b16.dat
C $ Print ctl>
C 2
C FIO999 ENDFIL end of file.
C
C MPPE : intercepted signal SIGILL(4) @ FE pc 0x42b170
C on behalf of /u2/cps615g/a.out
C MPPE : You have the DPU
C  ID            Function symbol SRC  SIDE   STAT     FE PC |   ACU PC
C   0               NB2MODE(480):333, BOTH, UFUNC 0x0403f9c | 0x006684
C No MPPEs registered.  Try again? (Y/N)[N]: n
C
C :/\/\/\/ I DON'T UNDERSTAND WHAT THIS MEANS
C --------------------------------------------------------------------


C
C
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
      k1=h*V0                                    ! X'=V
      k2=h*(V0+k1/2.0)
      k3=h*(V0+k2/2.0)
      k4=h*(V0+k3)
      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)+k1/2.0)
      k3=h*(Grav(X0)+k2/2.0)
      k4=h*(Grav(X0)+k3)
      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----------------
C===========================================================================
C===========================================================================
C mpfortran -pesize=8 -d_lines mpfrand.o sph.f 
C
C
C-----ROB SALGADO (cps615g)
C     sph.f  10/22/92 (PROBLEM 4.1)
C
C
C     CALCULATE THE VOLUME IN A UNIT n-BALL
C
C     why does this method work?
C     The idea came to me as follows:
C         One can map a family of concentric squares within the unit square
C         into a family of concentric circles within the unit disk, and vice
C         versa.
C         Given a random point (xr,yr) in the unit square, map it to a point on
C         a circle of radius max(xr,yr).  Although for the purpose of this
C         calculation, it is not necessary to specify exactly where on the
C         circle, I chose to map it to a point at obtained by "projection"
C         from the line y=yr (assuming yr>xr) to the circle.
C         So, now I see that a thin shell of points around this circle makes
C         a contribution to the area of the disk as a function of the
C         radius of this determined circle.  Thus, I should weight the
C         points randomly picked in the unit square by a function of the
C         radius.  I looked for an inverse map...
C         Although I have not explicitly derived it analytically for all n,
C         when I apply it to the problem of calculating the enclosed volume
C         it seems to work reasonably well.
C
C     P.S.
C         The program seems to bomb when I use fixed matrices of size
C         greater than (16384,~20)...  I'm not sure why this is so.
C         But, the program is simple enough in structure to either
C         use several matrices or use "two rows or so" of random numbers
C         to make up the components of my position vectors.
C         Since the program seems to run relatively fast, I would like to
C         try to calculate the volume of, say, a 100-ball.
C
C
 
      program mppsphere
 
      integer, parameter :: N=16384, maxdim=20
      real, parameter :: Two31=2147483648.0
      integer IR(N*maxdim),  SR(N*maxdim)
      real v(N,maxdim)
 
      integer dim
      real, parameter :: pi=3.1415926535
 
CMPF  mpl MPF_RANDOM, MPF_SRANDOM
      integer mptm                     ! timer

      integer samp,r
      real logv,truevol
 
      real, array(N) :: t, l, hit, oldhit
      real, array(1000) :: monterad,oldmonterad
      real mean, oldmean, sumdevsq, oldsumdevsq
 
 
      open (unit=3, file='d-sphere')
      write(3,*) 'Importance d-Sphere Monte-Carlo [ROB SALGADO]'
      print*, 'Importance d-Sphere Monte-Carlo [ROB SALGADO]'
 
C-----SEED THE RANDOM NUMBER GENERATOR
      SR= [1:maxdim*N]
      call MPF_SRANDOM(SR)
 
      print*, 'Dimension dim=?'
      read*, dim
C-----------------------------COMPUTE TRUE D-VOLUME
C     Use logarithms because partial factors for high dimensions
C     overflows the variables.
      if (2*int(dim/2).eq.dim) goto 8
C-------DIM IS ODD
      logv=(dim+1.0)/2.0*log(2.0)+(dim-1.0)/2.0*log(pi)
      goto 9
C-------DIM IS EVEN
 8    logv=(dim/2.0)*log(2.0*pi)
C-------COMPUTE PARITY-FACTORIAL n(n-2)(n-4)...
 9    do 5 i=dim,1,-2
           logv=logv-log(i*1.0)
 5    continue
      truevol=exp(logv)
 
C-------------------------------------------------|
      print*, 'Number of samples per run?'
      read*, samp
 10   print*, 'Number of runs?'
      read*, r
      if (r.gt.1000) goto 10
 
C------------------------------------------------------(r) RUNS
      write(3,*) '               Trial     Radius'
D     print*,    '               Trial     Radius'
      do 200 j=1,r
         hit=0
         oldhit=0
C---------------------------------(N*samp) SAMPLES

      call mpTimerStart() 
      do 100 i=1,samp
         l=0.0
         t=0.0
 
      call MPF_RANDOM(IR)
      v=reshape( SHAPE=[N,maxdim], SOURCE=IR )/Two31
      v(:,(dim+1):maxdim)=0.0
                                                     ! shape into dim-vectors
      t=maxval(v,2)                                  ! get maximum coordinate
      l=sum(v*v,2)                                   ! get square-norm
 
C------weight the randompoint in the unit d-cube
C      by (t/r)**dim, where t=max(x1,x2,x3,..,xdim) for d-dimensions
C                     where r=sqrt(l)=sqrt(x1^2+x2^2+...+xdim^2)
      hit=hit+(sqrt(t*t/l))**dim
 
C------use the same random numbers to Monte-Carlo the naive way
C
      where ( l < 1.0) oldhit=oldhit+1
 
 100  continue
C------------------------------------------|
      monterad(j)=2**dim*sum(hit)/(N*samp)
      oldmonterad(j)=2**dim*sum(oldhit)/(N*samp)
      write(3,*) j, monterad(j),oldmonterad(j)
D     print*, j, monterad(j),oldmonterad(j)
      mptm=mpTimerElapsed()
D     print*, 'time (sec): ', mptm/1000.0
 200  continue
C-----------------------------------------------------------|
 
C------------------------------STATISTICS
      mean=sum(monterad)/r
      oldmean=sum(oldmonterad)/r
      write(3,*) '----'
      write(3,*) r,' Runs, each with',N*samp,' samples'
      write(3,*) 'Spherical volume in ',dim,'-dimensions'
      sumdevsq=sum((mean-monterad)*(mean-monterad))
      oldsumdevsq=sum((oldmean-oldmonterad)*(oldmean-oldmonterad))
 300  continue
 
 
      write(3,*) oldmean,'=NAIVE      ',sqrt(oldsumdevsq)/n,
     *           ' std dev (mean)'
      write(3,*) mean,'=IMPORTANCE ',sqrt(sumdevsq)/n,' std dev (mean)'
      write(3,*) truevol,'=TRUE volume'
      write(3,*) '% Naive Error     =', 100.0*(oldmean-truevol)/truevol
      write(3,*) '% Importance Error=', 100.0*(mean-truevol)/truevol
 
      end
C===========================================================================
C===========================================================================
C mpfortran -pesize=8 sorshift.f
C
C
C----ROB SALGADO
C    sorshift.f (2.f) 11/13/92 (PROBLEM 6.1)
C
C
      subroutine showme(U)
      integer, parameter :: R=128,C=64
      real, array(R,C) :: U
      integer i,j
      print *,"------------------BOUNDARY POINTS + SAMPLE OF INTERIOR--"
      print 199, [1:C-1:int((C-1)/7)],C
      do, i=1, R-1, int((R-1)/11)
         print 200,i,U(i,1:C-1:int((C-1)/7)),U(i,C)
      end do
      print 200,R,U(R,1:C-1:int((C-1)/7)),U(R,C)
      print *,"-------------------------------------------------------/"
      print *," "
 199  FORMAT("    ",8I8)
 200  FORMAT(I4,8F8.4)
      end subroutine showme

C  Elementary solution of Laplace's equation in 2 dimensions
C  using the red/black SOR iteration.
C  Parameters R, C give array dimensions

C  Perform n cycles of SOR on U with relaxation parameter W
C  Returns updated function, writes error estimate in mxr
C  Error estimate is max abs change over all n iterations.

      function Sor(U, n, W, offset)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: Sor, U
      integer n
      real W, offset
      real maxU1, maxU0
      COMMON maxU1, maxU0

      real, array(R, C) ::  U1, U2, UU
      logical, array(R, C) :: Red
      integer k
      real OnemW                        ! 1 - W

      OnemW = 1.0 - W
C  Test for coord difference i - j being even
      Red = Iand( SPREAD([1:R], 2, C) - SPREAD([1:C], 1, R),
     *            1) == 0
      U1 = U
      U2 = U
      do, k = 1, n
         UU=U
       Where ( Red(2:R-1, 2:C-1) )
        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
     *                           + U1(3:R  , 2:C-1)
     *                           + U1(2:R-1, 1:C-2)
     *                           + U1(2:R-1, 3:C  ) )
        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
     *                    + OnemW*U1(2:R-1, 2:C-1)
       Elsewhere
        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
     *                           + U1(3:R  , 2:C-1)
     *                           + U1(2:R-1, 1:C-2)
     *                           + U1(2:R-1, 3:C  ) )
        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
     *                    + OnemW*U1(2:R-1, 2:C-1)
       end Where
      U1=U1-offset*UU
       maxU0=maxU1
       maxU1 = maxval(abs(U1))  ! maximum absolute value obtained
C  Re-normalize to the maximum value
      U1= U1/maxU1
      end do
C  Re-normalize to the maximum value
      Sor = U1
      print *,"Current:                    ",maxU1
      end function Sor


            
C  Test Redblack SOR

      program TestRedblack

      INTERFACE
        function Sor(U, n, W, offset)
        integer, parameter :: R = 128, C = 64
        real, array(R, C) :: Sor, U
CMPF    ONDPU Sor,U
        integer n
        real W, offset
      END INTERFACE

      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: U0, Z
      integer, array(R, C) :: I, J
      real bsum, eps, err, W, epsx, offset, maxU0, maxU1
      integer it, nits, notsure
      COMMON maxU1, maxU0
C--Ask for number of iterations
C
      print *, "How many NITS? "
      read *,nits

      I = SPREAD( [1:R], 2, C )         !Coordinates
      J = SPREAD( [1:C], 1, R )

C     Z is a solution of the difference equation

      Z = -1.4 + 0.12*I - 0.23*J + 0.001*(I**2 - J**2)
      U0 = Z    ! Install boundary values
C      print  *, "ORIGINAL ANALYTIC SOLUTION... ZERO BCs IMPOSED LATER"
C      call showme(U0)

       bsum = SUM(U0(1,:)) + SUM(U0(R,:)) ! Sum on boundary
     *     + SUM(U0(:,1)) + SUM(U0(:,C))
C
CC     Intialize interior to average on boundary
CC     The "-4" prevents counting corners twice
C
       U0(2:R-1, 2:C-1) = bsum/(2*R + 2*C - 4)
C       print *,"INITIAL GUESS: UNSOPHISTICATED"
C       call showme(U0)
      Z(1,1:C)=0.0
      Z(R,1:C)=0.0
      Z(1:R,1)=0.0
      Z(1:R,C)=0.0
      U0(1,1:C)=0.0
      U0(R,1:C)=0.0
      U0(1:R,1)=0.0
      U0(1:R,C)=0.0

C     Obtain SOR parameter W:
      print *, "SOR param: W >"
      read *, W
      print *, "epsilon: >"
      read *, epsx
      print *, "offset: >"
      read *,offset
      print *, "        ITS, EPS"
      it = 0
      notsure=0
      do while (notsure < 2)
      eps=1.0
       do while ( (abs(eps) > epsx))
        U0 = Sor(U0, nits, W, offset)    !Number of ITerationS=nits
         eps=maxU1-maxU0
        it = it + nits
        print *, it, eps
D        call showme(U0)
       end do
      notsure=notsure+1
      end do

      end program TestRedBlack

C===========================================================================
C===========================================================================
