hirect.f

Go to the documentation of this file.

      SUBROUTINE hirect(IM)
*
*
*       Rectification of hierarchical binary.
*       -------------------------------------
*
      include 'common6.h'
      common/binary/  cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
     &                hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
     &                namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
*
*
*       Include diagnostic error check.
      rb = 0.0
      upr2 = 0.0
      td2 = 0.0
      DO 5 k = 1,4
          rb = rb + um(k,im)**2
          upr2 = upr2 + umdot(k,im)**2
          td2 = td2 + 2.0*um(k,im)*umdot(k,im)
    5 CONTINUE
*
      zmb = cm(1,im) + cm(2,im)
      hi = (2.0*upr2 - zmb)/rb
      err = (hi - hm(im))/hi
      zmu = cm(1,im)*cm(2,im)/zmb
      db = zmu*(hi - hm(im))
      IF (abs(db).GT.1.0d-08) THEN
      semi = -0.5*zmb/hm(im)
      ecc2 = (1.0 - rb/semi)**2 + td2**2/(zmb*semi)
      ecc = sqrt(ecc2)
      ra = rb/semi
      WRITE (16,3) ttot, nameg(im), kstarm(im), ecc, ra, hm(im), db, err
    3 FORMAT (' HIRECT:   T NM K* E R/A H DB DH/H ',
     &                    f9.3,i6,i4,f8.4,f8.4,f9.2,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 rb = 0.0d0
      upr2 = 0.0
      DO 15 k = 1,4
          rb = rb + um(k,im)**2
          upr2 = upr2 + umdot(k,im)**2
   15 CONTINUE
*
*       Form KS scaling factors from energy and angular momentum relation.
      a1 = 0.25d0*zmb/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*hm(im)*rb/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 = zmb/sqrt(-8.0d0*hm(im)*rb*upr2)
          WRITE (6,20)  im, kstarm(im), rb, hm(im), upr2, a2, ck-1.0
   20     FORMAT (' WARNING!    HIRECT    IM K* R H UPR2 A2 CK-1 ',
     &                                    2i4,1p,5e10.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.
      DO 25 k = 1,4
          um(k,im) = c2*um(k,im)
          umdot(k,im) = c1*umdot(k,im)
   25 CONTINUE
*
*       Improve solution by second iteration in case of CK procedure.
      iter = iter + 1
      IF (iter.EQ.2) go to 10
*
      RETURN
*
      END