brake2.f

Go to the documentation of this file.

      SUBROUTINE brake2(IPAIR,ITERM)
*
*
*       Gravitational radiation 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)
      common/slow0/  range,islow(10)
      REAL*8 m0,aj,r1,lum1,mc,rcc,menv,renv,k2
      REAL*8 tscls(20),lums(10),gb(10),tm,tn
      REAL*8  m1,m2,rl,q
      EXTERNAL  rl,rtmsf
*
*
      i = n + ipair
      i1 = 2*ipair - 1
      im = 0
      DO 1 k = 1,nmerge
         IF (namem(k).EQ.name(i)) im = k
    1 CONTINUE
*       Find location of the secondary (i.e. ghost).
      CALL findj(i,i2,im)
      m1 = cm(1,im)*zmbar
      m2 = cm(2,im)*zmbar
*
*       Quit if active ROCHE or mass-losing star.
      iterm = 1
      IF(kstarm(im).EQ.11.OR.kstarm(im).EQ.13)THEN
          go to 30
      END IF
      IF (m1.GE.10.0.OR.(kstar(i1).GT.1.AND.kstar(i1).LT.10))THEN
          go to 30
      END IF
      IF (m2.GE.10.0.OR.(kstar(i2).GT.1.AND.kstar(i2).LT.10))THEN
          go to 30
      END IF
*
*       Skip braking if semi-major axis exceeds 10 solar radii.
      iterm = 0
      zmb = cm(1,im) + cm(2,im)
      semi = -0.5*zmb/hm(im)
      IF (semi*su.GT.10.0)THEN
          tgr = 1.0d+06
          tbr = 1.0d+06
          semi1 = semi
          go to 20
      END IF
*
*       Obtain current separation, square velocity and t'' = 2*U*U'.
      rb = 0.0
      v20 = 0.0
      td2 = 0.0
      DO 5 k = 1,4
          rb = rb + um(k,im)**2
          v20 = v20 + umdot(k,im)**2
          td2 = td2 + 2.0*um(k,im)*umdot(k,im)
    5 CONTINUE
*
*       Evaluate inner eccentricity and outer period.
      ecc2 = (1.0 - rb/semi)**2 + td2**2/(zmb*semi)
      ecc = sqrt(ecc2)
      semi0 = -0.5*body(i)/h(ipair)
*
*       Identify the most magnetically active star (#J2).
      IF (radius(i1).GE.radius(i2)) THEN
          IF ((kstar(i1).EQ.1.AND.m1.GT.1.25).OR.
     &        (kstar(i1).EQ.0.AND.m1.LT.0.35)) THEN
              j2 = i2
          ELSE
              j2 = i1
          END IF
      ELSE
          IF ((kstar(i2).EQ.1.AND.m2.GT.1.25).OR.
     &        (kstar(i2).EQ.0.AND.m2.LT.0.35)) THEN
              j2 = i1
          ELSE
              j2 = i2
          END IF
      END IF
*
*       Set mass and radius of star #J2 in solar units.
      IF (j2.EQ.i1) m2 = m1
      r2 = radius(j2)*su
*
*       Define binary period, solar radius & GMS in cgs units.
      tb = 3.147d+07*yrs*semi*sqrt(semi/zmb)
      rsun = 6.96d+10
      gms = gm*zmbar*zmb
*
*       Form time derivative of angular momentum (Regos & Tout 1995).
      zjdot = 3.8d-30*1.989d+33*m2*rsun**4*r2**3*(twopi/tb)**3
*
*       Determine new semi-major axis from angular momentum relation.
      zmu = cm(1,im)*cm(2,im)/zmb
      acm = 3.08d+18*rbar*semi
      adot = 2.0*sqrt(acm/gms)*zjdot/(1.989d+33*zmbar*zmu)
*
*       Define old primary and evaluate time scale for spin-down (yrs).
      tbr = acm/(3.147d+07*adot)
*       Evaluate timescale for 10% change in semi-major axis (N-body units). 
      tbr = 0.1d0*tbr/(1.0d+06*tstar)
*
*       Include gravitational radiation (64*GM**3/5*c**5; Ap.J. 418, 147).
*     FAC = 64.0*((6.67D-08*1.989D+33)/3.0D+10)**3/(5.0*9.0D+20)
      agdot = 1.23d+27*cm(1,im)*cm(2,im)*zmb*(zmbar/acm)**3
      tgr = acm/(3.147d+07*agdot)
      tgr = 0.1d0*tgr/(1.0d+06*tstar)
*
*       Suppress magnetic braking for massive MS/low-mass or evolved star.
      IF (((m2.GT.1.25.OR.m2.LT.0.35).AND.kstar(j2).LE.2).OR.
     &      kstar(j2).GE.10) THEN
          adot = 0.0
          tbr = 1.0d+06
      END IF
*
*       Specify time interval.
      dt = tev(i) - tev0(i)
      tk = twopi*abs(semi0)*sqrt(abs(semi0)/body(i))
*
*       Combine effects but include possible cutoff of magnetic braking.
      adot = adot + agdot
*
*       Convert from cgs to scaled units and update semi-major axis.
      adot = adot/(1.0d+05*vstar)
      IF (adot*dt.GT.0.1*semi) THEN
          semi1 = 0.9*semi
      ELSE
          semi1 = semi - adot*dt
      END IF
*
      WRITE (6,71)  semi, semi1, adot, dt
   71 FORMAT (' BRAKE2 !!    A A1 AD DT  ',1p,4e10.2)
*       Include safety test on new semi-major axis.
      IF (abs(semi1).LT.radius(j2).OR.semi1.LT.0.0) THEN
          semi1 = radius(j2)
      END IF
*
*       Update binding energy.
      hi = hm(im)
      hm(im) = -0.5*zmb/semi1
*
*       Correct EMERGE & ECOLL (consistency; no net effect).
      decorr = zmu*(hi - hm(im))
      emerge = emerge - decorr
      ecoll = ecoll + decorr
      egrav = egrav + decorr
*
*       Specify KS coordinate & velocity scaling factors at general point.
      c2 = sqrt(semi1/semi)
      v2 = 0.5*(zmb + hm(im)*rb*(semi1/semi))
      c1 = sqrt(v2/v20)
*
*       Re-scale KS variables to new energy with constant eccentricity.
      rb = 0.0d0
      DO 10 k = 1,4
          um(k,im) = c2*um(k,im)
          umdot(k,im) = c1*umdot(k,im)
          rb = rb + um(k,im)**2
   10 CONTINUE
*
*       Rectify the hierarchical KS variables.
      CALL hirect(im)
*
*       Determine indices for primary & secondary star (donor & accretor).
   20 j1 = i1
      j2 = i2
      q = cm(1,im)/cm(2,im)
      rl1 = rl(q)*semi1
*       Evaluate Roche radius for the second star.
      q = 1.0/q
      rl2 = rl(q)*semi1
*
*       Update radius for MS stars to current time. 
      tm0 = 100.d0
      rtms = 0.d0
      IF(kstar(i1).LE.1)THEN
         m0 = m1
         mc = 0.d0
         aj = time*tstar - epoch(i1)
         kw = kstar(i1)
         CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
         CALL hrdiag(m0,aj,m1,tm,tn,tscls,lums,gb,zpars,
     &               r1,lum1,kw,mc,rcc,menv,renv,k2)
         IF(kw.GT.1)THEN
            tm0 = 0.d0
            tev(i1) = min(tev(i1),time)
         ELSE
            tm = (tm + epoch(i1))/tstar - time
            tm0 = min(tm0,tm)
            radius(i1) = r1/su
            tev0(i1) = time
            rtms = rtmsf(m1)/su
            dtms = tm
         ENDIF
      ENDIF
      IF(kstar(i2).LE.1)THEN
         m0 = m2
         mc = 0.d0
         aj = time*tstar - epoch(i2)
         kw = kstar(i2)
         CALL star(kw,m0,m2,tm,tn,tscls,lums,gb,zpars)
         CALL hrdiag(m0,aj,m2,tm,tn,tscls,lums,gb,zpars,
     &               r1,lum1,kw,mc,rcc,menv,renv,k2)
         IF(kw.GT.1)THEN
            tm0 = 0.d0
            tev(i2) = min(tev(i2),time)
         ELSE
            tm = (tm + epoch(i2))/tstar - time
            tm0 = min(tm0,tm)
            radius(i2) = r1/su
            tev0(i2) = time
            rtms2 = rtmsf(m2)/su
            IF(rtms2.GT.rtms)THEN
               rtms = rtms2
               dtms = tm
            ENDIF
         ENDIF
      ENDIF
      IF(kstar(i1).GE.10.AND.tm0.GT.0.d0) tev0(i1) = time
      IF(kstar(i2).GE.10.AND.tm0.GT.0.d0) tev0(i2) = time
*
*       Compare scaled Roche radii when choosing the primary.
      IF (radius(j1)/rl1.LT.radius(j2)/rl2) THEN
          j1 = i2
          j2 = i1
          rl1 = rl2
      END IF
*
*       Check for Roche mass transfer at every stage (AM CVn formation).
      IF (radius(j1).GT.rl1.AND.kz(34).GT.0) THEN
          WRITE (6,25)  name(i1), nameg(im), kstar(i1), kstar(i2),
     &                  semi1*su, rl1*su, radius(j1)*su
   25     FORMAT (' BRAKE2 TERM    NM K* A RL R*  ',2i6,2i4,3f7.3)
          iterm = 1
*
*       Check for termination of MS stage. 
      ELSEIF (tm0.EQ.0.d0) THEN
          tev(i) = time + step(i1)
          iterm = 1
*
*       Increase TEV for standard case. 
      ELSE
          dt = min(tgr,tbr,tm0)
          IF(rtms.GT.rl1)THEN
*       Estimate time of Roche-lobe filling for MS star. 
             dt1 = dtms*(rl1 - radius(j1))/(rtms - radius(j1))
             dt = min(dt1,dt)
          ENDIF
          tev0(i) = tev(i)
          tev(i) = tev(i) + dt
          tev(i1) = tev(i) + step(i1)
          tev(i2) = tev(i1)
      END IF
*
   30 RETURN
*
      END