SUBROUTINE KSCORR(IPAIR,UI,UIDOT,FP,FD,TD2,TDOT4,TDOT5,TDOT6)
*
*
*       Corrector for KS motion.
*       ------------------------
*
      INCLUDE 'common6.h'
      COMMON/SLOW0/  RANGE,ISLOW(10)
      REAL*8  UI(4),UIDOT(4),FP(6),FD(6),FREG(4),FRD(4),A1(3,4),A(8),
     &        PR(4),PRD(4),U2(4),U3(4),U4(4),U5(4),FD1(3)
      PARAMETER  (ITLOW=2)
*
*
*       Convert from physical to regularized derivative using T' = R.
      RI = R(IPAIR)
      DO 1 K = 1,3
          FD1(K) = FD(K)
          FD(K) = RI*FD(K)
    1 CONTINUE
*
*       Include KS slow-down factor in the perturbation if ZMOD > 1.
      IF (KZ(26).GT.0) THEN
          IMOD = KSLOW(IPAIR)
          IF (IMOD.GT.1) THEN
              ZMOD = FLOAT(ISLOW(IMOD))
              DO 5 K = 1,3
                  FP(K) = ZMOD*FP(K)
                  FD(K) = ZMOD*FD(K)
                  FD1(K) = ZMOD*FD1(K)
    5         CONTINUE
          END IF
      END IF
*
*       Predict DH from Taylor series derivatives and save temporary value.
      DTU = DTAU(IPAIR)
      DT3 = ONE6*DTU
      DT4 = 0.25D0*DTU
      DH = (((HDOT4(IPAIR)*DT4 + HDOT3(IPAIR))*DT3 +
     &                     0.5D0*HDOT2(IPAIR))*DTU + HDOT(IPAIR))*DTU
*
*       Set time-step factors and copy Stumpff coefficients to scalars.
      DTSQ = DTU**2
      DT6 = 6.0/(DTU*DTSQ)
      DT2 = 2.0/DTSQ
      Z3 = SF(3,IPAIR)
      Z4 = SF(4,IPAIR)
      Z5 = SF(5,IPAIR)
*
*       Define optimized scalars outside iteration loop.
      DT5 = 0.2D0*DTU*Z5
      DTZ = DT4*Z4
      DH0 = (0.5D0*HDOT2(IPAIR)*DTU + HDOT(IPAIR))*DTU
*
*       Use extra iteration and new KSPERT for large perturbation.
*     IF (GAMMA(IPAIR).LT.1.0D-20) THEN
          ITMAX = ITLOW
*     ELSE
*         ITMAX = ITLOW + 1
*         I1 = 2*IPAIR - 1
*         NNB0 = LIST(1,I1)
*         I = N + IPAIR
*         BODYIN = 1.0/BODY(I)
*     END IF
*
*       Perform iteration with or without re-calculating perturbation.
      DO 40 ITER = 1,ITMAX
*
*       Obtain new transformation matrix.
          CALL MATRIX(UI,A1)
*
*       Form twice regularized force and half first derivative of H.
          HD = 0.0D0
          TD2 = 0.0D0
          DO 10 K = 1,4
              A(K) = A1(1,K)*FP(1) + A1(2,K)*FP(2) + A1(3,K)*FP(3)
              A(K+4) = A1(1,K)*FD(1) + A1(2,K)*FD(2) + A1(3,K)*FD(3)
              PR(K) = RI*A(K)
              FREG(K) = DH*UI(K) + PR(K)
              HD = HD + UIDOT(K)*A(K)
              TD2 = TD2 + UI(K)*UIDOT(K)
   10     CONTINUE
*
*       Set regularized velocity matrix (Levi-Civita matrix not required).
          CALL MATRIX(UIDOT,A1)
*
*       Include the whole (L*F)' term in explicit derivatives of FU & H'.
          HD2 = 0.0D0
          DO 15 K = 1,4
              AK4 = A(K+4) + A1(1,K)*FP(1) + A1(2,K)*FP(2) +
     &                                       A1(3,K)*FP(3)
              HD2 = HD2 + (H(IPAIR)*UI(K) + FREG(K))*A(K) +
     &                                      2.0D0*UIDOT(K)*AK4
              PRD(K) = 0.5D0*RI*AK4 + TD2*A(K)
              FRD(K) = HD*UI(K) + 0.5D0*DH*UIDOT(K) + PRD(K)
*       Form the regularized perturbation and modified force.
              PR(K) = 0.5D0*PR(K)
              FREG(K) = 0.5D0*FREG(K)
   15     CONTINUE
*
*       Determine new derivatives of U evaluated at beginning of the step.
          DO 20 K = 1,4
              DF = FP0(K,IPAIR) - FREG(K)
              SUM = FD0(K,IPAIR) + FRD(K)
	      BT2 = -3.0D0*DF - (SUM + FD0(K,IPAIR))*DTU
	      AT3 = 2.0D0*DF + SUM*DTU
*
              U2(K) = 0.5D0*H0(IPAIR)*U0(K,IPAIR) + FP0(K,IPAIR)
              U3(K) = 0.5D0*H0(IPAIR)*UDOT(K,IPAIR) + FD0(K,IPAIR)
              U4(K) = 0.5D0*H0(IPAIR)*U2(K) + BT2*DT2
              U5(K) = 0.5D0*H0(IPAIR)*U3(K) + AT3*DT6
   20     CONTINUE
*
*       Improve the solution of U and UDOT.
          DO 25 K = 1,4
              UI(K) = ((((U5(K)*DT5 + U4(K)*Z4)*DT4 + U3(K))*DT3 +
     &                 0.5*U2(K))*DTU + UDOT(K,IPAIR))*DTU + U0(K,IPAIR)
              UIDOT(K) = (((U5(K)*DTZ + U4(K)*Z3)*DT3 +
     &                       0.5*U3(K))*DTU + U2(K))*DTU + UDOT(K,IPAIR)
   25     CONTINUE
*
*       Update the physical distance for next iteration or final solution.
          RI = UI(1)**2 + UI(2)**2 + UI(3)**2 + UI(4)**2
*
*       Choose between old and new perturbation (skip last time).
*         IF (ITER.GT.ITMAX) THEN
*         IF (GAMMA(IPAIR).LT.1.0D-04) THEN
*             DO 30 K = 1,3
*                 FD(K) = RI*FD1(K)
*  30         CONTINUE
*         ELSE
*       Transform to improved coordinates & velocities.
*             CALL KSTRAN(IPAIR,I1,I,BODYIN,RI,UI,UIDOT,XI,VI)
*
*       Re-calculate the perturbing force & derivative.
*             CALL KSPERT(I1,NNB0,XI,VI,FP,FD)
*
*       Convert to regularized derivative (note slow-down not active).
*             DO 35  K = 1,3
*                 FD(K) = RI*FD(K)
*  35         CONTINUE
*         END IF
*         END IF
*
*       Re-evaluate DH by adding Hermite corrector.
          HD = 2.0D0*HD
          DHD = HDOT(IPAIR) - HD
          SUM = HDOT2(IPAIR) + HD2
          BT2 = -3.0D0*DHD - (SUM + HDOT2(IPAIR))*DTU
          AT3 = 2.0D0*DHD + SUM*DTU
          DH = DH0 + (0.25D0*AT3 + ONE3*BT2)*DTU
   40 CONTINUE
*
*       Copy final values and set higher derivatives.
      H(IPAIR) = H(IPAIR) + DH
      DO 50 K = 1,4
          U(K,IPAIR) = UI(K)
          U0(K,IPAIR) = UI(K)
          UDOT(K,IPAIR) = UIDOT(K)
          FUDOT2(K,IPAIR) = U4(K) + U5(K)*DTU
          FUDOT3(K,IPAIR) = U5(K)
          AK4 = A(K+4) + A1(1,K)*FP(1) + A1(2,K)*FP(2) + A1(3,K)*FP(3)
          PR(K) = 0.5D0*RI*A(K)
          PRD(K) = 0.5D0*RI*AK4 + TD2*A(K) + 0.5D0*HD*UI(K)
          FP0(K,IPAIR) = PR(K)
          FD0(K,IPAIR) = PRD(K)
   50 CONTINUE
*
*       Save new derivatives of H.
      HDOT(IPAIR) = HD
      HDOT2(IPAIR) = HD2
      HDOT3(IPAIR) = (3.0D0*AT3 + BT2)*DT2
      HDOT4(IPAIR) = AT3*DT6
*
*       Set new FU/2 & FUDOT/6 and form scalar terms for time derivatives.
      TD2 = 0.0D0
      TD3 = 0.0D0
      TDOT4 = 0.0D0
      TDOT5 = 0.0D0
      TDOT6 = 0.0D0
      H2 = 0.5D0*H(IPAIR)
      DO 60 K = 1,4
          U2K = PR(K) + H2*UI(K)
          U3K = PRD(K) + H2*UIDOT(K)
          FU(K,IPAIR) = 0.5D0*U2K
          FUDOT(K,IPAIR) = ONE6*U3K
          TD2 = TD2 + UI(K)*UIDOT(K)
          TD3 = TD3 + UIDOT(K)**2 + UI(K)*U2K
          TDOT4 = TDOT4 + UI(K)*U3K + 3.0D0*UIDOT(K)*U2K
          TDOT5 = TDOT5 + 0.5D0*FUDOT2(K,IPAIR)*UI(K) +
     &                    2.0D0*U3K*UIDOT(K) + 1.5D0*U2K**2
          TDOT6 = TDOT6 + U5(K)*UI(K) + 5.0D0*U4(K)*UIDOT(K) +
     &                    10.0D0*U3K*U2K
*       Note that TDOT4/2, TDOT5/4 and TDOT6/2 are calculated.
   60 CONTINUE
*
*       Save distance and second & third time derivatives.
      R(IPAIR) = RI
      TDOT2(IPAIR) = 2.0D0*TD2
      TDOT3(IPAIR) = 2.0D0*TD3
*
      RETURN
*
      END