sweep.f

Go to the documentation of this file.

      SUBROUTINE sweep(DTCL,RCL)
*
*
*       Enforced KS regularization of wide binaries.
*       --------------------------------------------
*
      include 'common6.h'
      SAVE i
*
*
*       Identify wider binaries giving rise to systematic errors.
      zmx = 100.0*bodym
      i = ifirst - 1
    1 i = i + 1
      IF (i.GT.n) go to 30
      IF (body(i).GT.zmx) go to 2
      IF (step(i).GT.dtcl.OR.body(i).EQ.0.0d0) go to 1
    2 CONTINUE
*
*       Search neighbour list of all KS candidates (STEP < DTCL).
      rx2 = 100.0
      nnb1 = list(1,i) + 1
      IF (nnb1.EQ.1) go to 1
      DO 10 l = 2,nnb1
          j = list(l,i)
          rij2 = 0.0
          DO 5 k = 1,3
              rij2 = rij2 + (x(k,i) - x(k,j))**2
    5     CONTINUE
          IF (rij2.LT.rx2) THEN
              rx2 = rij2
              jmin = j
          END IF
   10 CONTINUE
*
*       Skip any close c.m./chain body (small STEP treated by IMPACT).
      IF (jmin.GT.n.OR.name(jmin).LE.0) go to 1
*
*       Form inverse semi-major axis.
      vij2 = 0.0
      DO 15 k = 1,3
          vij2 = vij2 + (xdot(k,i) - xdot(k,jmin))**2
   15 CONTINUE
      ainv = 2.0/sqrt(rx2) - vij2/(body(i) + body(jmin))
*
      IF (ainv.GT.0.0) THEN
          semi = 1.0/ainv
          zmb = body(i) + body(jmin)
          erel = 0.5*zmb/semi
      ELSE
          erel = 0.0
      END IF
*       Initialize bound KS pairs outside standard parameters.
      IF (ainv.GT.1.0/rcl.OR.erel.GT.eclose) THEN
          icomp = min(i,jmin)
          jcomp = max(i,jmin)
*       Skip possible rare case of chain c.m. forming binary.
          IF (name(icomp).EQ.0.OR.name(jcomp).EQ.0) go to 1
*       Ensure most recent velocity used for new KS.
          DO 16 k = 1,3
              x0dot(k,icomp) = xdot(k,icomp)
              x0dot(k,jcomp) = xdot(k,jcomp)
   16     CONTINUE
          CALL ksreg
          newks = newks + 1
          ri2 = 0.0
          DO 20 k = 1,3
              ri2 = ri2 + (x(k,ntot) - rdens(k))**2
   20     CONTINUE
          semi = 1.0/ainv
          ecc2 = (1.0 - r(npairs)/semi)**2 +
     &                            tdot2(npairs)**2/(body(ntot)*semi)
          ecc = sqrt(ecc2)
          j1 = 2*npairs - 1
          IF (newks.LT.50.OR.erel.GT.eclose) THEN
             WRITE (6,25)   name(j1), name(j1+1), list(1,j1),
     &                      list(1,ntot), sqrt(ri2), ecc, semi,
     &                      gamma(npairs)
   25        FORMAT (' ENFORCED KS    NAM NP NNB r E A GAM ',
     &                                2i7,i4,i5,f7.2,f8.3,1p,e10.2,e9.1)
          END IF
      END IF
      IF (i.LT.n) go to 1
*
   30 RETURN
*
      END