eccmod.f

Go to the documentation of this file.

      SUBROUTINE eccmod(I,ITERM)
*
*
*       Eccentricity modulation 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)
      REAL*8  bodyi(2),w(2)
      DATA  itime,idelay /1,0/
      SAVE  itime,idelay
*
*
*       Determine merger & ghost index.
      iterm = 0
      CALL findj(i,ighost,im)
*
*       Quit for tidal dissipation or circular orbit (TMDIS set by IMPACT).
      IF (kstarm(im).EQ.-2.OR.kstarm(im).GE.10) THEN
          go to 50
      END IF
*
*       Initialize delay indicator and pair index.
      iq = 0
      ipair = i - n
*
*       Skip on hyperbolic outer orbit, double merger or circularized orbit.
      IF (h(ipair).GT.0.0.OR.name(i).LT.-2*nzero.OR.
     &    kstarm(im).EQ.10) THEN
          tmdis(im) = time + 10.0/tstar
          go to 50
      END IF
*
*       Resolve coordinates & velocities (first time only).
      CALL resolv(ipair,1)
*
*       Form inner distance, square KS velocity and radial velocity term.
      ri = 0.0
      v20 = 0.0
      td2 = 0.0
      DO 5 k = 1,4
          ri = ri + um(k,im)**2
          v20 = v20 + umdot(k,im)**2
          td2 = td2 + 2.0*um(k,im)*umdot(k,im)
    5 CONTINUE
*
*       Evaluate inner semi-major axis and eccentricity.
      zmb = cm(1,im) + cm(2,im)
      semi = -0.5*zmb/hm(im)
      IF (semi.LE.0.0) THEN
          tmdis(im) = time + 1.0
          WRITE(6,*)' ECCMOD ERROR ',semi,ecc,h(ipair),hm(im),zmb
          go to 50
      END IF
      ecc2 = (1.0 - ri/semi)**2 + td2**2/(zmb*semi)
      ecc = sqrt(ecc2)
*
*       Obtain growth time and modify KS elements from de/dt & dh/dt.
      i1 = 2*ipair - 1
      CALL higrow(i1,ighost,im,ecc,semi,emax,emin,tg,edav,zi,iq)
*
*       Check termination for new CHAOS & SPIRAL or collision.
      IF (iq.LT.0) THEN
          iterm = 1
          itime = 1
          go to 50
      END IF
*
*       Delay for long growth times or aged active SPIRAL.
      IF ((kstarm(im).GE.0.AND.tg.GT.20.0).OR.iq.GT.0.OR.
     &    (kstarm(im).EQ.-2.AND.max(ecc,emax).LT.0.9).OR.
     &    (kstarm(im).LT.0.AND.ecc.LT.0.1)) THEN
          rm = max(radius(i1),radius(ighost),1.0d-20)
          idelay = idelay + 1
          IF (emax.GT.0.99.AND.mod(idelay,10).EQ.0) THEN
              alph = 360.0*zi/twopi
              WRITE (6,10)  name(i1), iq, kstarm(im), list(1,i1), ecc,
     &                      emax, tg, semi*(1.0-ecc)/rm, alph
          END IF
   10     FORMAT (' ECCMOD DELAY    NAM IQ K* NP E EMAX TG QP/R* IN ',
     &                              i6,3i4,2f8.4,f8.3,2f7.1)
          dt = 10.0
          IF (list(1,i1).GT.0.AND.semi*(1.0-ecc).LT.5.0*rm) THEN
              dt = 1.0
          END IF
          tmdis(im) = time + max(tg,dt)/tstar
          itime = 1
          go to 50
      END IF
*
*       Estimate current t_{circ} and de/dt for relevant condition.
      pmin = semi*(1.0 - ecc)
      rm = max(radius(i1),radius(ighost),1.0d-20)
      IF (pmin.LT.50.0*rm) THEN
          bodyi(1) = cm(1,im)
          bodyi(2) = cm(2,im)
          CALL hicirc(pmin,ecc,i1,ighost,bodyi,tg,tc,ec,edt,w)
      ELSE
          tc = 1.0d+10
          edt = 1.0d-10
      END IF
*
*       Include diagnostics every 5000th time.
      IF (itime.EQ.1.OR.mod(itime,5000).EQ.0) THEN
          a1 = -0.5*body(i)/h(ipair)
          e2 = (1.0 - r(ipair)/a1)**2 + tdot2(ipair)**2/(a1*body(i))
          ecc1 = sqrt(e2)
          zid = 360.0*zi/twopi
          np = list(1,i1)
          yc = r0(ipair)/semi
          WRITE (6,15)  name(i1), np, ttot, ecc, emax, ecc1, pmin/rm,
     &                  yc, tg, tc, edav, zid
   15     FORMAT (' ECCMOD    NM NP T E EX E1 QP/R* PC/A TG TC EDA IN ',
     &                        i6,i4,f11.4,3f8.4,2f7.1,1p,3e9.1,0p,f9.3)
          CALL flush(3)
      END IF
*
      nemod = nemod + 1
      itime = itime + 1
*
*       Include termination for expected short circularization time.
      IF (kstarm(im).GE.0.AND.pmin.LT.3.0*rm) THEN
          zid = 360.0*zi/twopi
          WRITE (6,40)  itime, emax, ecc, semi*(1.0 - ecc)/rm, zid, tc
   40     FORMAT (' ECCMOD TERM    IT EX E QP/R IN TC ',
     &                             i5,2f9.5,f6.2,f8.2,f8.1)
          CALL flush(3)
          iterm = 1
          itime = 1
      END IF
*
   50 RETURN
*
      END