flyby.f

Go to the documentation of this file.

      SUBROUTINE flyby(I,ITERM)
*
*
*      Flyby check of KS orbit. 
*       -----------------------
*
      include 'common6.h'
*
*
*       Set index of KS pair & first component of c.m. body #I.
      ipair = i - n
      i1 = 2*ipair - 1
      rjmin2 = 1000.0
      iterm = 0
      jcomp = list(2,i1)
*
*       Find the closest body (JCOMP) and copy perturbers to JLIST.
      nnb = list(1,i1)
      DO 10 l = 2,nnb+1
          j = list(l,i1)
          jlist(l-1) = j
          rij2 = (x(1,i) - x(1,j))**2 + (x(2,i) - x(2,j))**2 +
     &                                  (x(3,i) - x(3,j))**2
          IF (rij2.LT.rjmin2) THEN 
              rjmin2 = rij2
              jcomp = j
          END IF
   10 CONTINUE
*
*       Skip rare case of merged binary or chain c.m. (denoted by NAME <= 0).
      IF (name(i).LE.0.OR.name(jcomp).LE.0) go to 20
*
      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))
*
      vrel2 = (xdot(1,i) - xdot(1,jcomp))**2 + 
     &        (xdot(2,i) - xdot(2,jcomp))**2 +
     &        (xdot(3,i) - xdot(3,jcomp))**2
 
*       Evaluate outer semi-major axis, eccentricity & impact parameter.
      rij = sqrt(rjmin2)
      semi1 = 2.0/rij - vrel2/(body(i) + body(jcomp))
      semi1 = 1.0/semi1
      ecc1 = sqrt((1.0d0 - rij/semi1)**2 +
     &                     rdot**2/(semi1*(body(i) + body(jcomp))))
      pmin = semi1*(1.0d0 - ecc1)
*
*       Set semi-major axis & eccentricity of inner binary.
      semi = -0.5d0*body(i)/h(ipair)
      ecc2 = (1.0d0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
      ecc = sqrt(ecc2)
*
      i2 = i1 + 1
      rji1 = (x(1,jcomp) - x(1,i1))**2 + (x(2,jcomp) - x(2,i1))**2 +
     &                                   (x(3,jcomp) - x(3,i1))**2
      rji2 = (x(1,jcomp) - x(1,i2))**2 + (x(2,jcomp) - x(2,i2))**2 +
     &                                   (x(3,jcomp) - x(3,i2))**2
*
*       Identify the dominant interaction (#I1 or #I2).
      fji1 = (body(i1) + body(jcomp))/rji1
      fji2 = (body(i2) + body(jcomp))/rji2
      nnb = nnb + 1
      IF (fji1.GT.fji2) THEN
          j = i1
          rj2 = rji1
*       Add the other component to the perturber list.
          jlist(nnb) = i2
      ELSE
          j = i2
          rj2 = rji2
          jlist(nnb) = i1
      END IF
*
*       Evaluate radial velocity (routine SEARCH requires RD < 0).
      rd = 0.0
      rjj = 0.0
      vjj = 0.0
      DO 15 k = 1,3
          rd = rd + (x(k,jcomp) - x(k,j))*(xdot(k,jcomp) - xdot(k,j))
          rjj = rjj + (x(k,jcomp) - x(k,j))**2
          vjj = vjj + (xdot(k,jcomp) - xdot(k,j))**2
   15 CONTINUE
*
*       Determine vectorial perturbation on intruder and closest component.
      CALL fpert(jcomp,j,nnb,pert)
      gj = pert*rj2/(body(jcomp) + body(j))
      rj = sqrt(rjj)
*     WRITE (6,18) NAME(J),NAME(JCOMP),GAMMA(IPAIR),GJ,R(IPAIR),RJ,RD,
*    &             RDOT
*  18 FORMAT (' NMJ NMJC GI GJ R RJ RD RDI ',2I5,2F6.2,1P,4E9.1)
*
*       Terminate if perturber & I1/I2 provides dominant force.
      apo = abs(semi)*(1.0 + ecc)
      rsum = r(ipair) + sqrt(rj2)
      xf = 2.0
*       Compare current radial velocity with #JCOMP & J.
      IF (tdot2(ipair)/r(ipair).LT.rd/rj) xf = 1.0
*
*       Increase the cross section for quadruples.
      IF (jcomp.GT.n) THEN
          semi2 = -0.5*body(jcomp)/h(jcomp-n)
          apo = apo + abs(semi2)
      END IF
*
      semij = 2.0/rj - vjj/(body(jcomp) + body(j))
      semij = 1.0/semij
*       Check for chain regularization test or standard termination.
      IF (pmin.LT.apo.AND.rsum.LT.xf*rmin.AND.rdot.LT.0.0) THEN
          IF (semi.LT.0.0.AND.gamma(ipair).LT.0.5) THEN
              iterm = 0
              go to 20
          END IF
          iterm = 1
      ELSE IF (gj.LT.0.7*gamma(ipair).AND.rd.LE.0.0) THEN
          iterm = 2
      END IF
*
   20 RETURN
*
      END