cmcorr.f

Go to the documentation of this file.

      SUBROUTINE cmcorr
*
*
*       Center of mass & total force corrections.
*       -----------------------------------------
*
      include 'common6.h'
*
*
*       Initialize centre of mass variables.
      DO 10 k = 1,3
          cmr(k) = 0.0d0
          cmrdot(k) = 0.0d0
   10 CONTINUE
*
*       Form c.m. coordinate & velocity displacements.
      DO 20 i = ifirst,ntot
          DO 15 k = 1,3
              cmr(k) = cmr(k) + body(i)*x(k,i)
              cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
   15     CONTINUE
   20 CONTINUE
*
      DO 30 k = 1,3
          cmr(k) = cmr(k)/zmass
          cmrdot(k) = cmrdot(k)/zmass
   30 CONTINUE
*
*       Include effect of c.m. motion in Plummer potential.
      IF (kz(14).EQ.4) THEN
          CALL plpot1(phi1)
      END IF
*
*       Apply c.m. corrections to X & XDOT and accumulate energy changes.
      errx = 0.0d0
      errv = 0.0d0
      DO 40 i = ifirst,ntot
          DO 35 k = 1,3
              xi2 = x(k,i)**2
              vi2 = xdot(k,i)**2
              x(k,i) = x(k,i) - cmr(k)
              xdot(k,i) = xdot(k,i) - cmrdot(k)
*       Note TIDAL(K) = 0 for KZ(14) = 3, hence no skip on ERRX is needed.
              errx = errx - tidal(k)*body(i)*(x(k,i)**2 - xi2)
              errv = errv + body(i)*(xdot(k,i)**2 - vi2)
   35     CONTINUE
   40 CONTINUE
*
*       Adjust the total energy to new kinetic energy & tidal potential.
      be(3) = be(3) + 0.5*(errx + errv)
      e(11) = e(11) - 0.5*(errx + errv)
*
*       Perform a consistent shift of the density centre.
      DO 50 k = 1,3
          rdens(k) = rdens(k) - cmr(k)
   50 CONTINUE
*
*       Subtract tidal corrections from total force & first derivative.
      IF (kz(14).GT.0.AND.kz(14).LE.2) THEN
          DO 60 i = ifirst,ntot
*       Skip ghosts to avoid spurious prediction inside 1.0E+10.
              IF (body(i).EQ.0.0d0) go to 60
              DO 55 k = 1,3
                  df = tidal(k)*cmr(k)
                  dd = tidal(k)*cmrdot(k)
                  fr(k,i) = fr(k,i) - df
                  f(k,i) = f(k,i) - 0.5*df
                  d1r(k,i) = d1r(k,i) - dd
                  frdot(k,i) = frdot(k,i) - dd
                  fdot(k,i) = fdot(k,i) - one6*dd
                  IF (k.EQ.1) THEN
                      fi(1,i) = fi(1,i) - tidal(4)*cmrdot(2)
                  ELSE IF (k.EQ.2) THEN
                      fi(2,i) = fi(2,i) + tidal(4)*cmrdot(1)
                  END IF
   55         CONTINUE
   60     CONTINUE
      END IF
*
*       Re-determine X0 & X0DOT consistently with current corrected X & XDOT.
      DO 70 i = ifirst,ntot 
          IF (body(i).EQ.0.0d0) go to 70
          dt = time - t0(i)
          dtr = time - t0r(i)
          a1 = 0.2*dt
          a2 = dt/24.0
*
          DO 65 k = 1,3
              f2dotk = d3r(k,i)*dtr + d2r(k,i) + d3(k,i)*dt + d2(k,i)
              f3dotk = d3r(k,i) + d3(k,i)
              dv0 = (((f3dotk*a2 + one6*f2dotk)*dt +
     &                                3.0*fdot(k,i))*dt + 2.0*f(k,i))*dt
              x0dot(k,i) = xdot(k,i) - dv0
              dx0 = ((((f3dotk*a1 + f2dotk)*a2 + fdot(k,i))*dt +
     &                                       f(k,i))*dt + x0dot(k,i))*dt
              x0(k,i) = x(k,i) - dx0
              x(k,i) = x0(k,i)
              xdot(k,i) = x0dot(k,i)
   65     CONTINUE
   70 CONTINUE
*
*       Check differential correction for Plummer potential.
      IF (kz(14).EQ.4) THEN
          CALL plpot1(phi2)
          emdot = emdot + (phi1 - phi2)
      END IF
*
*       Ensure consistent coordinates & velocities for binary components.
      DO 80 ipair = 1,npairs
          IF (body(n+ipair).GT.0.0d0) THEN
              CALL resolv(ipair,1)
          END IF
   80 CONTINUE
*
      RETURN
*
      END