shrink.f

Go to the documentation of this file.

      SUBROUTINE shrink(TMIN)
*
*
*       Shrinking of large time-steps.
*       ------------------------------
*
      include 'common6.h'
      parameter(dtm = 0.03125d0)
      REAL*8  rij(3),vij(3)
*
*
*       Determine impact to each high-velocity particle for large STEPR.
      nnb = listv(1)
      DO 20 l = 1,nnb
          i = listv(l+1)
          DO 10 j = ifirst,ntot
*       Form relative quantities for large regular time-steps.
              IF (stepr(j).LT.dtm.OR.j.EQ.i) go to 10
*       Check whether body #I is already a neighbour (accept LIST > I).
              nnb1 = list(1,j) + 1
              DO 1 k = 2,nnb1
                  IF (list(k,j).GT.i) go to 2
                  IF (list(k,j).EQ.i) go to 10
    1         CONTINUE
    2         rv = 0.0
              vv = 0.0
              DO 3 k = 1,3
                  rij(k) = x(k,i) - x(k,j)
                  vij(k) = xdot(k,i) - xdot(k,j)
                  rv = rv + rij(k)*vij(k)
                  vv = vv + vij(k)**2
    3         CONTINUE
*
*       Skip treatment for receding particles.
              IF (rv.GE.0.0d0) go to 10
*
*       Evaluate time of minimum approach and truncate to next time.
              dt = -rv/vv
              it = 0
    5         dt = min(t0r(j) + stepr(j) - time,dt)
*
*       Obtain minimum impact parameter for straight-line orbit.
              r2 = 0.0
              fj2 = 0.0
              DO 6 k = 1,3
                  r2 = r2 + (rij(k) + vij(k)*dt)**2
                  fj2 = fj2 + f(k,j)**2
    6         CONTINUE
*
*       Compare force at minimum distance with total force on body #J (> 0).
              fi2 = (body(i)/r2)**2
              IF (fi2.LT.0.04*fj2.OR.body(j).EQ.0.0d0) go to 10
*
*       See whether the regular time-step can be shortened.
              IF (t0r(j) + 0.5*stepr(j).GE.tmin.AND.it.LT.5) THEN
                  WRITE (29,8)  j, sqrt(r2), sqrt(fi2/fj2), dt, stepr(j)
    8             FORMAT (' SHRINK    J R FI/FJ DT STEP ',
     &                                i6,f7.3,f6.2,1p,2e10.2)
                  CALL flush(29)
                  stepr(j) = 0.5*stepr(j)
                  IF (step(j).GT.stepr(j)) THEN
                      step(j) = 0.5*step(j)
                      tnew(j) = t0(j) + step(j)
                  END IF
                  it = it + 1
                  go to 5
              END IF
   10     CONTINUE
   20 CONTINUE
*
      RETURN
*
      END