C
C     Simple Multiparticle Interaction Simulator
C
C     This program calculates the interactions of N particles, all 
C     interacting by an inverse-square force law.  
C
C
      REAL*4 POS(3, 20), VEL(3, 20), MASS(20)
      REAL*4 TIME/0.0/, FINISH/1.0/, DT/0.1/
      INTEGER*4 N
C
C     N is the number of particles; the program may take up to 20
C     before changing of the code will be necessary.
C
C     POS records the x, y, and z coordinates of each of N particles;
C     VEL the components of the velocity of each.  Acceleration is 
C     calculated within a subroutine and not globally stored.
C
C     MASS is a constant for each particle representing the strength 
C     of its force.  It is called mass to suggest gravitational 
C     forces; to simulate electrical forces MASS may take negative 
C     numbers.
C
C     TIME is the current system time, and begins initially at 0.0.
C     FINISH is the time at which the computation is to stop, and it 
C     defaults to 10.  DT is the time passing between each integration 
C     step, and it defaults to 0.1.
C

      N = 3
      CALL INCOND(POS, VEL, MASS, N)
      CALL PTINST(POS, VEL, MASS, N)
C
C
C     Below is the main integration loop.
C
C
 6510 CALL ADVNCE(POS, VEL, MASS, N, DT)
      CALL PTSTAT(POS, VEL, TIME, N)
      TIME = TIME + DT
      IF (TIME .LE. FINISH) GO TO 6510

      STOP
      END
C
C
C     Below are subroutines called by the program.
C
C
      SUBROUTINE PTSTAT(POS, VEL, TIME, N)
      REAL*4 POS(3, 20), VEL(3, 20), TIME
      INTEGER*4 N

      INTEGER*4 I
      WRITE (*, *)
      WRITE (*, 901) TIME
      DO 101, I = 1, N
         WRITE (*, 902) I, POS(1, I), POS(2, I), POS(3, I)
         WRITE (*, 903) VEL(1, I), VEL(2, I), VEL(3, I)
 101  CONTINUE
      
 901  FORMAT( ' ----- Time: ', F6.2,' ----- '  )
 902  FORMAT(I2, 1X, 'Position: ', 1X, 3(F12.3, 1X))
 903  FORMAT(4X, 'Velocity: ', 1X, 3(F12.3, 1X))
C
C     I2 means a two-character slot is allocated for an integer (in 
C     this case, the particle number).  1X allocates one blank space.
C     (Similarly, 3X allocates three blank spaces.)  'Position: ' 
C     and 'Velocity: ' appear as written.
C
C     3(F8.3, 1X) leaves space for three values, each placed into 
C     an eight-character wide slot, with three characters placed 
C     after the decimal.  One blank space follows each number.
C
C     Similarly, the F6.2 in the first FORMAT statement allocates a 
C     six-character block for the time, and reserves two spots for 
C     digits after the decimal.
C
      RETURN
      END


      SUBROUTINE PTINST(POS, VEL, MASS, N)
      REAL*4 POS(3, 20), VEL(3, 20), MASS(20)
      INTEGER*4 N
      
      INTEGER*4 I
      DO 111, I = 1, N
         WRITE (*, 905) I, POS(1, I), POS(2, I), POS(3, I)
         WRITE (*, 906) VEL(1, I), VEL(2,I), VEL(3, I)
         WRITE (*, 907) MASS(I)
 111  CONTINUE

 905  FORMAT(I2, 1X, 'Position: ', 1X, 3(F12.3, 1X))
 906  FORMAT(3X, 'Velocity: ', 1X, 3(F12.3, 1X))
 907  FORMAT(4X, 'Mass: ', 1X, F12.3)
C
C     Note the FORMAT commands are similar to those in the above 
C     PRINTSTATUS subroutine.
C
C
      RETURN
      END


      SUBROUTINE INCOND(POS, VEL, MASS, N)
      REAL*4 POS(3, 20), VEL(3, 20), MASS(20)
      INTEGER*4 N

      INTEGER*4 I, J
      REAL*4 COMPTS(6*20)

      CALL RNFSTR(109707)
      CALL RNFARR(MASS, N)
      DO 121, I = 1, N
         MASS(I) = AMOD(MASS(I), 1.0)
 121  CONTINUE
      CALL RNFSTR(109708)
      CALL RNFARR(COMPTS, 6*N)
      DO 122, I = 1, N
         DO 123, J = 1, 3
            POS(J, I) = AMOD(COMPTS(6*(I-1)+J), 1.0)
            VEL(J, I) = AMOD(COMPTS(6*(I-1)+J+3), 1.0)
 123     CONTINUE
 122  CONTINUE
      WRITE (*,*) MASS(1)
      RETURN
      END


      SUBROUTINE ADVNCE(POS, VEL, MASS, N, DT)
      REAL*4 POS(3, 20), VEL(3, 20), MASS(3, 20), DT
      INTEGER*4 N

      INTEGER*4 I
      REAL*4 ACC(3, 20)
      DO 131, I = 1, N
         CALL ACCEL(ACC, POS, VEL, MASS, N)
 131  CONTINUE
      CALL ADDVEC(VEL, VEL, ACC, DT, N)
      CALL ADDVEC(POS, POS, VEL, DT, N)
      RETURN
      END
C
C
C     The below subroutines are called by the ADVANCE subroutine.
C
C
      SUBROUTINE ACCEL(ACC, POS, VEL, MASS, N)
      REAL*4 ACC(3, 20), POS(3, 20), VEL(3, 20), MASS(20)
      INTEGER*4 N

      INTEGER*4 I, J, K
      REAL*4 DIST, FORCE

      DO 141, I = 1, N
         DO 142, K = 1, 3
            ACC(K, I) = 0.0
C            ACC(K, I) = -POS(K, I)
C     
C     For testing purposes: Remove the commenting on the above line, 
C     and add commenting at the line just before label 144, so that 
C     the system becomes simple harmonic oscillators in each of the 
C     coordinates of each of the particles.  
C
C
 142     CONTINUE
         DO 141, J = 1, N
            IF (I .EQ. J) GO TO 141
            DIST = 0.0
            DO 143, K = 1,3
               DIST = DIST + (POS(K, J) - POS(K, I))**2
 143        CONTINUE

            FORCE = MASS(J) / DIST
            DO 144, K = 1, 3
C
C     The below line may be commented out to test the system as an 
C     integrator for the simple harmonic oscillator, as mentioned 
C     in the comment block just before label 142 above.
C
C
               ACC(K, I) = FORCE * (POS(K, J) - POS(K, I))/SQRT(DIST)
 144        CONTINUE
 141     CONTINUE
      RETURN
      END
C
C
C     The ADDVECTORS routine is intended to add one vector to 
C     another.  As its use is for calculating v = v + a*dt, and 
C     x = x + v*dt, the multiplication by dt is folded into the 
C     procedure.  The multiplication by SCALAR (DT, in the uses 
C     this subroutine has) would need to be removed for most 
C     other programs.
C
C     
      SUBROUTINE ADDVEC(TOTAL, FIRST, SECOND, SCALAR, N)
      REAL*4 TOTAL(3, 20), FIRST(3, 20), SECOND(3, 20), SCALAR
      INTEGER*4 N

      INTEGER*4 I
      DO 151, I = 1, N
         DO 151, J = 1, 3
            TOTAL(J, I) = FIRST(J, I) + SECOND(J, I) * SCALAR
 151  CONTINUE
      RETURN
      END

         
C
C
C
C     Below is the RANARRAY program from Donald E. Knuth.  Its 
C     numbers are used for the initial conditions.
C
C

      SUBROUTINE RNFARR(AA,N)
C       FORTRAN 77 version of "ranf_array"
C       from Seminumerical Algorithms by D E Knuth, 3rd edition (1997)
C       ********* see the book for explanations and caveats! *********
C     Copyright (c) 2000 by the author;
C     copy it if you like, but don't change it!
      IMPLICIT DOUBLE PRECISION (A,R,X,Y)
      DIMENSION AA(*)
      PARAMETER (KK=100)
      PARAMETER (LL=37)
      COMMON /RSTATE/ RANX(KK)
      SAVE /RSTATE/
      DO 1 J=1,KK
 1    AA(J)=RANX(J)
      DO 2 J=KK+1,N
         Y=AA(J-KK)+AA(J-LL)
         AA(J)=Y-IDINT(Y)
 2    CONTINUE
      DO 3 J=1,LL
         Y=AA(N+J-KK)+AA(N+J-LL)
         RANX(J)=Y-IDINT(Y)
 3    CONTINUE
      DO 4 J=LL+1,KK
         Y=AA(N+J-KK)+RANX(J-LL)
         RANX(J)=Y-IDINT(Y)
 4    CONTINUE
      END

      SUBROUTINE RNFSTR(SEED)
      IMPLICIT DOUBLE PRECISION (A,R,U-Z)
      IMPLICIT INTEGER (T)
      PARAMETER (KK=100)
      PARAMETER (LL=37)
      PARAMETER (MM=2**30)
      PARAMETER (ULP=1D0/(2D0**52))
      PARAMETER (TT=70)
      PARAMETER (KKK=KK+KK-1)
      INTEGER SEED,S,SSEED
      DOUBLE PRECISION SS
      DIMENSION X(KKK), XL(KKK)
      COMMON /RSTATE/ RANX(KK)
      SAVE /RSTATE/
      IF (SEED .LT. 0) THEN
         SSEED=MM-1-MOD(-1-SEED,MM)
      ELSE
         SSEED=MOD(SEED,MM)
      END IF
      SS=2D0*ULP*(SSEED+2)
      DO 1 J=1,KK
         X(J)=SS
         XL(J)=0
         SS=SS+SS
         IF (SS .GE. 1D0) SS=SS-1D0+2*ULP
 1    CONTINUE
      X(2)=X(2)+ULP
      XL(2)=ULP
      DO 2 J=KK+1,KKK
      X(J)=0
 2    XL(J)=0
      S=SSEED
      T=TT-1
 10   DO 12 J=KK,2,-1
         XL(J+J-1)=XL(J)
 12      X(J+J-1)=X(J)
      DO 13 J=KKK,KK-LL+1,-2
         XL(KKK-J+2)=0
 13      X(KKK-J+2)=X(J)-XL(J)
      DO 14 J=KKK,KK+1,-1
         IF (XL(J) .NE. 0) THEN
            XL(J-(KK-LL))=ULP-XL(J-(KK-LL))
            Y=X(J-(KK-LL))+X(J)
            X(J-(KK-LL))=Y-IDINT(Y)
            XL(J-KK)=ULP-XL(J-KK)
            Y=X(J-KK)+X(J)
            X(J-KK)=Y-IDINT(Y)
         END IF
 14   CONTINUE
C     this part is redundant but useful while I'm debugging
      DO 15 J=1,KKK
      IF (DMOD(X(J),2*ULP) .NE. XL(J)) THEN
         PRINT '(I15)', J
      END IF
 15   CONTINUE
      IF (MOD(S,2) .EQ. 1) THEN
         DO 16 J=KK,1,-1
            XL(J+1)=XL(J)
 16         X(J+1)=X(J)
         X(1)=X(KK+1)
         XL(1)=XL(KK+1)
         IF (XL(KK+1) .NE. 0) THEN
            XL(LL+1)=ULP-XL(LL+1)
            Y=X(LL+1)+X(KK+1)
            X(LL+1)=Y-IDINT(Y)
         END IF
      END IF
      IF (S .NE. 0) THEN
         S=S/2
      ELSE
         T=T-1
      END IF
      IF (T .GT. 0) GO TO 10
      DO 20 J=1,LL
 20   RANX(J+KK-LL)=X(J)
      DO 21 J=LL+1,KK
 21   RANX(J-LL)=X(J)
      END