kslist.f

Go to the documentation of this file.

      SUBROUTINE kslist(IPAIR)
*
*
*       KS perturber selection.
*       -----------------------
*
      include 'common6.h'
*
*
*       Set component & c.m. index and form semi-major axis & eccentricity.
      i1 = 2*ipair - 1
      i = n + ipair
      semi = -0.5d0*body(i)/h(ipair)
      eb = -0.5*body(i1)*body(i1+1)/semi
      ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
      ecc = sqrt(ecc2)
      it = 0
*
*       Use semi-major axis and/or RMIN for perturber selection.
      IF (eb.LT.ebh) THEN
          rap = semi*(1.0 + ecc)
      ELSE
*       Include a tapered criterion depending on energy for soft binaries.
          IF (eb.LT.0.0) THEN
              zfac = 1.0 + abs(eb - ebh)/abs(ebh)
          ELSE
              zfac = 1.0
          END IF
*       Adopt actual apocentre for perturber selection if R > SEMI.
          IF (semi.GT.0.0d0.AND.r(ipair).GT.semi) THEN
              rap = semi*(1.0 + ecc)
          ELSE
*       Accept separation limit for hyperbolic cases.
              rap = max(zfac*semi,r(ipair))
          END IF
*       Ensure extra perturbers at new regularizations (R may be small).
          IF (iphase.GT.0.AND.semi.GT.0.0) THEN
              rap = semi*(1.0 + ecc)
          END IF
      END IF
*
*       Restrict perturber selection for massive eccentric binary.
      IF (body(i).GT.100.0*bodym.AND.
     &   (r(ipair).LT.0.1*semi.OR.semi.LT.0)) THEN
          rap = min(rmin,0.1*abs(semi))
      END IF
*
      rcrit2 = 2.0*rap**2/body(i)
      rcrit3 = rcrit2*rap/gmin
*       Base fast search on maximum binary mass (BODY1).
      rcrit2 = rcrit2*body1*cmsep2
*
*       Select new perturbers from the neighbour list.
    6 nnb1 = 1
      nnb2 = list(1,i) + 1
      DO 10 l = 2,nnb2
          j = list(l,i)
          w1 = x(1,j) - x(1,i)
          w2 = x(2,j) - x(2,i)
          w3 = x(3,j) - x(3,i)
          rsep2 = w1*w1 + w2*w2 + w3*w3
*       Include any merged c.m. or chain c.m. bodies in the fast test.
          IF (rsep2.LT.rcrit2.OR.name(j).LE.0) THEN
              rij3 = rsep2*sqrt(rsep2)
*       Estimate unperturbed distance from tidal limit approximation.
              IF (rij3.LT.body(j)*rcrit3) THEN
                  nnb1 = nnb1 + 1
                  list(nnb1,i1) = j
              ELSE IF (j.GT.n) THEN
*       Employ a more generous criterion for possible wide binary.
                  rj = min(10.0*abs(semi),-body(j)/h(j-n))
                  IF (rsep2.LT.cmsep2*rj**2) THEN
                      nnb1 = nnb1 + 1
                      list(nnb1,i1) = j
                  END IF
              END IF
          END IF
   10 CONTINUE
*
*       Ensure at least one perturber first time (max 5 tries except CHAOS).
      IF (nnb1.EQ.1.AND.iphase.GT.0.AND.nnb2.GT.1.AND.time.GT.0.0) THEN
          rcrit2 = 2.0*rcrit2
          rcrit3 = 2.0*rcrit3
          it = it + 1
*       Skip repeat for small size (next KSLIST requires many periods).
          IF ((semi*su.GT.10.0.AND.it.LE.5).OR.kstar(i).EQ.-1) go to 6
      END IF
*
*       Check case of no perturbers (dual purpose).
      IF (nnb1.EQ.1) THEN
*       Add distant perturber for hyperbolic orbit.
          IF (semi.LT.0.0) THEN
              nnb1 = 2
              list(2,i1) = list(2,i)
              go to 20
          END IF
*         ELSE IF (KZ(27).LE.0) THEN
          IF (kz(27).LE.0) THEN
*       Specify one unperturbed period at apocentre (NB! check STEP(I)).
              step(i1) = twopi*semi*sqrt(semi/body(i))
              step(i1) = min(step(i1),step(i))
          ELSE
*       Maintain perturbed motion during Chaos event (not ROCHE/SPIRAL).
              IF (kstar(i).EQ.-1) THEN
                  IF (list(1,i1).GT.0) THEN
                      nnb1 = 2
                      list(2,i1) = n
                  END IF
              ELSE
                  step(i1) = twopi*semi*sqrt(semi/body(i))
                  step(i1) = min(step(i1),step(i))
              END IF
          END IF
      END IF
*
*       Save perturber membership.
   20 list(1,i1) = nnb1 - 1
*
      RETURN
*
      END