ficorr.f

Go to the documentation of this file.

      SUBROUTINE ficorr(I,DM)
*
*
*       Local force corrections due to mass loss.
*       -----------------------------------------
*
      include 'common6.h'
      REAL*8  a(6)
*
*
*       Correct force & first derivative for neighbours of body #I.
      nnb1 = list(1,i) + 1
      DO 20 l = 2,nnb1
          j = list(l,i)
          nnb2 = list(1,j) + 1
*
*       Ensure that body #I is a neighbour of body #J.
          DO 1 k = 2,nnb2
              IF (list(k,j).EQ.i) go to 5
    1     CONTINUE
          go to 20
*
    5     rij2 = 0.0d0
          a7 = 0.0d0
*
          DO 10 k = 1,3
              a(k) = x(k,i) - x(k,j)
              a(k+3) = xdot(k,i) - xdot(k,j)
              rij2 = rij2 + a(k)**2
              a7 = a7 + a(k)*a(k+3)
   10     CONTINUE
*
          a7 = 3.0*a7/rij2
          a8 = dm/(rij2*sqrt(rij2))
*
          DO 15 k = 1,3
              a(k+3) = (a(k+3) - a7*a(k))*a8
              f(k,j) = f(k,j) - 0.5*a(k)*a8
              fi(k,j) = fi(k,j) - a(k)*a8
              fdot(k,j) = fdot(k,j) - one6*a(k+3)
              d1(k,j) = d1(k,j) - a(k+3)
              fidot(k,j) = fidot(k,j) - a(k+3)
   15     CONTINUE
   20 CONTINUE
*
*       Obtain the potential energy due to all particles.
      potj = 0.0d0
      DO 30 j = ifirst,ntot
          IF (j.EQ.i) go to 30
          rij2 = 0.0d0
          DO 25 k = 1,3
              rij2 = rij2 + (x(k,i) - x(k,j))**2
   25     CONTINUE
          potj = potj + body(j)/sqrt(rij2)
   30 CONTINUE
*
*       Update the potential energy loss.
      emdot = emdot - dm*potj
*
*       Include kinetic energy correction.
      vi2 = xdot(1,i)**2 + xdot(2,i)**2 + xdot(3,i)**2
      emdot = emdot + 0.5*dm*vi2
*
*       See whether linearized tidal terms should be included.
      IF (kz(14).GT.0.AND.kz(14).LT.3) THEN
          ecdot = ecdot - 0.5*dm*(tidal(1)*x(1,i)**2 +
     &                            tidal(3)*x(3,i)**2)
      END IF
*
*       Check optional Plummer potential.
      IF (kz(14).EQ.4.OR.kz(14).EQ.3) THEN
          ri2 = ap2
          DO 50 k = 1,3
              ri2 = ri2 + x(k,i)**2
  50      CONTINUE
          ecdot = ecdot - dm*mp/sqrt(ri2)
      END IF
*
      RETURN
*
      END