insert.f

Go to the documentation of this file.

      SUBROUTINE insert(I,LI)
*
*
*       Insert particle index in KS time-step list.
*       -------------------------------------------
*
      include 'common6.h'
      REAL*8  ti,ti1,ti2,tj
*
*
*       Avoid increasing KBLIST if body #I can be exchanged with next member.
      li2 = min(li + 2,nntb)
      i2 = kblist(li2)
      ti = t0(i) + step(i)
      ti2 = t0(i2) + step(i2)
*
*       Check swapping condition (#I due before #I2 but after #I1 at LI + 1).
      IF (ti.LE.ti2) THEN
          li1 = min(li+1,nntb)
          i1 = kblist(li1)
          ti1 = t0(i1) + step(i1)
          IF (ti.GT.ti1) THEN
              kblist(li) = kblist(li1)
              kblist(li1) = i
          END IF
*       Reduce pointer so that current location will be selected again.
          li = li - 1
          go to 50
*       Also swap if body #I is among the two last members and TI > TI2.
      ELSE IF (li.GT.nntb - 2) THEN
          kblist(li) = kblist(li2)
          kblist(li2) = i
          li = li - 1
          go to 50
      END IF
*
*       Estimate the insert index from the remaining time interval.
      fac = step(i)/(tblist - time)
      lstar = li + float(nntb - li)*fac
*
*       Improve insert index by another iteration (check LI < LSTAR <= NNTB).
      j = kblist(lstar)
      tj = t0(j) + step(j)
*       Avoid division by zero (unperturbed steps may also be quantized).
      IF (tblist.NE.tj) THEN
          fac = (ti - tj)/(tblist - tj)
          lstar = lstar + float(nntb - lstar)*fac
          lstar = max(lstar,li+1)
          lstar = min(lstar,nntb)
          j = kblist(lstar)
      END IF
*
*       Determine correct index by comparing neighbouring KBLIST members.
      IF (ti.GE.t0(j) + step(j)) THEN
          l1 = lstar + 1
          lstar = l1
          DO 10 l = l1,nntb
              j = kblist(l)
              tj = t0(j) + step(j)
*       Advance index until TI < TJ or last member.
              IF (ti.LT.tj) go to 30 
              lstar = l + 1
   10     CONTINUE
      ELSE
   20     lstar = lstar - 1
          j = kblist(lstar)
          IF (j.EQ.i) THEN
              li = li - 1
              go to 50
          END IF
          IF (ti.LT.t0(j) + step(j)) go to 20
          lstar = lstar + 1
*       Ensure that current index LI is not chosen due to round-off.
          lstar = max(lstar,li+1)
      END IF 
*
*       Create free location at LSTAR by moving all subsequent members down.
   30 DO 40 l = nntb,lstar,-1
          kblist(l+1) = kblist(l)
   40 CONTINUE 
*
*       Insert body #I and update memberships.
      kblist(lstar) = i
      nntb = nntb + 1
*
   50 RETURN
*
      END