Nbody6
induce.f
Go to the documentation of this file.
1  SUBROUTINE induce(IPAIR,EMAX,EMIN,ICIRC,TC2,ANGLE,TG,EDAV)
2 *
3 *
4 * Induced eccentricity of hierarchical binary.
5 * --------------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 a1(3),a2(3),xrel(3),vrel(3),ei(3),hi(3),ho(3),bhat(3)
9 *
10 *
11 * Define c.m. & KS indices and evaluate semi-major axis & eccentricity.
12  i = n + ipair
13  i1 = 2*ipair - 1
14  i2 = i1 + 1
15  semi = -0.5d0*body(i)/h(ipair)
16  ecc2 = (1.d0 - r(ipair)/semi)**2 + tdot2(ipair)**2
17  & /(body(i)*semi)
18  ecc = sqrt(ecc2)
19 *
20 * Ensure that KS components are resolved (perturbed not always done).
21  CALL resolv(ipair,1)
22 *
23 * Obtain inner & outer angular momentum by cyclic notation.
24  rij2 = 0.d0
25  vij2 = 0.d0
26  rdot = 0.d0
27  a12 = 0.d0
28  a22 = 0.d0
29  a1a2 = 0.d0
30  ri2 = 0.d0
31  vi2 = 0.d0
32  rvi = 0.d0
33  DO 5 k = 1,3
34  rij2 = rij2 + (x(k,i) - x(k,jcomp))**2
35  vij2 = vij2 + (xdot(k,i) - xdot(k,jcomp))**2
36  rdot = rdot + (x(k,i) - x(k,jcomp))*(xdot(k,i) -xdot(k,jcomp))
37  k1 = k + 1
38  IF (k1.GT.3) k1 = 1
39  k2 = k1 + 1
40  IF (k2.GT.3) k2 = 1
41  a1(k) = (x(k1,i1) - x(k1,i2))*(xdot(k2,i1) - xdot(k2,i2))
42  & - (x(k2,i1) - x(k2,i2))*(xdot(k1,i1) - xdot(k1,i2))
43  a2(k) = (x(k1,jcomp) - x(k1,i))*(xdot(k2,jcomp) - xdot(k2,i))
44  & - (x(k2,jcomp) - x(k2,i))*(xdot(k1,jcomp) - xdot(k1,i))
45  a12 = a12 + a1(k)**2
46  a22 = a22 + a2(k)**2
47  a1a2 = a1a2 + a1(k)*a2(k)
48 * Form relative vectors and scalars for inner binary.
49  xrel(k) = x(k,i1) - x(k,i2)
50  vrel(k) = xdot(k,i1) - xdot(k,i2)
51  ri2 = ri2 + xrel(k)**2
52  vi2 = vi2 + vrel(k)**2
53  rvi = rvi + xrel(k)*vrel(k)
54  5 CONTINUE
55 *
56 * Evaluate orbital parameters for outer orbit.
57  rij = sqrt(rij2)
58  zmb = body(i) + body(jcomp)
59  semi1 = 2.d0/rij - vij2/zmb
60  semi1 = 1.d0/semi1
61  ecc1 = sqrt((1.d0 - rij/semi1)**2 + rdot**2/(semi1*zmb))
62 * PMIN1 = SEMI1*(1.d0 - ECC1)
63 * Determine inclination (8 bins of 22.5 degrees).
64  fac = a1a2/sqrt(a12*a22)
65  fac = acos(fac)
66  angle = fac
67 *** ANGLE = FAC*360.d0/TWOPI
68  in = 1 + fac*360.0/(twopi*22.5)
69 *
70 * Exit on hyperbolic orbit or outer eccentricity near 1.
71  IF (semi1.LT.0.0.OR.ecc1.GT.0.9999) THEN
72  emax = 0.d0
73  tc2 = 999.d0
74  tg = 1.0d+04
75  edav = 1.0d+04
76  go to 30
77  END IF
78 *
79 * Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
80  ei2 = 0.d0
81  DO 10 k = 1,3
82  ei(k) = (vi2*xrel(k) - rvi*vrel(k))/body(i) -
83  & xrel(k)/sqrt(ri2)
84  ei2 = ei2 + ei(k)**2
85  10 CONTINUE
86 *
87 * Define unit vectors for inner eccentricity and angular momenta.
88  cosj = 0.d0
89  sjsg = 0.d0
90  DO 15 k = 1,3
91  ei(k) = ei(k)/sqrt(ei2)
92  hi(k) = a1(k)/sqrt(a12)
93  ho(k) = a2(k)/sqrt(a22)
94  cosj = cosj + hi(k)*ho(k)
95  sjsg = sjsg + ei(k)*ho(k)
96  15 CONTINUE
97 *
98 * Form unit vector BHAT and scalars AH & BH (Douglas Heggie, 10/9/96).
99  ah = 0.d0
100  bh = 0.d0
101  DO 20 k = 1,3
102  k1 = k + 1
103  IF (k1.GT.3) k1 = 1
104  k2 = k1 + 1
105  IF (k2.GT.3) k2 = 1
106  bhat(k) = hi(k1)*ei(k2) - hi(k2)*ei(k1)
107  ah = ah + ei(k)*ho(k)
108  bh = bh + bhat(k)*ho(k)
109  20 CONTINUE
110 *
111 * Evaluate the expressions A & Z.
112  a = cosj*sqrt(1.d0 - ei2)
113  z = (1.d0 - ei2)*(2.d0 - cosj**2) + 5.d0*ei2*sjsg**2
114 *
115 * Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
116  z2 = z**2 + 25.d0 + 16.d0*a**4 - 10.d0*z -
117  & 20.d0*a**2 - 8.d0*a**2*z
118  emax = one6*(z + 1.d0 - 4.d0*a**2 + sqrt(z2))
119  emax = max(emax,0.0001d0)
120  emax = sqrt(emax)
121  zi = fac
122  em = sqrt(sin(zi)**2 + ei2*cos(zi)**2)
123 *
124 * Form minimum eccentricity (Douglas Heggie, Sept. 1996).
125  az = a**2 + z - 2.d0
126  IF (az.GE.0.0) THEN
127  az1 = 1.d0 + z - 4.d0*a**2
128  emin2 = one6*(az1 - sqrt(az1**2 - 12.d0*az))
129  ELSE
130  emin2 = 1.d0 - 0.5d0*(a**2 + z)
131  END IF
132  emin2 = max(emin2,0.0001d0)
133  emin = sqrt(emin2)
134 *
135 * Set impact parameters and estimate eccentricity growth time-scale.
136  pmin = semi*(1.d0 - ecc)
137  pmin2 = semi*(1.d0 - emax)
138  tk = twopi*semi*sqrt(semi/body(i))
139  tk1 = twopi*semi1*sqrt(semi1/zmb)
140  tg = tk1**2*zmb*(1.d0 - ecc1**2)**1.5/(body(jcomp)*tk)
141 *
142  const = pfac(a,z)
143  const = const*4.d0/(1.5d0*twopi*sqrt(6.d0))
144 * Convert growth time to units of 10**6 yrs.
145  tg = const*tg*tstar
146 *
147 * Form doubly averaged eccentricity derivative (Douglas Heggie 9/96).
148  yfac = 15.d0*body(jcomp)/(4.d0*zmb)*twopi*tk/tk1**2
149  yfac = yfac*ecc*sqrt(1.d0 - ecc**2)/
150  & (1.d0 - ecc1**2)**(1.5)
151  edav = yfac*ah*bh
152 *
153 * Determine circularization time for current & smallest pericentre.
154  tc2 = 2000.d0
155  tc = 2000.d0
157  & kz(27).EQ.2.AND.ecc.GT.0.002) THEN
158  CALL tcirc(pmin,ecc,i1,i2,icirc,tc)
159  IF (icirc.GE.0) icirc = -1
160 * Note: only use PMIN call because ICIRC may become >0 after first!
161  CALL tcirc(pmin2,emax,i1,i2,icirc,tc2)
162  ELSE IF (ecc.LT.0.002) THEN
163  tc = 0.0
164  tc2 = 0.0
165  END IF
166 *
167 * Print diagnostics for high inclinations and TC2 < 10**7 yrs.
168  IF (in.GE.4.AND.in.LE.5.AND.tc2.LT.10.0.AND.kstar(i).NE.-2.AND.
169  & ei2.GT.4.0d-06) THEN
170  WRITE (44,25) sqrt(ei2), em, emax, kstar(i1), kstar(i2),
171  & kstar(i), semi, pmin2, in, tg, tc, tc2, tphys
172  25 FORMAT (' INDUCE: E EMX EMAX K* SEMI PM2 IN TG TC TC2 TP ',
173  & 3f8.4,3i3,1p,2e9.1,0p,i3,f7.3,3f7.1)
174  CALL flush(44)
175  END IF
176 *
177 * Use first value for now.
178  tc2 = tc
179 *
180  30 RETURN
181 *
182  END