jacobi.f

Go to the documentation of this file.

      SUBROUTINE jacobi(NESC)
*
*
*       Jacobi escape criterion.
*       ------------------------
*
      include 'common6.h'
*
*
*       Specify escape energy (tidal field or isolated system).
      IF (kz(14).GT.0.AND.kz(14).LE.2) THEN
          ecrit = -1.5*(tidal(1)*zmass**2)**0.333
      ELSE
          ecrit = 0.0
      END IF
*       Define current mean square velocity (= 0.5 initially).
      v2m = 0.5*zmass/rscale
*
*       Count all escapers.
      nesc = 0
      DO 60 i = ifirst,ntot
          ri2 = (x(1,i) - rdens(1))**2 + (x(2,i) - rdens(2))**2 +
     &                                   (x(3,i) - rdens(3))**2
          vi2 = xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
          IF (ri2.GT.rscale**2) THEN
              vc2 = zmass/sqrt(ri2) + ecrit
              IF (vi2.LT.vc2) go to 60
          ELSE
*       Skip velocities below sqrt(2) times equilibrium value.
              IF (vi2.LT.2.0*v2m + ecrit) go to 60
          END IF
          poti = 0.0
          DO 55 j = ifirst,ntot
              IF (j.EQ.i) go to 55
              rij2 = 0.0
              DO 52 k = 1,3
                  rij2 = rij2 + (x(k,i) - x(k,j))**2
   52         CONTINUE
              poti = poti + body(j)/sqrt(rij2)
   55     CONTINUE
          ei = 0.5*vi2 - poti
          IF (ei.GT.ecrit) nesc = nesc + 1
   60 CONTINUE
*
      RETURN
*
      END