brake.f

Go to the documentation of this file.

      SUBROUTINE brake(IPAIR,DSEP,ECC1)
*
*
*       Orbital changes (GR, mass loss and/or tides).
*       ---------------------------------------------
*
      include 'common6.h'
      REAL*8  m1,m2
      SAVE ndiag
      DATA ndiag /0/
*
*
*       Set indices of KS components & c.m.
      i1 = 2*ipair - 1
      i2 = i1 + 1
      i = n + ipair
      dsep = dsep/su 
*
*       Set mass and radius of each star in solar units.
      m1 = body(i1)*zmbar
      m2 = body(i2)*zmbar
      r1 = radius(i1)*su
      r2 = radius(i2)*su
*
*       Obtain current semi-major axis & eccentricity and define new semi.
      semi = -0.5d0*body(i)/h(ipair)
      ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
      ecc = sqrt(ecc2)
      semi1 = semi - dsep 
*
*       Include circularization check (otherwise unperturbed KS missed).
      IF (kz(27).GT.0.AND.ecc.GT.0.01) THEN
          qperi = semi*(1.0 - ecc)
          icirc = -1
          CALL tcirc(qperi,ecc,i1,i2,icirc,tc)
          IF (tc.LT.100.0) THEN
              kx = max(kstar(i1),kstar(i2))
              WRITE (6,3)  icirc, kx, ecc, qperi*su, tc, max(r1,r2)
    3         FORMAT (' TCIRC CHECK   IC K* E QP TC R* ',2i5,1p,3e10.2)
              CALL flush(6)
          END IF
          IF (icirc.GT.0.AND.tc.LT.100.0) THEN
              CALL kstide(ipair,qperi)
              go to 50
          END IF
      END IF
*
*      Check for circularized case. 
      IF(kstar(i).LT.10.AND.ecc1.LT.0.002)THEN
         ecc1 = 0.001d0
         kstar(i) = 10
      ENDIF
*
*       Check collision/coalescence (IQCOLL=-2 skips ECOLL updating in COAL).
      rp = semi1*su*(1.d0 - ecc1)
*
*       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*body(i)/clight**2
      ELSE
         rcoal = r1 + r2
      END IF
      IF (kz(28).GT.1.AND.dsep.GT.0.0) THEN
         WRITE (6,5)  ttot, name(i1), ksx, ecc, semi*su, dsep/semi
    5    FORMAT (' BRAKE    T NAM KX* E A DA/A ',
     &                      f8.2,i6,i4,f9.5,1p,2e10.2)
      END IF
*
*       Check collision condition for stars or degenerate objects.
      IF (rp.LT.rcoal) THEN
         CALL ksperi(ipair)
         iqcoll = -2
         kspair = ipair
         CALL cmbody(r(ipair),2)
         goto 50
      ENDIF
      IF(dsep.EQ.0.d0) goto 50
*
*       Include safety test on new semi-major axis.
      rchck = min(radius(i1),radius(i2)) 
      IF(semi1.LT.rchck) semi1 = rchck 
*
*       Transform to pericentre (R > A & unperturbed).
*     IF(R(IPAIR).GT.SEMI.AND.LIST(1,I1).EQ.0)THEN
*        CALL KSAPO(IPAIR)
*     ENDIF
*
*       Form square of regularized velocity.
      v20 = 0.0
      DO 10 k = 1,4
         v20 = v20 + udot(k,ipair)**2
   10 CONTINUE
*
*       Update binding energy and collision energy.
      hi = h(ipair)
      h(ipair) = -0.5*body(i)/semi1
      zmu = body(i1)*body(i2)/body(i) 
      ecoll = ecoll + zmu*(hi - h(ipair))
      egrav = egrav + zmu*(hi - h(ipair))
*
*       Distinguish between update of eccentric binaries and standard case.
      IF(abs(ecc-ecc1).GT.tiny)THEN
*       Change KS variables at original ecc and modify ECC at H = const.
         CALL expand(ipair,semi)
         CALL ksperi(ipair)
         CALL deform(ipair,ecc,ecc1)
      ELSE
*       Specify KS coordinate & velocity scaling factors at arbitrary point.
         c2 = sqrt(semi1/semi)
*        V2 = 0.5*(BODY(I) + H(IPAIR)*SEMI1*(1.0 - ECC))
         v2 = 0.5*(body(i) + h(ipair)*r(ipair)*(semi1/semi))
         c1 = sqrt(v2/v20)
*
*       Re-scale KS variables to new energy with constant eccentricity.
         r(ipair) = 0.0d0
*        TDOT2(IPAIR) = 0.0D0
         DO 20 k = 1,4
            u(k,ipair) = c2*u(k,ipair)
            udot(k,ipair) = c1*udot(k,ipair)
            u0(k,ipair) = u(k,ipair)
            r(ipair) = r(ipair) + u(k,ipair)**2
*           TDOT2(IPAIR) = TDOT2(IPAIR) + 2.0*U(K,IPAIR)*UDOT(K,IPAIR)
   20    CONTINUE
      END IF
*
*       Transform back to apocentre for standard unperturbed motion.
*     IF (LIST(1,I1).EQ.0) THEN
*        CALL KSAPO(IPAIR)
*     END IF
*
*       Include new initialization for perturbed orbit.
      IF (list(1,i1).GT.0) THEN
         imod = kslow(ipair)
         CALL kspoly(ipair,imod)
      ENDIF
*
*       Include some diagnostic output.
      IF(kstar(i).EQ.13)THEN
         ndiag = ndiag + 1
         IF(ndiag.LT.100.OR.mod(ndiag,100).EQ.100)THEN
            rcoll = radius(i1) + radius(i2)
            WRITE (6,25)  ttot, ipair, m1, m2, r1, r2, r(ipair),
     &                    rcoll
   25       FORMAT (' BRAKE    T KS M12 R12 R RCOLL ',
     &                         f10.4,i4,2f6.2,2f7.3,1p,2e10.2)
         ENDIF
      ENDIF
*
   50 RETURN
*
      END