brake3.f

Go to the documentation of this file.

      SUBROUTINE brake3(IPAIR,ITERM)
*
*
*       Gravitational radiation of hierarchical BH/NS 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)
      parameter(aursun=214.95d0)
      REAL*8  m1,m2,jorb
*
*
*       Define scalars and determine merger index.
      i = n + ipair
      i1 = 2*ipair - 1
      iterm = 0
      im = 0
      DO 1 k = 1,nmerge
         IF (namem(k).EQ.name(i)) im = k
    1 CONTINUE
*
*       Find location of secondary (i.e. ghost) and form mass factors & SEMI.
      CALL findj(i,i2,im)
      IF (i2.LT.0) go to 30
      IF(kstar(i2).LT.10)THEN
         iterm = 1
         go to 30
      ENDIF
      m1 = cm(1,im)*zmbar
      m2 = cm(2,im)*zmbar
      zmb = cm(1,im) + cm(2,im)
      zmu = cm(1,im)*cm(2,im)/zmb
      semi = -0.5*zmb/hm(im)
*
*       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
*
*       Form inner eccentricity and angular momentum.
      sep = semi*su
      ecc2 = (1.0 - rb/semi)**2 + td2**2/(zmb*semi)
      ecc0 = sqrt(ecc2)
      oorb = twopi*sqrt((m1 + m2)/(sep/aursun)**3)
      jorb = m1*m2/(m1 + m2)*sqrt(1.d0-ecc2)*sep**2*oorb
*
*       Obtain angular momentum and eccentricity derivatives (SEMI < 10 R_s).
      CALL grrad(m1,m2,sep,ecc0,jorb,djgr,delet)
*
*       Allow 2% angular momentum change (expressions adapted from MDOT).
      IF (abs(djgr).GT.0.d0) THEN
          dtgr = 0.02d0*jorb/abs(djgr)
          IF(delet.GT.tiny.AND.ecc0.GT.0.0011d0)THEN
             dtgr = min(dtgr,0.05d0*ecc0/delet)
          ENDIF
          dtgr = max(dtgr,100.d0)
*       Note that time interval is in units of yr (cf. grrad.f).
          djorb = djgr*dtgr
*       Terminate on large change of JORB which occurs at GR coalescence.
          IF (djorb.GT.0.1*jorb) THEN
              iterm = 1
              go to 30
          END IF
          jorb = jorb - djorb
          jorb = max(jorb,1.d0)
          ecc = max(ecc0 - delet*dtgr,0.001d0)
*       Convert to N-body units and update both TEV to same time (NB!).
          dtgr = dtgr/(1.0d+06*tstar)
          tev0(i1) = time
          tev0(i2) = time
          tev(i1) = time + dtgr
          tev(i2) = time + dtgr
      ELSE
          dt = min(tev(i1)-tev0(i1),tev(i2)-tev0(i2))
          dt = max(dt,0.1d0)
          tev0(i1) = time
          tev0(i2) = time
          tev(i1) = time + dt
          tev(i2) = time + dt
          go to 30
      END IF
*
*       Obtain new semi-major axis (in N-body units) from angular momentum.
      semi1 = (m1 + m2)*jorb*jorb/
     &                  ((m1*m2*twopi)**2*aursun**3*(1.d0-ecc**2))
      semi1 = semi1/su
*
*       Update inner 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
*
*       Specify new peri- or apo-centre at turning-point.
      IF (rb.LT.semi1) THEN
          rnew = semi1*(1.0 - ecc)
          efac = (1.0 - ecc)/(1.0 - ecc0)
      ELSE
          rnew = semi1*(1.0 + ecc)
          efac = (1.0 + ecc)/(1.0 + ecc0)
      END IF
      c1 = sqrt(efac)
      v2 = 0.5*(zmb + hm(im)*rnew)
      IF (v2.LE.0.d0) THEN
         c2 = 1.0d-06
      ELSE
         c2 = sqrt(v2/v20)
      END IF
*
*       Re-scale KS variables at constant energy to new eccentricity (ECC).
      rb = 0.0d0
      DO 20 k = 1,4
          um(k,im) = c1*um(k,im)
          umdot(k,im) = c2*umdot(k,im)
          rb = rb + um(k,im)**2
   20 CONTINUE
*
*       Rectify the hierarchical KS variables.
      CALL hirect(im)
*
*       Include GR coalescence criterion for compact objects.
      ksx = max(kstar(i1),kstar(i2))
      IF (ksx.GE.13.AND.kz(28).GT.0) THEN
         rcoal = 6.0*zmb/clight**2
      ELSE
         rcoal = 0.0
      END IF
*
      WRITE (6,25)  name(i1), toff+time, ecc, semi1, dtgr, rcoal
   25 FORMAT (' BRAKE3    NM T E A DTGR RCOAL',i8,f9.2,f8.4,1p,3e10.2)
*
*       Check termination criterion for coalescence.
      IF (rb.LT.rcoal) THEN
          iterm = 1
      END IF
   30 CONTINUE
*
      RETURN
*
      END