nbtide.f

Go to the documentation of this file.

      SUBROUTINE nbtide(I,J,RJMIN2)
*
*
*       Close two-body interaction.
*       ---------------------------
*
      include 'common6.h'
      REAL*8  xi(3),xj(3),vi(3),vj(3),vcm(3),vr(3)
*
*
*       Copy coordinates of close bodies and predict velocity for body #J.
      rdot = 0.0
      vrel2 = 0.0
      dt = time - t0(j)
      DO 10 k = 1,3
          xi(k) = x(k,i)
          xj(k) = x(k,j)
          vi(k) = x0dot(k,i)
          vj(k) = (3.0*fdot(k,j)*dt + 2.0*f(k,j))*dt + x0dot(k,j)
          rdot = rdot + (xi(k) - xj(k))*(vi(k) - vj(k))
          vrel2 = vrel2 + (vi(k) - vj(k))**2
   10 CONTINUE
*
*       Only consider approaching bodies.
      IF (rdot.LE.0.0) go to 100
*
*       Predict coordinates & velocities at beginning of step for body #I.
      rdot0 = 0.0
      dti = time - t0(i) - step(i)
      dtj = time - t0(j) - step(i)
      DO 20 k = 1,3
          xi(k) = ((fdot(k,i)*dti + f(k,i))*dti + x0dot(k,i))*dti +
     &                                                           x0(k,i)
          xj(k) = ((fdot(k,j)*dtj + f(k,j))*dtj + x0dot(k,j))*dtj +
     &                                                           x0(k,j)
          vi(k) = (3.0*fdot(k,i)*dti + 2.0*f(k,i))*dti + x0dot(k,i)
          vj(k) = (3.0*fdot(k,j)*dtj + 2.0*f(k,j))*dtj + x0dot(k,j)
          rdot0 = rdot0 + (xi(k) - xj(k))*(vi(k) - vj(k))
   20 CONTINUE
*
*       Check pericentre condition (old radial velocity < 0).
      IF (rdot0.GT.0.0) go to 100
*
      zm = body(i) + body(j)
      erel = 0.5*vrel2 - zm/sqrt(rjmin2)
*
*       Specify the energy loss (experimental).
      dh = 0.1*vrel2/2.0
      hnew = erel - dh
      semi0 = -0.5*zm/erel
      semi = -0.5*zm/hnew
      WRITE (6,50)  name(i), name(j), semi0, semi, sqrt(rjmin2), dh
   50 FORMAT (5x,'NBTIDE:  NAMES A0 A RIJ DH  ',2i5,2f10.5,f8.4,f8.3)
*
*       Skip energy loss treatment unless final orbit is significantly bound.
      IF (semi.LT.0.0.OR.semi.GT.4.0*rmin) go to 100
*
*       Predict coordinates & velocity for body #J to highest order.
      CALL xvpred(j,0)
*
*       Set current velocities and form c.m. velocity.
      DO 30 k = 1,3
          vi(k) = xdot(k,i)
          vj(k) = xdot(k,j)
          vcm(k) = (body(i)*xdot(k,i) + body(j)*xdot(k,j))/zm
   30 CONTINUE
*
*       Introduce velocity change for each component and initialize X0DOT.
      fac = sqrt((vrel2 - 2.0d0*dh)/vrel2)
      DO 40 k = 1,3
          vr(k) = fac*(vi(k) - vj(k))
          xdot(k,i) = vcm(k) + body(j)*vr(k)/zm
          xdot(k,j) = vcm(k) - body(i)*vr(k)/zm
          x0dot(k,i) = vi(k)
          x0dot(k,j) = vj(k)
   40 CONTINUE
*
      de = body(i)*body(j)*dh/zm
*       Increase event counter and update total energy loss.
      ndiss = ndiss + 1
      ecoll = ecoll + de
*
*       Set components and phase indicator for new KS regularization.
      icomp = min(i,j)
      jcomp = max(i,j)
      iphase = 1
*
      WRITE (6,60)  name(icomp), name(jcomp), semi, be(3) + de , de
   60 FORMAT (' TIDAL CAPTURE   NM A E DE ',2i5,f8.4,2f11.6)
*
  100 RETURN
*
      END