ksres.f

Go to the documentation of this file.

      SUBROUTINE ksres(J,J1,J2,RIJ2)
*
*
*       Coordinates of regularized pair.
*       --------------------------------
*
      include 'common6.h'
      common/slow0/  range,islow(10)
      REAL*8  ui(4)
*
*
*       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)
   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)
*       Ignore the expansion for collision (already at pericentre).
          IF (list(1,j1).EQ.0.OR.iphase.EQ.9) dtu = 0.0
*
*       Predict regularized coordinates to fourth order.
          dtu1 = dtu/24.0d0
          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)
   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
*
      RETURN
*
      END