ntint.f

Go to the documentation of this file.

      SUBROUTINE ntint(I)
*
*
*       Single star integration.
*       ------------------------
*
      include 'common6.h'
      REAL*8  xi(3),xidot(3),firr(3),fd(3)
*
*
*       Predict position and velocity up to end of STEP.
      iter = 0
    1 j = i
      s = step(j)
      s1 = 1.5*s
      s2 = 2.0*s
      x(1,j) = ((fdot(1,j)*s + f(1,j))*s +x0dot(1,j))*s + x0(1,j)
      x(2,j) = ((fdot(2,j)*s + f(2,j))*s +x0dot(2,j))*s + x0(2,j)
      x(3,j) = ((fdot(3,j)*s + f(3,j))*s +x0dot(3,j))*s + x0(3,j)
      xdot(1,j) = (fdot(1,j)*s1 + f(1,j))*s2 + x0dot(1,j)
      xdot(2,j) = (fdot(2,j)*s1 + f(2,j))*s2 + x0dot(2,j)
      xdot(3,j) = (fdot(3,j)*s1 + f(3,j))*s2 + x0dot(3,j)
*
*       Copy X & XDOT to scalars and initialize FIRR & FD.
      DO 5 k = 1,3
          xi(k) = x(k,i)
          xidot(k) = xdot(k,i)
          firr(k) = 0.0d0
          fd(k) = 0.0d0
    5 CONTINUE
*
*       Evaluate the galactic force and first derivative.
      CALL xtrnlt(xi,xidot,firr,fd)
*
*       Include the corrector and set new T0, F, FDOT, D1, D2 & D3.
      dt = step(i)
      dtsq = dt**2
      dt6 = 6.0d0/(dt*dtsq)
      dt2 = 2.0d0/dtsq
      dtsq12 = one12*dtsq
      dt13 = one3*dt
      t0(i) = t0(i) + step(i)
*
      DO 10 k = 1,3
          df = fi(k,i) - firr(k)
          fid = fidot(k,i)
          sum = fid + fd(k)
          at3 = 2.0d0*df + dt*sum
          bt2 = -3.0d0*df - dt*(sum + fid)
*
          x0(k,i) = xi(k) + (0.6d0*at3 + bt2)*dtsq12
          x0dot(k,i) = xidot(k) + (0.75d0*at3 + bt2)*dt13
*
*       Update the corrected values (OK for test particles).
          x(k,i) = x0(k,i)
          xdot(k,i) = x0dot(k,i)
*
          fi(k,i) = firr(k)
          fidot(k,i) = fd(k)
          f(k,i) = 0.5d0*firr(k)
          fdot(k,i) = one6*fd(k)
*
*       Form derivatives even though not needed for commensurate times.
          d1(k,i) = fd(k)
          d2(k,i) = (3.0d0*at3 + bt2)*dt2
          d3(k,i) = at3*dt6
*       NOTE: These are real derivatives!
   10 CONTINUE
*
*       Specify new time-step by standard criterion.
      ttmp = tstep(firr,fd,d2(1,i),d3(1,i),etai)
      dt0 = ttmp
      ttime = t0(i)
*
*       Select discrete value (increased by 2, decreased by 2 or unchanged).
      IF (ttmp.GT.2.0*step(i)) THEN
          ttmp = min(2.0*step(i),1.0d0)
      ELSE IF (ttmp.LT.step(i)) THEN
          ttmp = 0.5*step(i)
      ELSE
          ttmp = step(i)
      END IF
*
*       Do not permit current TIME to be exceeded.
      IF (t0(i) + ttmp.GT.time) THEN
          ttmp = time - t0(i)
      END IF
*
*       Set new block step and update next time.
      step(i) = ttmp
      tnew(i) = step(i) + t0(i)
*
*       Increase step counter.
      nstail = nstail + 1
*
*       See whether to continue until end of large block-step.
      IF (tnew(i).LT.time) THEN
          iter = iter + 1
          IF (iter.LT.10) go to 1
          WRITE (6,20)  i, time, step(i), firr
   20     FORMAT (' SMALL TIDAL STEP    I T DT F ',i7,f7.1,1p,4e10.2)
      END IF
*
      RETURN
*
      END