Revision dd8daa98413988b11c594e5a3bff8f5659d7ecf7 authored by Florent Renaud on 26 June 2014, 09:41:19 UTC, committed by Florent Renaud on 26 June 2014, 09:41:19 UTC
Originally released on 18 March 2013
Based on Nbody6 downloaded on 29 January 2013
1 parent edab1b7
Raw File
flyby.f
      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
      JCL = LIST(2,I1)
*
*       Find the closest body (JCL) 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
              JCL = 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(JCL).LE.0) GO TO 20
*
      RDOT = (X(1,I) - X(1,JCL))*(XDOT(1,I) - XDOT(1,JCL)) +
     &       (X(2,I) - X(2,JCL))*(XDOT(2,I) - XDOT(2,JCL)) +
     &       (X(3,I) - X(3,JCL))*(XDOT(3,I) - XDOT(3,JCL))
*
      VREL2 = (XDOT(1,I) - XDOT(1,JCL))**2 + 
     &        (XDOT(2,I) - XDOT(2,JCL))**2 +
     &        (XDOT(3,I) - XDOT(3,JCL))**2
 
*       Evaluate outer semi-major axis, eccentricity & impact parameter.
      RIJ = SQRT(RJMIN2)
      SEMI1 = 2.0/RIJ - VREL2/(BODY(I) + BODY(JCL))
      SEMI1 = 1.0/SEMI1
      ECC1 = SQRT((1.0D0 - RIJ/SEMI1)**2 +
     &                     RDOT**2/(SEMI1*(BODY(I) + BODY(JCL))))
      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,JCL) - X(1,I1))**2 + (X(2,JCL) - X(2,I1))**2 +
     &                                 (X(3,JCL) - X(3,I1))**2
      RJI2 = (X(1,JCL) - X(1,I2))**2 + (X(2,JCL) - X(2,I2))**2 +
     &                                 (X(3,JCL) - X(3,I2))**2
*
*       Identify the dominant interaction (#I1 or #I2).
      FJI1 = (BODY(I1) + BODY(JCL))/RJI1
      FJI2 = (BODY(I2) + BODY(JCL))/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,JCL) - X(K,J))*(XDOT(K,JCL) - XDOT(K,J))
          RJJ = RJJ + (X(K,JCL) - X(K,J))**2
          VJJ = VJJ + (XDOT(K,JCL) - XDOT(K,J))**2
   15 CONTINUE
*
*       Determine vectorial perturbation on intruder and closest component.
      CALL FPERT(JCL,J,NNB,PERT)
      GJ = PERT*RJ2/(BODY(JCL) + BODY(J))
      RJ = SQRT(RJJ)
*     WRITE (6,18) NAME(J),NAME(JCL),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 #JCL & J.
      IF (TDOT2(IPAIR)/R(IPAIR).LT.RD/RJ) XF = 1.0
*
*       Increase the cross section for quadruples.
      IF (JCL.GT.N) THEN
          SEMI2 = -0.5*BODY(JCL)/H(JCL-N)
          APO = APO + ABS(SEMI2)
      END IF
*
      SEMIJ = 2.0/RJ - VJJ/(BODY(JCL) + 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
back to top