hmdot2.f

Go to the documentation of this file.

      SUBROUTINE hmdot2(J,IGHOST,IMERGE,M1,KW,MC,DMS,RNEW,ITERM)
*
*
*       Mass loss from outer 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)
      REAL*8  bodyi(2),w(2)
      REAL*8  m1,mc
      LOGICAL iquit
*
*
*       Define merger termination index and relevant KS/c.m. indices.
      iterm = 0
      ipair = kspair
      icm = n + ipair
      jpair = ighost - n
      i1 = 2*jpair - 1
      i2 = i1 + 1
*
*       Include quitting conditions to avoid large changes at end-point.
      iquit = .false.
*     IF (RADIUS(J)*SU.GT.500.0) IQUIT = .TRUE.
*     IF (BODY0(J)*ZMBAR.GT.15.0.AND.KSTAR(J).GE.4) IQUIT = .TRUE.
*     IF (ABS(M1 - MC).LT.0.1*M1.OR.MC.GT.1.4) IQUIT = .TRUE.
*     IF (BODY0(J)*ZMBAR - M1.GT.0.5*M1) IQUIT = .TRUE.
*       Note that tidal dissipation needs full KS for rectification in MDOT.
*     IF (KSTAR(IGHOST).LT.0.OR.BODY(ICM).EQ.0.0D0) IQUIT = .TRUE.
      IF (body(icm).EQ.0.0d0) iquit = .true.
*
*       Quit for mis-identification or advanced evolution.
      IF (ighost.LE.n.OR.iquit) THEN
          IF (ighost.LE.n) THEN
              WRITE (6,5)  name(i1), ighost
    5         FORMAT (' DANGER!    HMDOT2    NAM IG ',2i5)
              stop
          END IF
          iterm = 1
          go to 50
      END IF
*
*       Update radius for #J and copy #I2 for MDOT & HCORR in case J = I2.
      rstar = rnew
      IF (j.GT.i1) THEN
          radius(j) = rnew
          rnew = radius(i2)
*       Return RNEW for copying back to RADIUS(I2) on DM/M < 0.01.
      END IF
*
*       Skip further modifications on zero mass loss of same type.
      IF (abs(dms/m1).EQ.0.0d0.AND.kw.EQ.kstar(j)) go to 50
*
*       Set two-body elements for outer orbit and inner semi-major axis.
      semi = -0.5*body(icm)/h(ipair)
      ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(icm)*semi)
      ecc1 = sqrt(ecc2)
      zmb0 = cm(3,imerge) + cm(4,imerge)
*       Note old H(JPAIR) used with new inner SEMI0 (small effect).
      semi0 = -0.5*zmb0/h(jpair)
*
*       Form modified mass ratio and semi-major axes from M*A = const.
      dm = dms/zmbar
      zmb = zmb0 - dm
      semi1 = semi*(zmb0 + body(2*ipair-1))/(zmb + body(2*ipair-1))
      semi2 = semi0*zmb0/zmb
      pmin = semi1*(1.0 - ecc1)
*
*       Evaluate old separation, square regularized velocity, t'' & ECC.
      ri = 0.0d0
      v20 = 0.0
      td2 = 0.0
      DO 10 k = 1,4
          ri = ri + u(k,jpair)**2
          v20 = v20 + udot(k,jpair)**2
          td2 = td2 + 2.0*u(k,jpair)*udot(k,jpair)
   10 CONTINUE
      ecc2 = (1.0 - ri/semi0)**2 + td2**2/(semi0*zmb0)
      ecc = sqrt(ecc2)
*
*       Obtain stability parameters for the new configuration.
      IF (ecc1.LT.1.0) THEN
          nst = nstab(semi2,semi1,ecc,ecc1,zi,cm(1,imerge),
     &                                     cm(2,imerge),body(2*ipair))
          IF (nst.EQ.0) THEN
              pcrit = 0.99*pmin
          ELSE
              pcrit = 1.01*pmin
          END IF
      ELSE
          pcrit = stability(cm(3,imerge),cm(4,imerge),body(2*ipair-1),
     &                                          ecc,ecc1,0.0d0)*semi2
      END IF
*
*       Update pericentre distance on successful stability test or exit.
      IF (pmin.GT.pcrit.AND.h(jpair).LT.-eclose) THEN
          r0(ipair) = pcrit
      ELSE
          iterm = 1
          go to 50
      END IF
*
*       Check for possible circularization (tidal or sequential).
      rp = semi2*(1.0 - ecc)
      IF (kz(27).GT.0.AND.rp.LT.4.0*radius(j).AND.ecc.GT.0.002) THEN
          tc = 1.0d+10
          IF (kz(27).EQ.2) THEN
              j1 = 2*jpair - 1
              j2 = j1 + 1
              DO 14 k = 1,2
                  bodyi(k) = cm(k+2,imerge)
   14         CONTINUE
              tg = 1.0
*       Evaluate the circularization time.
              CALL hicirc(rp,ecc,j1,j2,bodyi,tg,tc,ec,ed,w)
              IF (tc.GT.2000) go to 18
          END IF
          WRITE (6,15)  name(j), kstar(j), kstarm(imerge), ecc, dms,
     &                  rp, radius(j), tc
   15     FORMAT (' HMDOT2 TERM    NM K* E DM RP R* TC ',
     &                             i6,2i4,f8.4,f7.3,1p,3e10.2)
          iterm = 1
          go to 50
      END IF
*
*       Obtain energy change from M*A = const and H = -M/(2*A) (DCH 8/96).
   18 dh = dm/semi0*(1.0 - 0.5*dm/zmb0)
      hi = h(jpair)
      h(jpair) = h(jpair) + dh
*
*       Form KS coordinate & velocity scaling factors (general point is OK).
      semi2 = -0.5*zmb/h(jpair)
      c2 = sqrt(semi2/semi0)
      v2 = 0.5*(zmb + h(jpair)*ri*(semi2/semi0))
      c1 = sqrt(v2/v20)
*
*       Re-scale KS variables to new energy (H < 0: constant eccentricity).
      r(jpair) = 0.0d0
      DO 20 k = 1,4
          u(k,jpair) = c2*u(k,jpair)
          udot(k,jpair) = c1*udot(k,jpair)
          r(jpair) = r(jpair) + u(k,jpair)**2
   20 CONTINUE
*
*       Reduce mass of relevant component (ghost is original second member).
      zmu0 = cm(3,imerge)*cm(4,imerge)/zmb0
      km = 3
      IF (j.EQ.2*jpair) km = 4
      cm(km,imerge) = cm(km,imerge) - dm
*
*       Include corrections to EMERGE & EMDOT (consistency; no net effect!).
      zmu = cm(3,imerge)*cm(4,imerge)/zmb
      decorr = zmu*h(jpair) - zmu0*hi
      emerge = emerge + decorr
      emdot = emdot - decorr
      egrav = egrav - decorr
*
*       Correct outer orbit for mass loss (use #I2 for consistency).
      CALL hcorr(i2,dm,rstar)
*
*       Print some diagnostics on non-zero mass loss.
      IF (dms.GT.1.0d-03) THEN
          WRITE (6,30)  name(j), kstar(j), imerge, km, m1, dms, pmin,
     &                  pcrit, semi2, semi1, r(jpair), h(jpair), decorr
   30     FORMAT (' HMDOT2    NAM K* IM KM M1 DM PM PC A A1 R H DE ',
     &                        i6,3i4,f6.2,1p,7e10.2,0p,f10.6)
      END IF
*
   50 RETURN
*
      END