search.f

Go to the documentation of this file.

      SUBROUTINE search(I,IKS)
*
*
*       Close encounter search.
*       -----------------------
*
      include 'common6.h'
      common/clump/   bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
     &                names(ncmax,5),isys(5)
*
*
*       Increase counter for regularization attempts.
      nkstry = nkstry + 1
*
      fmax = 0.0
      nclose = 0
*       Find dominant neighbour by selecting all STEP(J) < 2*DTMIN.
      l = list(1,i) + 2
    2 l = l - 1
      IF (l.LT.2) go to 10
*
*       First see whether any c.m. with small step is within 2*RMIN.
      IF (list(l,i).LE.n) go to 4
*
      j = list(l,i)
      IF (step(j).GT.smin) go to 2
      a1 = x(1,j) - x(1,i)
      a2 = x(2,j) - x(2,i)
      a3 = x(3,j) - x(3,i)
      rij2 = a1*a1 + a2*a2 + a3*a3
      IF (rij2.GT.rmin22) go to 2
*
      fij = body(j)/rij2
      IF (fmax.LT.fij) fmax = fij
*       Abandon further search if c.m. force exceeds half the total force.
      IF (fmax**2.LT.f(1,i)**2 + f(2,i)**2 + f(3,i)**2) THEN
          nclose = nclose + 1
          jlist(nclose) = j
          go to 2
      ELSE
          go to 10
      END IF
*
*       Continue searching single particles with current value of FMAX.
    4 jcomp = 0
*
      DO 6 k = l,2,-1
          j = list(k,i)
          IF (step(j).GT.smin) go to 6
          a1 = x(1,j) - x(1,i)
          a2 = x(2,j) - x(2,i)
          a3 = x(3,j) - x(3,i)
          rij2 = a1*a1 + a2*a2 + a3*a3
          IF (rij2.LT.rmin22) THEN
              nclose = nclose + 1
              jlist(nclose) = j
*       Remember index of every single body with small step inside 2*RMIN.
              fij = (body(i) + body(j))/rij2
              IF (fij.GT.fmax) THEN
                  fmax = fij
*       Save square distance and global index of dominant body.
                  rjmin2 = rij2
                  jcomp = j
              END IF
          END IF
    6 CONTINUE
*
*       See whether dominant component is a single particle inside RMIN.
      IF (jcomp.LT.ifirst.OR.jcomp.GT.n) go to 10
*       Accept one single candidate inside 2*RMIN (which makes PERT = 0).
      IF (rjmin2.GT.rmin2.AND.nclose.GT.1) go to 10
*
      rdot = (x(1,i) - x(1,jcomp))*(xdot(1,i) - xdot(1,jcomp)) +
     &       (x(2,i) - x(2,jcomp))*(xdot(2,i) - xdot(2,jcomp)) +
     &       (x(3,i) - x(3,jcomp))*(xdot(3,i) - xdot(3,jcomp))
*
*       Only select approaching particles (include nearly circular case).
      rijmin = sqrt(rjmin2)
      IF (rdot.GT.0.02*sqrt((body(i) + body(jcomp))*rijmin)) go to 10
*
*       Ensure a massive neighbour is included in perturbation estimate.
      bcm = body(i) + body(jcomp)
      IF (body1.GT.10.0*bcm) THEN
          jbig = 0
          big = bcm
          nnb1 = list(1,i) + 1
          DO 20 l = 2,nnb1
              j = list(l,i)
              IF (body(j).GT.big) THEN
                  jbig = j
                  big = body(j)
              END IF
   20     CONTINUE
*       Check whether already present, otherwise add to JLIST.
          DO 25 l = 1,nclose
              IF (jlist(l).EQ.jbig) THEN
                  jbig = 0
              END IF
   25     CONTINUE
          IF (jbig.GT.0) THEN
              nclose = nclose + 1
              jlist(nclose) = jbig
          END IF
      END IF
*
*       Evaluate vectorial perturbation due to the close bodies.
      CALL fpert(i,jcomp,nclose,pert)
*
*       Accept #I & JCOMP if the relative motion is dominant (GI < 0.25).
      gi = pert*rjmin2/bcm
      IF (gi.GT.0.25) THEN
*         IF (KZ(4).GT.0.AND.TIME-TLASTT.GT.4.44*TCR/FLOAT(N))
*    &                                             CALL EVOLVE(JCOMP,0)
          go to 10
      END IF
*
*       Exclude any c.m. body of compact subsystem (TRIPLE, QUAD or CHAIN).
      DO 8 isub = 1,nsub
          namei = names(1,isub)
          IF (namei.EQ.name(i).OR.namei.EQ.name(jcomp)) go to 10
    8 CONTINUE
*
*       Also check possible c.m. body of chain regularization (NAME = 0).
      IF (nch.GT.0) THEN
          IF (name(i).EQ.0.OR.name(jcomp).EQ.0) go to 10
      END IF
*
*       Save index and increase indicator to denote new regularization.
      icomp = i
      iks = iks + 1
*
   10 RETURN
*
      END