hcorr.f

Go to the documentation of this file.

      SUBROUTINE hcorr(I,DM,RNEW)
*
*
*       Mass loss correction of KS orbit.
*       ---------------------------------
*
      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),iflag(mmax)
*
*
*       Copy pair index and determine secondary component (note I1 = I here).
      ipair = kspair
      i1 = i
      i2 = 2*ipair - 1
      IF (i2.EQ.i1) THEN
          i2 = i2 + 1
      END IF
*
*       Define c.m. index and old semi-major axis.
      j = n + ipair
      semi0 = -0.5d0*body(j)/h(ipair)
*
*       Distinguish different cases (KS, triple or outer ghost of quadruple).
      IF (body(i).GT.0.0d0) THEN
          zmu0 = body(i1)*body(i2)/body(j)
      ELSE
          im = 0
          jj = i
          IF (i.LT.ifirst) jj = n + kvec(i)
*       Determine merger index for inner and outer binary masses.
          DO 1 k = 1,nmerge
              IF (nameg(k).EQ.name(jj)) im = k
    1     CONTINUE
          IF (im.EQ.0) THEN
              WRITE (6,5) name(i), name(jj), (nameg(k),k=1,nmerge)
    5         FORMAT (' DANGER HCORR!    NMI NMJJ NMG ',12i6)
              stop
          END IF
*       Form old reduced mass from inner and outer binary.
          zm1 = cm(1,im) + cm(2,im)
          zm2 = cm(3,im) + cm(4,im)
          zmu0 = zm1*zm2/body(j)
      END IF
*
*       Obtain energy change from MK*A = const and H = -M/(2*A) (DCH 8/96).
      dh = dm/semi0*(1.0 - 0.5*dm/body(j))
*
*       Reduce mass of c.m. body and subtract energy correction terms.
      body(j) = body(j) - dm
      IF (body(i).GT.0.0d0) THEN
          zmu1 = (body(i1) - dm)*body(i2)/body(j)
      ELSE
          zmu1 = zm1*(zm2 - dm)/body(j)
      END IF
      emdot = emdot - zmu1*dh - (zmu1 - zmu0)*h(ipair)
*       Retain binding energy gained in past interactions.
      egrav = egrav - zmu1*dh - (zmu1 - zmu0)*h(ipair)
*
*       Include diagnostic warning for large relative energy change.
      IF (dh.GT.0.2*abs(h(ipair))) THEN
          WRITE (6,10)  name(i), dh, h(ipair), r(ipair)/semi0, dm*smu
   10     FORMAT (' WARNING!    LARGE CORRECTION    NM DH H R/A DMS ',
     &                                              i6,2f8.2,2f7.2)
      END IF
*
*       Update the binding energy due to mass loss DM and set new radius.
      h(ipair) = h(ipair) + dh
      radius(i) = rnew
*
*       Modify KS variables at constant eccentricity.
      CALL expand(ipair,semi0)
*
      RETURN
*
      END