ksres2.f

Go to the documentation of this file.

      SUBROUTINE ksres2(J,J1,J2,RIJ2)
*
*
*       Coordinates & velocities of KS pair.
*       ------------------------------------
*
      include 'common6.h'
      common/slow0/  range,islow(10)
      REAL*8  ui(4),rdot(3),v(4),a1(3,4)
*
*
*       Skip merged ghost binary (defined by zero total mass).
      IF (body(n+j).LE.0.0d0) go to 40
*
*       Resolve components of pair #J at curent time.
      j2 = j + j
      j1 = j2 - 1
      a2 = 1.0/r(j)
      a3 = a2*(time - t0(j1))
      IF (kslow(j).GT.1) THEN
          imod = kslow(j)
          a3 = a3/float(islow(imod))
      END IF
*
*       Decide appropriate order for interpolation & prediction.
      IF (rij2.GT.625.0*r(j)**2) THEN
*       Convert physical interval to regularized time (second order only).
          dtu = (1.0 - 0.5d0*tdot2(j)*a2*a3)*a3
          IF (abs(dtu).GT.dtau(j)) dtu = dtau(j)
*
*       Predict regularized coordinates of distant pair to second order.
          DO 10 k = 1,4
              ui(k) = (fu(k,j)*dtu + udot(k,j))*dtu + u0(k,j)
              v(k) = (3.0*fudot(k,j)*dtu + 2.0*fu(k,j))*dtu + udot(k,j)
   10     CONTINUE
      ELSE
*        Expand regularized time interval to third order.
          a4 = 3.0d0*tdot2(j)**2*a2 - tdot3(j)
          dtu = ((one6*a4*a3 - 0.5d0*tdot2(j))*a2*a3 + 1.0)*a3
*       Apply safety test near small pericentre or for unperturbed motion.
          IF (dtu.GT.dtau(j)) dtu = 0.8*dtau(j)
          IF (list(1,j1).EQ.0) dtu = 0.0
*
*       Predict regularized coordinates to third order.
          dtu1 = dtu/24.0d0
          dtu2 = one6*dtu
          DO 20 k = 1,4
              ui(k) = (((fudot2(k,j)*dtu1 + fudot(k,j))*dtu +
     &                           fu(k,j))*dtu + udot(k,j))*dtu + u0(k,j)
              v(k) = ((fudot2(k,j)*dtu2 + 3.0d0*fudot(k,j))*dtu +
     &                                    2.0d0*fu(k,j))*dtu + udot(k,j)
   20     CONTINUE
      END IF
*
*       Employ KS transformation.
      q1 = ui(1)**2 - ui(2)**2 - ui(3)**2 + ui(4)**2
      q2 = ui(1)*ui(2) - ui(3)*ui(4)
      q3 = ui(1)*ui(3) + ui(2)*ui(4)
      q2 = q2 + q2
      q3 = q3 + q3
      j3 = n + j
      a2 = body(j2)/body(j3)
*
*       Set global coordinates of regularized components.
      x(1,j1) = x(1,j3) + a2*q1
      x(2,j1) = x(2,j3) + a2*q2
      x(3,j1) = x(3,j3) + a2*q3
      x(1,j2) = x(1,j1) - q1
      x(2,j2) = x(2,j1) - q2
      x(3,j2) = x(3,j1) - q3
*
*       Set current transformation matrix and two-body separation.
      CALL matrix(ui,a1)
      ri = ui(1)**2 + ui(2)**2 + ui(3)**2 + ui(4)**2
      rinv = 2.0/ri
*
*       Obtain relative velocities from KS transformation.
      DO 30 l = 1,3
          rdot(l) = 0.0d0
          DO 25 k = 1,4
              rdot(l) = rdot(l) + a1(l,k)*v(k)*rinv
   25     CONTINUE
   30 CONTINUE
*
*       Set global velocities of KS components.
      xdot(1,j1) = xdot(1,j3) + a2*rdot(1)
      xdot(2,j1) = xdot(2,j3) + a2*rdot(2)
      xdot(3,j1) = xdot(3,j3) + a2*rdot(3)
      xdot(1,j2) = xdot(1,j1) - rdot(1)
      xdot(2,j2) = xdot(2,j1) - rdot(2)
      xdot(3,j2) = xdot(3,j1) - rdot(3)
*
   40 RETURN
*
      END