hmdot.f

Go to the documentation of this file.

      SUBROUTINE hmdot(J,IMERGE,M1,KW,MC,DMS,RNEW,ITERM)
*
*
*       Mass loss from inner 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
      i1 = 2*ipair - 1
      icm = n + ipair
*
*       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 (KSTARM(IMERGE).LT.0.OR.BODY(ICM).EQ.0.0D0) IQUIT = .TRUE.
      IF (body(icm).EQ.0.0d0) iquit = .true.
*
*       Quit for mis-identification, advanced evolution or double merger.
      IF (j.LE.0.OR.iquit) THEN
          iterm = 1
          go to 50
      END IF
*
*       Update radius for #J and copy #I1 for MDOT & HCORR in case of ghost.
      rstar = rnew
      IF (j.GT.i1) THEN
          radius(j) = rnew
          rstar = radius(i1)
*       Return RNEW for copying back to RADIUS(I1) on ghost & 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(1,imerge) + cm(2,imerge)
      semi0 = -0.5*zmb0/hm(imerge)
*
*       Form modified mass ratio and semi-major axes from M*A = const.
      dm = dms/zmbar
      zmb = zmb0 - dm
      semi1 = semi*(zmb0 + body(2*ipair))/(zmb + body(2*ipair))
      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 + um(k,imerge)**2
          v20 = v20 + umdot(k,imerge)**2
          td2 = td2 + 2.0*um(k,imerge)*umdot(k,imerge)
   10 CONTINUE
      ecc2 = (1.0 - ri/semi0)**2 + td2**2/(semi0*zmb0)
      ecc = sqrt(ecc2)
*
*       Determine inclination (use #I1 as first KS component).
      CALL himax(i1,imerge,ecc,semi0,emax,emin,zi,tg,edav)
*
*       Evaluate the general stability function (Mardling 2008).
      IF (ecc1.LT.1.0) THEN
          eout = ecc1
*       Increase tolerance near sensitive stability boundary (RM 10/2008).
          IF (eout.GT.0.9) THEN
              de = 0.5*(1.0 - eout)
              de = min(de,0.01d0)
*       Add extra amount 0.011 to avoid switching.
              IF (ecc1.GT.0.9) de = de + 0.011
              eout = eout - de
              pmin = semi1*(1.0 - eout)
          END IF
          nst = nstab(semi2,semi1,ecc,eout,zi,cm(1,imerge),
     &                                     cm(2,imerge),body(2*ipair))
          IF (nst.EQ.0) THEN
              pcrit = 0.99*pmin
              pcr = stability(cm(1,imerge),cm(2,imerge),body(2*ipair),
     &                                          ecc,ecc1,zi)*semi2
          ELSE
              pcrit = 1.01*pmin
          END IF
      ELSE
          pcrit = stability(cm(1,imerge),cm(2,imerge),body(2*ipair),
     &                                          ecc,ecc1,zi)*semi2
      END IF
*
*       Update pericentre distance on successful stability test or exit.
      IF (pmin.GT.pcrit) 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.AND.
     &    kstarm(imerge).GE.0) THEN
          tc = 1.0d+10
          IF (kz(27).EQ.2) THEN
              i2 = 0
*       Identify the ghost by searching single bodies.
              DO 12 k = ifirst,n
                  IF (body(k).EQ.0.0d0.AND.
     &                name(k).EQ.nameg(imerge)) THEN
                      i2 = k
                  END IF
   12         CONTINUE
              IF (i2.GT.0) THEN
                  DO 14 k = 1,2
                      bodyi(k) = cm(k,imerge)
   14             CONTINUE
                  tg = 1.0
*       Evaluate the circularization time.
                  CALL hicirc(rp,ecc,i1,i2,bodyi,tg,tc,ec,ed,w)
                  IF (tc.GT.2000) go to 18
              END IF
          END IF
          WRITE (6,15)  name(j), kstar(j), kstarm(imerge), ecc, dms,
     &                  rp, radius(j), tc
   15     FORMAT (' HMDOT TERM    NAM K* E DM RP R* TC ',
     &                            i6,2i4,f8.4,f7.3,1p,3e10.2)
          iterm = 1
          go to 50
      ELSE IF (ecc.GT.0.002) THEN
          IF (rp.LT.2.0*radius(j)) THEN
              WRITE (6,15)  name(j), kstar(j), kstarm(imerge), ecc, dms,
     &                      rp, radius(j)
              iterm = 1
              go to 50
          END IF
      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)
      hm0 = hm(imerge)
      hm(imerge) = hm(imerge) + dh
*
*       Form KS coordinate & velocity scaling factors (general point is OK).
      semi2 = -0.5*zmb/hm(imerge)
      c2 = sqrt(semi2/semi0)
      v2 = 0.5*(zmb + hm(imerge)*ri*(semi2/semi0))
      c1 = sqrt(v2/v20)
*
*       Re-scale KS variables to new energy (H < 0: constant eccentricity).
      DO 20 k = 1,4
          um(k,imerge) = c2*um(k,imerge)
          umdot(k,imerge) = c1*umdot(k,imerge)
*       Note no need to update XREL & VREL for RESTART (c.m. error cancels).
   20 CONTINUE
*
*       Reduce mass of relevant component (ghost is original second member).
      zmu0 = cm(1,imerge)*cm(2,imerge)/zmb0
      km = 1
      IF (j.GE.ifirst) km = 2
      cm(km,imerge) = cm(km,imerge) - dm
*
*       Include corrections to EMERGE & EMDOT (consistency; no net effect!).
      zmu = cm(1,imerge)*cm(2,imerge)/zmb
      decorr = zmu*hm(imerge) - zmu0*hm0
      emerge = emerge + decorr
      emdot = emdot - decorr
*       Compensate for EMERGE contribution in DEGRAV (cf. EVENTS).
      egrav = egrav - decorr
*
*       Correct outer orbit for mass loss (use #I1 in case #J is a ghost).
      CALL hcorr(i1,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, decorr
   30     FORMAT (' HMDOT    NAM K* IM KM M1 DM PM PC A A1 DE ',
     &                       i6,3i4,f6.2,1p,5e10.2,0p,f10.6)
      END IF
*
   50 RETURN
*
      END