Nbody6
 All Files Functions Variables
eccmod.f
Go to the documentation of this file.
1  SUBROUTINE eccmod(I,ITERM)
2 *
3 *
4 * Eccentricity modulation of hierarchical binary.
5 * -----------------------------------------------
6 *
7  include 'common6.h'
8  common/binary/ cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
9  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
10  & namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
11  REAL*8 bodyi(2),w(2)
12  DATA itime,idelay /1,0/
13  SAVE itime,idelay
14 *
15 *
16 * Determine merger & ghost index.
17  iterm = 0
18  CALL findj(i,ighost,im)
19 *
20 * Quit for tidal dissipation or circular orbit (TMDIS set by IMPACT).
21  IF (kstarm(im).EQ.-2.OR.kstarm(im).GE.10) THEN
22  go to 50
23  END IF
24 *
25 * Initialize delay indicator and pair index.
26  iq = 0
27  ipair = i - n
28 *
29 * Skip on hyperbolic outer orbit, double merger or circularized orbit.
30  IF (h(ipair).GT.0.0.OR.name(i).LT.-2*nzero.OR.
31  & kstarm(im).EQ.10) THEN
32  tmdis(im) = time + 10.0/tstar
33  go to 50
34  END IF
35 *
36 * Resolve coordinates & velocities (first time only).
37  CALL resolv(ipair,1)
38 *
39 * Form inner distance, square KS velocity and radial velocity term.
40  ri = 0.0
41  v20 = 0.0
42  td2 = 0.0
43  DO 5 k = 1,4
44  ri = ri + um(k,im)**2
45  v20 = v20 + umdot(k,im)**2
46  td2 = td2 + 2.0*um(k,im)*umdot(k,im)
47  5 CONTINUE
48 *
49 * Evaluate inner semi-major axis and eccentricity.
50  zmb = cm(1,im) + cm(2,im)
51  semi = -0.5*zmb/hm(im)
52  IF (semi.LE.0.0) THEN
53  tmdis(im) = time + 1.0
54  WRITE(6,*)' ECCMOD ERROR ',semi,ecc,h(ipair),hm(im),zmb
55  go to 50
56  END IF
57  ecc2 = (1.0 - ri/semi)**2 + td2**2/(zmb*semi)
58  ecc = sqrt(ecc2)
59 *
60 * Obtain growth time and modify KS elements from de/dt & dh/dt.
61  i1 = 2*ipair - 1
62  CALL higrow(i1,ighost,im,ecc,semi,emax,emin,tg,edav,zi,iq)
63 *
64 * Check termination for new CHAOS & SPIRAL or collision.
65  IF (iq.LT.0) THEN
66  iterm = 1
67  itime = 1
68  go to 50
69  END IF
70 *
71 * Delay for long growth times or aged active SPIRAL.
72  IF ((kstarm(im).GE.0.AND.tg.GT.20.0).OR.iq.GT.0.OR.
73  & (kstarm(im).EQ.-2.AND.max(ecc,emax).LT.0.9).OR.
74  & (kstarm(im).LT.0.AND.ecc.LT.0.1)) THEN
75  rm = max(radius(i1),radius(ighost),1.0d-20)
76  idelay = idelay + 1
77  IF (emax.GT.0.99.AND.mod(idelay,10).EQ.0) THEN
78  alph = 360.0*zi/twopi
79  WRITE (6,10) name(i1), iq, kstarm(im), list(1,i1), ecc,
80  & emax, tg, semi*(1.0-ecc)/rm, alph
81  END IF
82  10 FORMAT (' ECCMOD DELAY NAM IQ K* NP E EMAX TG QP/R* IN ',
83  & i6,3i4,2f8.4,f8.3,2f7.1)
84  dt = 10.0
85  IF (list(1,i1).GT.0.AND.semi*(1.0-ecc).LT.5.0*rm) THEN
86  dt = 1.0
87  END IF
88  tmdis(im) = time + max(tg,dt)/tstar
89  itime = 1
90  go to 50
91  END IF
92 *
93 * Estimate current t_{circ} and de/dt for relevant condition.
94  pmin = semi*(1.0 - ecc)
95  rm = max(radius(i1),radius(ighost),1.0d-20)
96  IF (pmin.LT.50.0*rm) THEN
97  bodyi(1) = cm(1,im)
98  bodyi(2) = cm(2,im)
99  CALL hicirc(pmin,ecc,i1,ighost,bodyi,tg,tc,ec,edt,w)
100  ELSE
101  tc = 1.0d+10
102  edt = 1.0d-10
103  END IF
104 *
105 * Include diagnostics every 5000th time.
106  IF (itime.EQ.1.OR.mod(itime,5000).EQ.0) THEN
107  a1 = -0.5*body(i)/h(ipair)
108  e2 = (1.0 - r(ipair)/a1)**2 + tdot2(ipair)**2/(a1*body(i))
109  ecc1 = sqrt(e2)
110  zid = 360.0*zi/twopi
111  np = list(1,i1)
112  yc = r0(ipair)/semi
113  WRITE (6,15) name(i1), np, ttot, ecc, emax, ecc1, pmin/rm,
114  & yc, tg, tc, edav, zid
115  15 FORMAT (' ECCMOD NM NP T E EX E1 QP/R* PC/A TG TC EDA IN ',
116  & i6,i4,f11.4,3f8.4,2f7.1,1p,3e9.1,0p,f9.3)
117  CALL flush(3)
118  END IF
119 *
120  nemod = nemod + 1
121  itime = itime + 1
122 *
123 * Include termination for expected short circularization time.
124  IF (kstarm(im).GE.0.AND.pmin.LT.3.0*rm) THEN
125  zid = 360.0*zi/twopi
126  WRITE (6,40) itime, emax, ecc, semi*(1.0 - ecc)/rm, zid, tc
127  40 FORMAT (' ECCMOD TERM IT EX E QP/R IN TC ',
128  & i5,2f9.5,f6.2,f8.2,f8.1)
129  CALL flush(3)
130  iterm = 1
131  itime = 1
132  END IF
133 *
134  50 RETURN
135 *
136  END