https://github.com/florentrenaud/nbody6tt
Raw File
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
ntint.f
      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
back to top