hivel.f

Go to the documentation of this file.

      SUBROUTINE hivel(IH)
*
*
*       High-velocity particle search.
*       ------------------------------
*
      include 'common6.h'
*
*
*       Copy membership and add 1 for analogy with HARP/GRAPE.
      nhi = listv(1) + 1
*
*       Specify square velocity limit in terms of current state.
      vmax2 = 16.0*eclose
*
*       Check for removal of distant high-velocity particle.
      IF (ih.LT.0) THEN
    1     lh = 0
          DO 2 l = 2,nhi
              i = listv(l)
              vi2 = x0dot(1,i)**2 + x0dot(2,i)**2 + x0dot(3,i)**2
              ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
     &                                       (x(3,i) - rdens(3))**2
*       Save index of fast particle outside 3*<R> or VI < VMAX/4.
              IF (ri2.GT.9.0*rscale**2.OR.vi2.LT.vmax2/16.0) THEN
                  lh = l
                  li = i
                  rl2 = ri2
                  vl2 = vi2
              END IF
    2     CONTINUE
*       Reduce membership and remove any distant member from the list.
          IF (lh.GT.0) THEN
              WRITE (29,3)  time+toff, li, name(li), sqrt(rl2),
     &                      sqrt(vl2)
    3         FORMAT (' HIVEL REMOVE    T I NAM R V ',f12.4,2i6,2f6.1)
              CALL flush(29)
              nhi = nhi - 1
              listv(1) = listv(1) - 1
              DO 4 l = lh,nhi
                  listv(l) = listv(l+1)
    4         CONTINUE
              go to 1
          END IF
      END IF
*
*       Set index to known neutron star or terminated KS/chain.
      IF (ih.GT.0) THEN
          i1 = ih
          i2 = ih
      ELSE
*       Include two first single particles or up to five chain members.
          i1 = ifirst
          i2 = ifirst + 1
          IF (iphase.EQ.8) THEN
              i2 = ifirst + 4
*       Search all particles after escaper removal (defined by IPHASE = -1).
          ELSE IF (ih.EQ.0.AND.iphase.EQ.-1) THEN
              i2 = ntot
              nhi = 1
              listv(1) = 0
          ELSE IF (iphase.EQ.1) THEN
*       See whether the first few locations may have been exchanged.
              DO 5 l = 2,nhi
                  IF (listv(l).LT.ifirst + 3) THEN
                      i2 = ntot
                  END IF
    5         CONTINUE
          END IF
      END IF
*
*       Add any new high-velocity particles (save F**2 > N & STEP < DTMIN).
      nhv = 0
      DO 10 i = i1,i2
          fi2 = f(1,i)**2 + f(2,i)**2 + f(3,i)**2
          vi2 = x0dot(1,i)**2 + x0dot(2,i)**2 + x0dot(3,i)**2
          ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
     &                                   (x(3,i) - rdens(3))**2
*       Form a list of recently ejected candidates.
          IF (fi2.GT.float(n).OR.step(i).LT.dtmin.OR.
     &        stepr(i).LT.20.0*dtmin) THEN
              IF (vi2.GT.vmax2) THEN
                  nhv = nhv + 1
                  jlist(nhv) = i
              END IF
              go to 10
          END IF
*      Include all VMAX2 members inside 2*RSCALE or with negative velocity.
          rd = x(1,i)*xdot(1,i) + x(2,i)*xdot(2,i) + x(3,i)*xdot(3,i)
          IF (vi2.GT.vmax2.AND.(ri2.LT.4.0*rscale**2.OR.rd.LT.0.)) THEN
              DO 8 l = 2,nhi
                  IF (i.EQ.listv(l)) go to 10
    8         CONTINUE
*       Check maximum membership and possible ghost particle.
              IF (nhi.GE.mlv.OR.step(i).GT.1.0) go to 10
              nhi = nhi + 1
              listv(1) = listv(1) + 1
              nfast = nfast + 1
              listv(nhi) = i
              WRITE (29,9)  time+toff, nhi-1, i, name(i), kstar(i),
     &                      sqrt(vi2), sqrt(ri2), step(i)
    9         FORMAT (' HIVEL ADD    T NHI I NM K* VI R DT ',
     &                               f12.4,i4,2i6,i4,2f6.2,1p,e10.2)
              CALL flush(29)
          END IF
   10 CONTINUE
*
*       Consider single fast particle or hyperbolic two-body motion.
      IF ((nhv.EQ.1.OR.nhv.EQ.2).AND.nhi.LE.mlv-nhv) THEN
*       Compare any candidates with existing members.
          DO 20 k = 1,nhv
              DO 15 l = 2,nhi
                  IF (jlist(k).EQ.listv(l)) go to 30
   15         CONTINUE
   20     CONTINUE
          i1 = jlist(1)
*       Include single particles without further tests.
          IF (nhv.EQ.1) THEN
              nhi = nhi + 1
              listv(1) = listv(1) + 1
              nfast = nfast + 1
              listv(nhi) = i1
              ri2 = (x(1,i1) - rdens(1))**2 + (x(2,i1) - rdens(2))**2 +
     &                                        (x(3,i1) - rdens(3))**2
              vi2 = x0dot(1,i1)**2 + x0dot(2,i1)**2 + x0dot(3,i1)**2
              WRITE (29,22)  ttot, nhi-1, name(i1), iphase, sqrt(vi2),
     &                       sqrt(ri2), step(i1)
   22         FORMAT (' HIVEL ADD    T NHI NM IPH VI R DT ',
     &                               f12.4,i4,i6,i4,2f6.2,1p,e10.2)
              go to 30
          END IF
*       Evaluate two-body energy.
          i2 = jlist(2)
          rij2 = 0.0
          vij2 = 0.0
          rdot = 0.0
          DO 25 k = 1,3
              rij2 = rij2 + (x(k,i1) - x(k,i2))**2
              vij2 = vij2 + (xdot(k,i1) - xdot(k,i2))**2
              rdot = rdot + (x(k,i1) - x(k,i2))*(xdot(k,i1)-xdot(k,i2))
   25     CONTINUE
          rij = sqrt(rij2)
          semi = 2.0/rij - vij2/(body(i1) + body(i2))
*       Accept outwards hyperbolic motion arising from recent interaction.
          IF (semi.LT.0.0.AND.rij.LT.10.0*rmin.AND.rdot.GT.0.0) THEN
              nhi = nhi + 1
              listv(nhi) = i1
              nhi = nhi + 1
              listv(nhi) = i2
              listv(1) = listv(1) + 2
              nfast = nfast + 2
              WRITE (29,28) ttot, nhi-1, name(i1), name(i2), iphase, rij
   28         FORMAT (' HIVEL ADD    T NHI NM IPH RIJ ',
     &                               f12.4,i4,2i6,i4,1p,e10.2)
          END IF
      END IF
*
   30 CONTINUE
*       Increase counter if loop is over all particles.
*     IF (I2.EQ.NTOT) THEN
*         NHIVEL = NHIVEL + 1
*     END IF
*
      RETURN
*
      END