fpert.f

Go to the documentation of this file.

      SUBROUTINE fpert(I,J,NP,PERT)
*
*
*       Perturbing force on dominant two-body system.
*       ---------------------------------------------
*
      include 'common6.h'
      REAL*8  fp(3)
*
*
      DO 1 k = 1,3
          fp(k) = 0.0d0
    1 CONTINUE
*
*       Obtain perturbation on body #I & J due to NP members of JLIST.
      rcrit2 = 1.0d+04*rmin2
      iter = 0
    2 DO 10 kk = 1,np
          k = jlist(kk)
          IF (k.EQ.i.OR.k.EQ.j) go to 10
          l = i
    5     a1 = x(1,k) - x(1,l)
          a2 = x(2,k) - x(2,l)
          a3 = x(3,k) - x(3,l)
          rij2 = a1**2 + a2**2 + a3**2
*       Restrict evaluation RCRIT2 (neighbour list may be used).
          IF (l.EQ.i.AND.rij2.GT.rcrit2) go to 10
          a4 = body(k)/(rij2*sqrt(rij2))
          IF (l.EQ.j) a4 = -a4
          fp(1) = fp(1) + a1*a4
          fp(2) = fp(2) + a2*a4
          fp(3) = fp(3) + a3*a4
          IF (l.EQ.i) THEN
              l = j
              go to 5
          END IF
   10 CONTINUE
*
*       Increase search distance on zero contributions (< 10 times).
      IF (fp(1).EQ.0.0d0) THEN
          rcrit2 = 4.0*rcrit2
          iter = iter + 1
          IF (iter.LT.10) go to 2
      END IF
*
      pert = sqrt(fp(1)**2 + fp(2)**2 + fp(3)**2)
*
      RETURN
*
      END