higrow.f

Go to the documentation of this file.

      SUBROUTINE higrow(I,IG,IM,ECC,SEMI,EMAX,EMIN,TG,EDAV,ZI,IQ)
*
*
*       Induced change of hierarchical binary.
*       --------------------------------------
*
      include 'common6.h'
      common/binary/  cm(4,mmax),yrel(3,mmax),zrel(3,mmax),
     &                hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
     &                namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
      common/modes/  eb0(ntmax),zj0(ntmax),ecrit(ntmax),ar(ntmax),
     &               br(ntmax),eosc(4,ntmax),edec(ntmax),tosc(ntmax),
     &               rp(ntmax),es(ntmax),zm(2,ntmax),iosc(ntmax),
     &               namec(ntmax)
      common/slow0/  range,islow(10)
      common/tidal/  cq(2),ct(2),cgr,dedt
      common/rksave/  coeff,hohat(3),e0,a0,hh,xmb
      REAL*8  a1(3),a2(3),xrel(3),vrel(3),ei(3),hi(3),ho(3),bhat(3),
     &        ui(4),v(4),evec(3),xr0(3),vr0(3),ei0(3),ww(6),w(4)
      REAL*8  bodyi(2),qg(2),wg(2)
      LOGICAL icoll
      DATA  ww  /2.119,3.113,8.175,3.742,4.953,9.413/
      SAVE  icall,itime,itry,namei,dtprev,tcheck,zfac,pmin1
      DATA  icall,itry,namei,tcheck  /0,0,0,0.0d0/
*
*
*       Copy KS variables to local scalars.
      DO 1 k = 1,4
          ui(k) = um(k,im)
          v(k) = umdot(k,im)
    1 CONTINUE
*
*       Transform to physical variables and multiply by 4 (momentum formula).
      CALL ksphys(ui,v,xrel,vrel)
      DO 2 k = 1,3
          vrel(k) = 4.0*vrel(k)
    2 CONTINUE
*
*       Specify index for outer body (second KS component) & pair index.
      jcomp = i + 1
      ipair = kvec(i)
      dtsum = 0.0
      itime = 0
      icall = 0
      ecc0 = ecc
      semi0 = semi
      a0 = semi
      hm0 = hm(im)
      zm3 = body(jcomp)
      bodyi(1) = cm(1,im)
      bodyi(2) = cm(2,im)
*
*       Obtain inner & outer angular momentum by cyclic notation.
    4 a12 = 0.0
      a22 = 0.0
      a1a2 = 0.0
      ri2 = 0.0
      vi2 = 0.0
      rvi = 0.0
      DO 5 k = 1,3
          k1 = k + 1
          IF (k1.GT.3) k1 = 1
          k2 = k1 + 1
          IF (k2.GT.3) k2 = 1
          a1(k) = xrel(k1)*vrel(k2) - xrel(k2)*vrel(k1)
          a2(k) = (x(k1,jcomp) - x(k1,i))*(xdot(k2,jcomp) - xdot(k2,i))
     &          - (x(k2,jcomp) - x(k2,i))*(xdot(k1,jcomp) - xdot(k1,i))
          a12 = a12 + a1(k)**2
          a22 = a22 + a2(k)**2
          a1a2 = a1a2 + a1(k)*a2(k)
          ri2 = ri2 + xrel(k)**2
          vi2 = vi2 + vrel(k)**2
          rvi = rvi + xrel(k)*vrel(k)
    5 CONTINUE
*
*       Evaluate orbital parameters for outer orbit from KS elements.
      zmb = body(i) + body(jcomp)
      semi1 = -0.5*zmb/h(ipair)
      ecc2 = (1.0 - r(ipair)/semi1)**2 + tdot2(ipair)**2/(semi1*zmb)
      ecc1 = sqrt(ecc2)
*       Determine inclination in radians.
      fac = a1a2/sqrt(a12*a22)
      IF(fac.LT.-1.d0.OR.fac.GT.1.d0)THEN
         IF(fac.LT.-1.05d0.OR.fac.GT.1.05d0)THEN
            WRITE (6,9)  ipair, i, jcomp, ecc1, semi1*su, fac 
   9        FORMAT (' WARNING! HIGROW FAC: IP I J E A FAC ',
     &                                   3i6,f7.3,f8.1,f9.5)
         ENDIF
         fac = max(fac,-1.d0) 
         fac = min(fac,1.d0) 
      ENDIF
      zi = acos(fac)
      cc = (1.0 - ecc**2)*cos(zi)**2
*
*       Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
      ei2 = 0.0
      DO 10 k = 1,3
          ei(k) = (vi2*xrel(k) - rvi*vrel(k))/body(i) -
     &                                                 xrel(k)/sqrt(ri2)
          ei2 = ei2 + ei(k)**2
          evec(k) = ei(k)
   10 CONTINUE
      ei2 = min(ei2,0.9999d0)
*
*       Define unit vectors for inner eccentricity and angular momenta.
      cosj = 0.0
      sjsg = 0.0
      DO 15 k = 1,3
          ei(k) = ei(k)/sqrt(ei2)
          hi(k) = a1(k)/sqrt(a12)
          ho(k) = a2(k)/sqrt(a22)
          cosj = cosj + hi(k)*ho(k)
          sjsg = sjsg + ei(k)*ho(k)
   15 CONTINUE
*
*       Form unit vector BHAT and scalars AH & BH (Douglas Heggie, 10/9/96).
      ah = 0.0
      bh = 0.0
      DO 16 k = 1,3
          k1 = k + 1
          IF (k1.GT.3) k1 = 1
          k2 = k1 + 1
          IF (k2.GT.3) k2 = 1
          bhat(k) = hi(k1)*ei(k2) - hi(k2)*ei(k1)
          ah = ah + ei(k)*ho(k)
          bh = bh + bhat(k)*ho(k)
   16 CONTINUE
*
*       Evaluate the expressions A & Z.
      a = cosj*sqrt(1.0 - ei2)
      z = (1.0 - ei2)*(2.0 - cosj**2) + 5.0*ei2*sjsg**2
*
*       Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
      z2 = z**2 + 25.0 + 16.0*a**4 - 10.0*z - 20.0*a**2 - 8.0*a**2*z
      emax = one6*(z + 1.0 - 4.0*a**2 + sqrt(z2))
      emax = max(emax,0.0001d0)
      emax = sqrt(emax)
*
*       Form minimum eccentricity (Douglas Heggie, Sept. 1996).
      az = a**2 + z - 2.0
      IF (az.GE.0.0) THEN
          az1 = 1.0 + z - 4.0*a**2
          emin2 = one6*(az1 - sqrt(az1**2 - 12.0*az))
      ELSE
          emin2 = 1.0 - 0.5*(a**2 + z)
      END IF
      emin2 = max(emin2,0.0001d0)
      emin = sqrt(emin2)
*
*       Estimate eccentricity growth time-scale (N-body units).
      tk = twopi*semi*sqrt(semi/body(i))
      tk1 = twopi*abs(semi1)*sqrt(abs(semi1)/zmb)
      tg = tk1**2*zmb*(1.0 - ecc1**2)**1.5/(body(jcomp)*tk)
*
*       Evaluate numerical precession factor (involves elliptic integral).
      const = pfac(a,z)
      const = const*4.0/(1.5*twopi*sqrt(6.0))
*
*       Convert growth time to units of 10**6 yrs.
      tg = const*tg*tstar
*
*       Form doubly averaged eccentricity derivative (Douglas Heggie 9/96).
      yfac = 15.0*body(jcomp)/(4.0*zmb)*twopi*tk/tk1**2
      yfac = yfac*ecc*sqrt(1.0 - ecc**2)/(1.0 - ecc1**2)**(1.5)
      edav = yfac*ah*bh
*
*       Skip small EMAX or large pericentre distance (TMDIS increased).
      qperi = semi*(1.0 - ecc)
      pmin = semi*(1.0 - emax)
      rm = max(radius(i),radius(ig),1.0d-20)
      IF (kstarm(im).GE.0.AND.emax.LT.0.9) THEN
          iq = 1
          go to 40
      END IF
*
*       Delay Kozai cycle for EMAX > 0.9 and TC > 2000.
      IF (emax.GE.0.9) THEN
          CALL hicirc(pmin,emax,i,ig,bodyi,tg,tc,ec,edt,w)
          IF (tc.GT.2000.OR.(ecc.LT.0.1.AND.emax.LT.0.95)) THEN
              iq = 2
              go to 40
          END IF
      END IF
*
*       Check termination time for marginal stability criterion.
      IF (name(i).EQ.namei.AND.(time+toff).GT.tcheck) THEN
          WRITE (6,17)  name(i), ecc, emax, alph, zfac, ecc1, pmin1,
     &                  pcrit
   17     FORMAT (' ECCMOD UNSTAB    NM E EX IN YF E1 PM PC ',
     &                               i6,2f8.4,f7.1,f6.2,f7.3,1p,2e10.2)
          iq = -4
          go to 40
      END IF
*
*       Ensure evaluation of EDT for large eccentricity and EDAV > 0.
      IF (ecc.GT.0.9.AND.edav.GT.0.0) THEN
          CALL hicirc(qperi,ecc,i,ig,bodyi,tg,tc,ec,edt,w)
          IF (tc.LT.2000.0) THEN
              icoll = .true.
          ELSE
              icoll = .false.
              edt = edav
          END IF
      ELSE
          icoll = .false.
          edt = edav
      END IF
*
*       Activate indicator on small QPERI to prevent collision before EMAX.
      IF (qperi.LT.2.0*rm) icoll = .true.
*
*       Check circularization time near turning point of normal binary.
      IF ((kstarm(im).GE.0.AND.edav.GT.0.0.AND.
     &   (dedt.LT.0.0.OR.ecc.GT.abs(emax-0.00005))).OR.icoll) THEN
          IF (.NOT.icoll) THEN
              CALL hicirc(qperi,ecc,i,ig,bodyi,tg,tc,ec,edt,w)
          END IF
          IF (tc.LT.7.0d+04.AND.qperi.LT.6.0*rm.AND.itry.LT.4.OR.
     &        tc.LT.1.0d+04.AND.qperi.LT.5.0*rm.AND.itry.LT.8) THEN
              WRITE (6,18)  ecc, emax, emin, semi*su, tg, tc, edav,
     &                      edt, qperi/rm
   18         FORMAT (' HICIRC TRY:    E EX EM A TG TC EDA EDT QP/R ',
     &                                 2f9.5,f7.3,f7.1,1p,5e9.1)
              itry = itry + 1
          END IF
          IF (itry.GT.8) itry = 0
*
*       Check termination for TC < 2000 to be consistent with routine TCIRC.
          IF (tc.GT.2000.0) go to 25
*
*       Evaluate chaos boundary parameters.
          eb = -0.5*cm(1,im)*cm(2,im)/a0
          cj = cm(1,im)*cm(2,im)/body(i)*sqrt(body(i))
          zj = cj*sqrt(qperi*(1.0 + ecc))
          ic = nchaos + 1
          idis = 0
          itry = 0
*
*       Define oscillation period (dimensionless time).
          DO 20 k = 1,2
              IF (k.EQ.1) THEN
                  ik = i
              ELSE
                  ik = ig
              END IF
*       Specify polytropic index for each star (n = 3, 2 or 3/2).
              IF (kstar(ik).EQ.3.OR.kstar(ik).EQ.5) THEN
                  bodi = cm(k,im)
                  CALL giant3(ik,bodi,wg,qg,xn,ql)
                  w(k) = wg(1)
              ELSE
                  ip = 3
                  IF (kstar(ik).GE.3) ip = 2
                  IF (kstar(ik).EQ.4.OR.kstar(ik).EQ.6) ip = 3
                  IF (kstar(ik).EQ.0) ip = 1
                  w(k) = ww(ip)
              END IF
   20     CONTINUE
*
          CALL chaos0(qperi,ecc,eb,zj,cm(1,im),cm(2,im),radius(i),
     &                      radius(ig),w,ecrit(ic),ar(ic),br(ic),idis)
*
*       Begin chaos/spiral stage if chaos boundary has been crossed.
          IF (idis.EQ.0) THEN
              WRITE (6,22)  ttot, name(i), ic, ecc, emax, semi, qperi,
     &                      edav, qperi/rm
   22         FORMAT (' ECCMOD CHAOS    T NAM IC E EX A QP EDAV QP/R ',
     &                                  f9.2,i6,i4,2f9.5,1p,4e10.2)
              CALL flush(6)
              iq = -1
*       Activate chaos indicator for calling KSTIDE from RESET.
              kstarm(im) = -1
              go to 40
      ELSE IF (idis.EQ.-1) THEN
              WRITE (6,23)  ttot, name(i), ic, ecc, emax, semi, qperi,
     &                      edav, qperi/rm
   23         FORMAT (' ECCMOD SPIRAL    T NAM IC E EX A QP EDAV QP/R ',
     &                                   f9.2,i6,i4,2f9.5,1p,4e10.2)
              CALL flush(6)
              iq = -2
*       Note that KSTARM(IM) = -2 would be treated as existing spiral.
              kstarm(im) = -1
              go to 40
          ELSE IF (idis.EQ.1) THEN
*       Check collision condition to be sure.
              IF (qperi.LT.radius(i) + radius(ig)) THEN
              WRITE (6,24)  ttot, name(i), ic, ecc, semi, qperi, edav,
     &                      qperi/rm
   24         FORMAT (' ECCMOD DISRUPT    T NAM IC E A QP EDAV QP/R  ',
     &                                    f9.2,i6,i4,f8.4,1p,4e10.2)
                  iq = -3
                  go to 40
              END IF
          END IF
      END IF
*
*       Form quantities for the outer orbit.
   25 rr = 0.0
      rv0 = 0.0
      v20 = 0.0
      DO 30 k = 1,3
          xr0(k) = x(k,jcomp) - x(k,i)
          vr0(k) = xdot(k,jcomp) - xdot(k,i)
          rv0 = rv0 + xr0(k)*vr0(k)
          v20 = v20 + vr0(k)**2
          rr = rr + xr0(k)**2
   30 CONTINUE
*
*       Construct the outer Runge-Lenz vector.
      DO 35 k = 1,3
          ei0(k) = (v20*xr0(k) - rv0*vr0(k))/zmb - xr0(k)/sqrt(rr)
   35 CONTINUE
*
*       Specify total time interval (unperturbed or perturbed case).
      IF (list(1,i).EQ.0) THEN
          dt0 = time - t0(i)
*         DT = MIN(DT0,0.1*(1.0 - ECC)/(ABS(EDAV) + 1.0D-15))
      ELSE
*       Ensure that interval extends to next apocentre (including slow-down).
          imod = kslow(ipair)
          dt0 = 0.98*float(islow(imod))*tk1
*         DT = MIN(DT0,0.1*(1.0 - ECC)*TG/TSTAR)
      END IF
*
*       Set partial interval and check remaining time.
      dt1 = min(dt0,0.1*tg*sqrt(1.0 - ecc**2)/tstar)
      IF (dtsum + dt1.GT.dt0) dt1 = dt0 - dtsum + 1.0d-12
*
*       Obtain quadrupole and tidal constants for each new binary.
      dtm = time - max(tev0(i),tev0(ig))
      IF (name(i).NE.namei.OR.dtm.LE.dt0) THEN
          CALL qtides(i,ig,im,semi,ecc)
          namei = name(i)
          dedt = 0.0
          dtprev = 1.0
          nst = nstab(semi,semi1,ecc,ecc1,zi,bodyi(1),bodyi(2),zm3)
          IF (nst.EQ.0) THEN
              pcrit = 0.98*qperi*(1.0 - pert)
              pcr = stability(bodyi(1),bodyi(2),zm3,ecc,ecc1,zi)
              pcr = pcr*semi
*       Specify reduced peri if old criterion < PCRIT/2 (avoids switching).
              IF (pcr.LT.0.5*pcrit) THEN
                  pcrit = 0.75*pcrit
              END IF
              IF (pcrit.LT.pcr.AND.pert.LT.0.01.AND.
     &            itime.LT.20) THEN
                  alph = 360.0*zi/twopi
                  fail = qperi*(1-pert) - pcr
                  WRITE (6,42)  ttot, alph, ecc, ecc1, qperi, fail, pert
   42             FORMAT (' NEWSTAB    T INC EI EO QP FAIL PERT ',
     &                                 f7.1,f7.2,2f8.4,1p,3e10.2)
              END IF
          ELSE
              pcrit = 1.01*qperi
          END IF
          pmin1 = semi1*(1.0 - ecc1)
          IF (pmin1.LT.pcrit) THEN
              nk = 1 + 10.0*ecc1/(1.0 - ecc1)
              tcheck = time + toff + nk*tk1
          ELSE
              tcheck = 1.0d+10
          END IF
      END IF
*
*       Determine time-step for integration from several criteria.
      ge=30*ecc+45*ecc**3+3.75*ecc**5
*       Apply factor of 2 correction for agreement with classical result.
      ge = 0.5*ge
      slr=semi*(1.0 - ecc**2)
      zq=ge/(slr**5*semi*sqrt(semi))
*       Correct for wrong eccentricity dependence.
      zq = zq*(1.0 - ecc**2)
      edq = zq*(cq(1) + cq(2))
      tq = ecc/edq
      zgr = cgr/(slr*semi*sqrt(semi))
*       Note cgr is angular velocity (hence multiply by 2*pi; 12/03).
      tgr = twopi*ecc/zgr
*
*       Delay first time if apsidal motion or GR are important (ECC < 0.9).
      IF (ecc.LT.0.90.AND.itime.EQ.0) THEN
          IF (min(tq,tgr)*tstar.LT.tg) THEN
              iq = 2
              go to 40
          END IF
      END IF
*
*       Include check on tidal dissipation for large eccentricity.
      IF (kstarm(im).GE.0.AND.ecc.GT.0.9) THEN
          tt = ecc/abs(edt)
          tq = min(tt,tq)
*       Note that EDT from HICIRC agrees well with actual value from RKINT.
      END IF
*
*       Adopt harmonic mean of inner period and growth/quadrupole time.
      dt = 0.5*sqrt(tk*min(tg/tstar,tq))
      IF (ecc.GT.0.98) dt = 0.1*dt
      dt = min(dt,0.001/(abs(edav) + 1.0d-20))
*       Ensure conservative step first time and limit increase to factor 2.
      IF (itime.EQ.0) THEN
          dt = min(dt,1.0d-03*tg/tstar)
      ELSE
          dt = min(dt,2.0d0*dtprev)
      END IF
*
*       Include safety check to avoid problems.
      etry = ecc + edav*dt1
      IF (ecc.GT.0.96.AND.etry.GT.0.98) THEN
          etry = min(etry,emax)
          qptry = semi*(1.0 - etry)
          CALL hicirc(qptry,etry,i,ig,bodyi,tg,tc,ec,edt,w)
          IF (tc.LT.10.0) THEN
              dt1 = 0.5*dt1
              dt = 0.5*dt
          END IF
      END IF
*
*       Ensure frequent termination check inside 5*RM to avoid overshooting.
      IF (qperi.LT.5.0*rm.AND.edav.GT.0.0.AND.kstarm(im).GE.0) THEN
*       Restrict HIMOD interval so E + EDAV*DT1 gives new peri outside 4*RM.
          dtt = (qperi - 4.0*rm)/(semi*edav)
          dt1 = min(dtt,dt1)
          dt1 = max(dt,dt1)
      END IF 
*
*       Modify the inner orbit using averaging of de/dt & dh/dt.
      dtprev = dt
      CALL himod(evec,a1,zm3,zmb,ei0,a2,icall,dt1,dt,ecc,xrel,vrel)
*
      IF (e0.ge.1.0) THEN
          iq = 1
          go to 40
      END IF
      semi = a0
      ecc = e0
      itime = itime + 1
      dtsum = dtsum + dt1
      IF (dtsum.LT.dt0) go to 4
*
*       Update step counter and next time for checking.
      neint = neint + itime
      tmdis(im) = time + dtsum
*
*       Transform back to KS variables after integration.
   40 IF (itime.GT.0) THEN
          CALL physks(xrel,vrel,ui,v)
          DO 45 k = 1,4
              um(k,im) = ui(k)
              umdot(k,im) = 0.25*v(k)
              IF (k.LT.4) THEN
                  yrel(k,im) = xrel(k)
                  zrel(k,im) = vrel(k)
              END IF
   45     CONTINUE
          hm(im) = -0.5d0*body(i)/a0
          CALL hirect(im)
          zmu = cm(1,im)*cm(2,im)/body(i)
          decorr = zmu*(hm0 - hm(im))
          emerge = emerge - decorr
          ecoll = ecoll + decorr
          egrav = egrav + decorr
*       Define circularization if e0 < 0.002 (end of Kozai cycle).
          IF (e0.LT.0.002) THEN
              km = kstarm(im)
              kstarm(im) = 10
              IF (km.LT.0) ncirc = ncirc + 1
              tmdis(im) = 1.0d+10
              tev(n+ipair) = min(tev(i),tev(ig)) - 2.0*stepx
              tb = days*semi*sqrt(semi/body(i))
              alph = zi*360.0/twopi
              WRITE (6,48)  name(i), name(ig), km, e0, a0*su, tb, alph
   48         FORMAT (' END KOZAI    NM K* E A P INC ',
     &                               2i6,i4,f7.3,f6.1,2f7.1)
          END IF
*       Include rectification of active SPIRAL (suppressed).
          ic = 0
          DO 50 k = 1,nchaos
              IF (namec(k).EQ.nzero - namem(im)) ic = k
*             IF (NAMEC(K).EQ.NZERO + NAME(I)) IC = K
   50     CONTINUE
          IF (ic.EQ.-3) THEN
              rp(ic) = a0*(1.0 - ecc)
              es(ic) = ecc
              tosc(ic) = time
              dh = hm0 - hm(im)
              WRITE (6,60)  ttot, tmdis(im), ecc, rp(ic), dh
   60         FORMAT (' UPDATE:    T TMD E RP DH ',
     &                             2f10.2,f8.4,1p,2e10.2)
          END IF
      END IF
*
      RETURN
*
      END