ksrect.f

Go to the documentation of this file.

      SUBROUTINE ksrect(IPAIR)
*
*
*       Rectification of KS orbit.
*       --------------------------
*
      include 'common6.h'
*
*
*       Include some diagnostic output.
      upr2 = 0.0
      DO 1 k = 1,4
          upr2 = upr2 + udot(k,ipair)**2
    1 CONTINUE
*
*       Skip rectification for small eccentricity or large perturbation.
      i = n + ipair
      semi = -0.5*body(i)/h(ipair)
      ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(semi*body(i))
      ecc = sqrt(ecc2)
      IF (ecc.LE.0.01) go to 50
      IF (gamma(ipair).GT.0.1) go to 50
*
*       Include Roche rectification for large ECC (additional to roche.f).
      IF (kstar(i).GT.10.AND.ecc.GT.0.2) THEN
          qp = semi*(1.0 - ecc)
          icirc = -1
          CALL tcirc(qp,ecc,2*ipair-1,2*ipair,icirc,tc)
          WRITE (6,2)  name(2*ipair-1), kstar(i), ecc, qp, tc,
     &                 gamma(ipair)
    2     FORMAT (' RECTIFY K*    NM K* E QP TC G',i8,i4,f8.4,1p,3e10.2)
          kstar(i) = 0
      END IF
*
      hi = (2.0*upr2 - body(n+ipair))/r(ipair)
      err = (hi - h(ipair))/hi
      zmu = body(2*ipair)*body(2*ipair-1)/body(n+ipair)
      db = zmu*(hi - h(ipair))
      IF (abs(db).GT.1.0d-08) THEN
          semi = -0.5*body(n+ipair)/h(ipair)
          ra = r(ipair)/semi
          IF (semi.LT.0.0) ra = r(ipair)
          WRITE (16,3)  time+toff, ipair, ra, h(ipair), gamma(ipair),
     &                  db, err
    3     FORMAT (' KSRECT:   T # R/A H G DB DH/H ',
     &                        f8.2,i4,f8.4,f8.1,f7.3,1p,2e10.1)
          CALL flush(16)
      END IF
*
*       Initialize iteration counter for difficult case (SJA 10/97).
      iter = 0
*
*       Form square regularized velocity for the explicit binding energy.
   10 upr2 = 0.0
      DO 15 k = 1,4
          upr2 = upr2 + udot(k,ipair)**2
   15 CONTINUE
*
*       Form KS scaling factors from energy and angular momentum relation.
      a1 = 0.25d0*body(n+ipair)/upr2
*       Solve for C1 from H = (2*U'*U'*C1**2 - M)/(U*U*C2**2) with C2 = 1/C1.
      a2 = a1**2 + 0.5d0*h(ipair)*r(ipair)/upr2
*
*       Avoid negative round-off value on second try (NB! no change in CK).
      IF (iter.EQ.2.AND.a2.LT.0.0) a2 = 0.0d0
*
*       Check for undefined case (circular orbit or eccentric anomaly = 90).
      IF (a2.GE.0.0d0) THEN
          IF (a1.LT.1.0) THEN
*       Choose square root sign from eccentric anomaly (e*cos(E) = 1 - R/a).
              c1 = sqrt(a1 + sqrt(a2))
          ELSE
              c1 = sqrt(a1 - sqrt(a2))
          END IF
          ck = 1.0
      ELSE
*       Adopt C1*C2 = CK for difficult case (Seppo's suggestion of 1991).
          c1 = 1.0
          ck = body(n+ipair)/sqrt(-8.0d0*h(ipair)*r(ipair)*upr2)
          WRITE (6,20)  ipair, kstar(n+ipair), ecc, r(ipair), h(ipair),
     &                  gamma(ipair), upr2, a2, ck-1.0
   20     FORMAT (' WARNING!    KSRECT    KS K* E R H G UPR2 A2 CK-1 ',
     &                                    2i4,f8.4,1p,6e10.2)
          iter = iter + 1
      END IF
*
*       Specify KS coordinate scaling from angular momentum conservation.
      c2 = ck/c1
*
*       Transform KS variables to yield the prescribed elements.
      r(ipair) = 0.0d0
*     UPR2 = 0.0D0
      td2 = 0.0d0
      DO 25 k = 1,4
          u(k,ipair) = c2*u(k,ipair)
          udot(k,ipair) = c1*udot(k,ipair)
          u0(k,ipair) = u(k,ipair)
          r(ipair) = r(ipair) + u(k,ipair)**2
          td2 = td2 + u(k,ipair)*udot(k,ipair)
*         UPR2 = UPR2 + UDOT(K,IPAIR)**2
   25 CONTINUE
      tdot2(ipair) = 2.0*td2
*
*       Include diagnostic output.
*     HI = (2.0*UPR2 - BODY(N+IPAIR))/R(IPAIR)
*     WRITE (16,30)  IPAIR, IPHASE, KSTAR(N+IPAIR), R(IPAIR),
*    &               GAMMA(IPAIR), (HI-H(IPAIR))/HI
*  30 FORMAT (' KSRECT:    KS IPH K* R G DH/H ',3I4,1P,3E10.2)
*     CALL FLUSH(16)
*
*       Improve solution by second iteration in case of CK procedure.
      iter = iter + 1
      IF (iter.EQ.2) go to 10
*
   50 RETURN
*
      END