xvpred.f

Go to the documentation of this file.

      SUBROUTINE xvpred(I1,NNB)
*
*
*       Prediction of coordinates & velocities.
*       ---------------------------------------
*
      include 'common6.h'
*
*
*       Perform full N prediction to highest order (defined by NNB >= N).
      IF (nnb.GE.n) go to 20
*       Use high order for special singles (NNB < 0) or c.m. bodies (I1 > N).
      IF (nnb.LT.0.OR.(nnb.EQ.0.AND.i1.GT.n)) go to 20
*
*       Adopt low order for standard single particles (I1 <= N & NNB = 0).
      IF (nnb.EQ.0) THEN
          i = i1
          s = time - t0(i)
          IF (s.EQ.0.0d0) go to 40
          s1 = 1.5*s
          s2 = 2.0*s
          x(1,i) = ((fdot(1,i)*s + f(1,i))*s + x0dot(1,i))*s + x0(1,i)
          x(2,i) = ((fdot(2,i)*s + f(2,i))*s + x0dot(2,i))*s + x0(2,i)
          x(3,i) = ((fdot(3,i)*s + f(3,i))*s + x0dot(3,i))*s + x0(3,i)
          xdot(1,i) = (fdot(1,i)*s1 + f(1,i))*s2 + x0dot(1,i)
          xdot(2,i) = (fdot(2,i)*s1 + f(2,i))*s2 + x0dot(2,i)
          xdot(3,i) = (fdot(3,i)*s1 + f(3,i))*s2 + x0dot(3,i)
          go to 40
      END IF
*
*       Predict all neighbours to order FDOT.
      nnb1 = nnb + 1
      DO 5 l = 2,nnb1
          i = list(l,i1)
          s = time - t0(i)
          IF (s.EQ.0.0d0) go to 5
          s1 = 1.5*s
          s2 = 2.0*s
          x(1,i) = ((fdot(1,i)*s + f(1,i))*s + x0dot(1,i))*s + x0(1,i)
          x(2,i) = ((fdot(2,i)*s + f(2,i))*s + x0dot(2,i))*s + x0(2,i)
          x(3,i) = ((fdot(3,i)*s + f(3,i))*s + x0dot(3,i))*s + x0(3,i)
          xdot(1,i) = (fdot(1,i)*s1 + f(1,i))*s2 + x0dot(1,i)
          xdot(2,i) = (fdot(2,i)*s1 + f(2,i))*s2 + x0dot(2,i)
          xdot(3,i) = (fdot(3,i)*s1 + f(3,i))*s2 + x0dot(3,i)
    5 CONTINUE
*
*       Resolve the components of any perturbed pairs (last neighbours).
   10 i = list(nnb1,i1)
      IF (i.GT.n) THEN
          jpair = i - n
          IF (list(1,2*jpair-1).GT.0) THEN
              CALL resolv(jpair,1)
          END IF
          nnb1 = nnb1 - 1
          IF (nnb1.GT.1) go to 10
      END IF
      go to 40
*
*       Set index of first body (#I1 if NNB = 0 or NNB >= N).
   20 i = i1
*       Predict coordinates & velocities of body #I to order F3DOT.
   25 dt = time - t0(i)
      IF ((dt.EQ.0.0d0.AND.iphase.LT.4).OR.body(i).EQ.0.0d0) go to 35
      a1 = 0.05*dt
      a2 = 0.25*dt
      a3 = t0(i) - t0r(i)
      a4 = 0.5*a3
*
      DO 30 k = 1,3
          fk = ((one6*d3r(k,i)*a3 + 0.5*d2r(k,i))*a3 + d1r(k,i))*a3
     &                                              + d0r(k,i) + d0(k,i)
          f1dotk = (d3r(k,i)*a4 + d2r(k,i))*a3 + d1r(k,i) + d1(k,i)
          f2dotk = 0.5*(d3r(k,i)*a3 + d2r(k,i) + d2(k,i))
          f3dotk = one6*(d3r(k,i) + d3(k,i))
          x(k,i) = ((((f3dotk*a1 + one12*f2dotk)*dt +
     &                   one6*f1dotk)*dt + 0.5*fk)*dt + x0dot(k,i))*dt +
     &                   x0(k,i)
          xdot(k,i)  = (((f3dotk*a2 + one3*f2dotk)*dt +
     &                            0.5*f1dotk)*dt + fk)*dt + x0dot(k,i)
   30 CONTINUE
*
*       Resolve the components of perturbed pairs.
   35 IF (i.GT.n) THEN
          jpair = i - n
          IF (list(1,2*jpair-1).GT.0) THEN
              CALL resolv(jpair,1)
          END IF
      END IF
*
*       Check whether prediction of single particle or full N.
      IF (nnb.GT.0) THEN
          i = i + 1
          IF (i.LE.nnb) go to 25
      END IF
*
   40 RETURN
*
      END