nblist.f

Go to the documentation of this file.

      SUBROUTINE nblist(I,RS0)
*
*
*       Neighbour list & radius.
*       ------------------------
*
      include 'common6.h'
*
*
*       Form square neighbour radius.
      rs2 = rs0**2
      ilist(1) = 0
*
*       Modify initial radius for Plummer or King-Michie models.
      IF (kz(5).GT.0.OR.kz(22).GE.2) THEN
          ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
     &                                   (x(3,i) - rdens(3))**2
*       Note possible large RS0 set in routine MERGE at later times.
          IF (time.EQ.0.0d0) THEN
              rs2 = rs2*(1.0 + ri2)
*         ELSE
*             RSX = 0.3*(1000.0/FLOAT(N - NPAIRS))**0.3333
*             RS2 = (RSX*RSCALE)**2*(1.0 + RI2)
          END IF
      END IF
*
      vi2 = xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
      IF (vi2.GT.0.0) THEN
          dtr = 0.1*rs0/sqrt(vi2)
          vrfac = -0.1*rs2/dtr
*       Estimated radial velocity factor for outer sphere.
      ELSE
          vrfac = 0.0
      END IF
*
*       Search all single particles & c.m. bodies.
    1 rcrit2 = 1.59*rs2
      nnb = 1
      DO 10 j = ifirst,ntot
          rij2 = (x(1,i) - x(1,j))**2 + (x(2,i) - x(2,j))**2 +
     &                                  (x(3,i) - x(3,j))**2
          IF (rij2.GT.rcrit2) go to 10
          IF (rij2.GT.rs2) THEN
*       Skip velocity test for primordial binaries.
              IF (time.LE.0.0d0.AND.j.LE.2*nbin0) go to 6
              rijdot = 0.0
              DO 4 k = 1,3
                  rijdot = rijdot +
     &                         (x(k,i) - x(k,j))*(xdot(k,i) - xdot(k,j))
    4         CONTINUE
*       Accept member if maximum penetration factor is significant.
              IF (rijdot.GE.vrfac) go to 10
          ELSE
              IF (j.EQ.i) go to 10
          END IF
*       Include both components in case of primordial binaries.
    6     IF (j.LE.2*nbin0.AND.time.EQ.0.0d0) THEN
              jpair = kvec(j)
*       Consider odd and even locations and skip body $I.
              IF (j.EQ.2*jpair-1) THEN
                  nnb = nnb + 1
                  ilist(nnb) = j
                  IF (i.NE.j + 1) THEN
                      nnb = nnb + 1
                      ilist(nnb) = j + 1
                  END IF
*       Exclude previously saved component.
              ELSE IF (ilist(nnb).LT.j-1) THEN
                  IF (i.NE.j - 1) THEN
                      nnb = nnb + 1
                      ilist(nnb) = j - 1
                  END IF
                  nnb = nnb + 1
                  ilist(nnb) = j
              END IF
          ELSE
              nnb = nnb + 1
              ilist(nnb) = j
          END IF
   10 CONTINUE
*
*       Check membership (NNB > 2 for primordial binaries).
      IF (nnb.EQ.1.OR.(i.LE.2*nbin0.AND.nnb.EQ.2)) THEN
*       Double the neighbour sphere and try again.
          rs2 = 1.59*rs2
          nbvoid = nbvoid + 1
          go to 1
      ELSE IF (nnb.GE.nnbmax) THEN
*       Reduce neighbour sphere but avoid possible looping.
          a1 = znbmax/float(nnb)
          IF (abs(a1 - 0.5).LT.0.05) a1 = 1.2*a1
          rs2 = rs2*a1**0.66667
          nbfull = nbfull + 1
          go to 1
      END IF
*
*       Set neighbour radius and copy list membership.
      rs(i) = sqrt(rs2)
      list(1,i) = nnb - 1
      DO 20 l = 2,nnb
          list(l,i) = ilist(l)
   20 CONTINUE
*
      RETURN
*
      END