brake4.f

Go to the documentation of this file.

      SUBROUTINE brake4(I1,I2,DT)
*
*
*       GR analytical orbit shrinkage.
*       ------------------------------
*
      include 'common6.h'
      common/postn/  cvel,taugr,rz1,gammaz,tkoz,emax,tsp,kz24,igr,ipn
      common/kspar/  istat(kmax)
      REAL*8 m1,m2,ui(4),uidot(4)
      SAVE iter
      DATA iter /0/
*
*
*       Check relativistic conditions (at least one >= NS).
      IF (max(kstar(i1),kstar(i2)).LT.13) go to 100
*
*       See whether CLIGHT has been initialized in ARchain.
      IF (iter.EQ.0) THEN
          READ (5,*)  clight
          iter = 1
      END IF
*
*       Specify the basic elements from BH or KS treatments.
      m1 = body(i1)
      m2 = body(i2)
      ipair = kvec(i1)
      i = n + ipair
      semi = -0.5*body(i)/h(ipair)
      e2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
      ecc = sqrt(e2)
*
*       Replace the BH value of CVEL by CLIGHT for standard binaries.
      IF (cvel.EQ.0.0d0) cvel = clight
*
*       Form da/dt & de/dt according to Peters 1964.
      adot =  64.0/5.0*m1*m2*(m1+m2)/(cvel**5*semi**3*(1.0 - e2)**3.5)
      adot = adot*(1.0 + 73.0/24.0*e2 + 37.0/96.0*e2**2)
      edot = 304.0/15.0*ecc*m1*m2*(m1 + m2)/(cvel**5*semi**4)
      edot = edot/(1.0 - e2)**2.5*(1.0 + 121.0/304.0*e2)
      tgr = semi/adot
*
*       Skip for long time-scale (low probability of unperturbed orbit).
      IF (tgr.GT.500.0) go to 100
*
*       Set local KS vectors.
      DO 2 k = 1,4
          ui(k) = u0(k,ipair)
          uidot(k) = udot(k,ipair)
    2 CONTINUE
*
*       Evaluate the Einstein shift per orbit and check time-step.
      dw = 3.0*twopi*(body(i1) + body(i2))/(semi*cvel**2*(1.0-e2))
      tk = twopi*semi*sqrt(semi/(body(i1) + body(i2)))
*       Adopt time-scale of 2 % relative change (subject to c.m. step).
      dt = min(0.02*semi/adot,step(i))
      theta = dw*dt/tk
*       Impose limit of time-step if THETA > TWOPI.
      IF (theta.GT.twopi) THEN
          dt = twopi*tk/dw
      END IF
      step(i1) = dt
*     THETA = DMOD(THETA,TWOPI)   ! original condition.
*
*       Rotate KS orbit by THETA/2 (period halving).
      CALL ksrot(ui,uidot,theta)
*
*       Copy back to KS common variables.
      DO 15 k = 1,4
          u(k,ipair) = ui(k)
          u0(k,ipair) = ui(k)
          udot(k,ipair) = uidot(k)
   15 CONTINUE
*
*       Define RZ and obtain the incremental change in SEMI & ECC.
      rz = 8.0*(m1 + m2)/cvel**2
      semi1 = max(semi - adot*dt,rz)
      ecc1 = max(ecc - edot*dt,0.0d0)
      IF (semi1.LT.0.5*semi) THEN
          WRITE (6,20)  name(i1), semi, semi1, adot*dt
   20     FORMAT (' PN WARNING    NM A A1 ADOT*DT ',i6,1p,3e10.2)
      END IF
*
*       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))
*
      IF (iter.LT.10000) THEN
          WRITE (93,30)  name(i1), time+toff, ecc1, semi1
   30     FORMAT (' ',i7,f10.3,f9.5,1p,e12.4)
          CALL flush(93)
      END IF
*
*       Change KS variables at original ecc and modify ECC at H = const.
      CALL expand2(ipair,semi)
      CALL ksperi(ipair)
*       Note simplified versions of standard routines (new STEP(I1) given).
      CALL deform2(ipair,ecc,ecc1)
*
      iter = iter + 1
      IF (iter.LT.1000.OR.mod(iter,1000).EQ.0) THEN
          WRITE (94,40)  time+toff, ecc, theta, dt, tgr, semi
   40     FORMAT (' GR SHRINK    T E TH DT TGR A ',
     &                           f11.4,f9.5,1p,3e9.1,e12.4)
          CALL flush(94)
      END IF
*
*       Check KS termination with added perturber to activate PN.
      IF (tgr.LT.0.5) THEN
*       Note that first order Peters formulation is not valid for strong GR.
          jp = list(2,i)
          list(1,i1) = 1
          list(2,i1) = jp
*       Set PN indicator for ARCHAIN (TGR limit means small TZ).
          ipn = 3
          WRITE (6,44)  jp, name(jp), step(i1), step(i), semi, tgr
   44     FORMAT (' ENFORCED PERTURB    JP NM S1 SI A TZ',2i6,1p,5e10.2)
          go to 100
      END IF
*
*       Activate coalescence condition using local index.
      jphase = 0
      IF (semi1.LT.100.0*rz.AND.tgr.LT.0.1) THEN
          WRITE (6,45)  kstar(i1), kstar(i2), radius(i1), radius(i2),
     &                  semi1
   45     FORMAT (' PN COAL    K* R1 R2 A  ',2i4,1p,3e10.2)
          CALL flush(6)
          iqcoll = -2
          jphase = -1
*         KSPAIR = IPAIR
*         CALL CMBODY(SEMI1,2)
*
*       Include optional kick velocity of 3*VRMS km/s for coalescence recoil.
          IF (kz(43).GT.0) THEN
*       Initialize basic variables at start of new step.
              vi20 = 0.0
              DO 48 k = 1,3
                  x0(k,i) = x(k,i)
                  x0dot(k,i) = xdot(k,i)
                  vi20 = vi20 + xdot(k,i)**2
   48         CONTINUE
*
              vf = 3.0*(vrms/vstar)/sqrt(vi20)
              DO 50 k = 1,3
                  xdot(k,i) = vf*xdot(k,i)
                  x0dot(k,i) = xdot(k,i)
   50         CONTINUE
              ecd0 = ecdot
              ecdot = ecdot + 0.5*body(i)*vi20*(1.0 - vf**2)
              vesc = 3.0*vrms
              WRITE (6,60)  vf, ecd0-ecdot, vesc
   60         FORMAT (' COALESCENCE KICK    VF ECDOT VESC ',
     &                                      f7.3,f10.6,f6.1)
*       Form neighbour list and new polynomials.
              rs0 = rs(i)
              CALL nblist(i,rs0)
              CALL fpoly1(i,i,0)
              CALL fpoly2(i,i,0)
          END IF
      END IF
*
      IF (istat(kcase).EQ.0) istat(kcase) = jphase
*
  100 RETURN
*
      END