Nbody6
 All Files Functions Variables
kozai.f
Go to the documentation of this file.
1  SUBROUTINE kozai(IPAIR,IM,ECC1,SEMI1,ITERM)
2 *
3 *
4 * Kozai relations & decision-making.
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  common/koztmp/ emax
12  REAL*8 ei(3),hi(3),ho(3),a1(3),a2(3)
13  SAVE it
14  DATA it /0/
15 *
16 *
17 * Define KS & c.m. indices for outer orbit with ECC1 & SEMI1.
18  i1 = 2*ipair - 1
19  i2 = i1 + 1
20  i = n + ipair
21  iterm = 0
22 *
23 * Evaluate inner semi-major axis and eccentricity.
24  ri = 0.0
25  td2 = 0.0
26  DO 1 k = 1,4
27  ri = ri + um(k,im)**2
28  td2 = td2 + 2.0*um(k,im)*umdot(k,im)
29  1 CONTINUE
30  zmb = cm(1,im) + cm(2,im)
31  semi = -0.5*zmb/hm(im)
32  ecc2 = (1.0 - ri/semi)**2 + td2**2/(zmb*semi)
33  ecc = sqrt(ecc2)
34 *
35 * Skip circularized binaries.
36  IF (ecc.LT.0.1) THEN
37  IF (tev(i).GE.1.0d+10) tev(i) = time + 1.0
38  go to 50
39  END IF
40 *
41 * Form useful scalars (XREL & VREL consistent with UM & UMDOT).
42  a12 = 0.0
43  a22 = 0.0
44  a1a2 = 0.0
45  ri2 = 0.0
46  vi2 = 0.0
47  rv = 0.0
48  DO 5 k = 1,3
49  k1 = k + 1
50  IF (k1.GT.3) k1 = 1
51  k2 = k1 + 1
52  IF (k2.GT.3) k2 = 1
53  a1(k) = xrel(k1,im)*vrel(k2,im) - xrel(k2,im)*vrel(k1,im)
54  a2(k) = (x(k1,i2) - x(k1,i1))*(xdot(k2,i2) - xdot(k2,i1))
55  & - (x(k2,i2) - x(k2,i1))*(xdot(k1,i2) - xdot(k1,i1))
56  a12 = a12 + a1(k)**2
57  a22 = a22 + a2(k)**2
58  a1a2 = a1a2 + a1(k)*a2(k)
59  ri2 = ri2 + xrel(k,im)**2
60  vi2 = vi2 + vrel(k,im)**2
61  rv = rv + xrel(k,im)*vrel(k,im)
62  5 CONTINUE
63 *
64 * Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
65  ei2 = 0.0
66  DO 10 k = 1,3
67  ei(k) = (vi2*xrel(k,im) - rv*vrel(k,im))/body(i1) -
68  & xrel(k,im)/sqrt(ri2)
69  ei2 = ei2 + ei(k)**2
70  10 CONTINUE
71  ei2 = min(ei2,0.999999d0)
72 *
73 * Define unit vectors for inner eccentricity and angular momenta.
74  cosj = 0.0
75  sjsg = 0.0
76  DO 15 k = 1,3
77  ei(k) = ei(k)/sqrt(ei2)
78  hi(k) = a1(k)/sqrt(a12)
79  ho(k) = a2(k)/sqrt(a22)
80  cosj = cosj + hi(k)*ho(k)
81  sjsg = sjsg + ei(k)*ho(k)
82  15 CONTINUE
83 *
84 * Evaluate the expressions A & Z.
85  a = cosj*sqrt(1.0 - ei2)
86  z = (1.0 - ei2)*(2.0 - cosj**2) + 5.0*ei2*sjsg**2
87 *
88 * Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
89  z2 = z**2 + 25.0 + 16.0*a**4 - 10.0*z - 20.0*a**2 - 8.0*a**2*z
90  emax = one6*(z + 1.0 - 4.0*a**2 + sqrt(z2))
91  emax = max(emax,0.0001d0)
92  emax = sqrt(emax)
93 * Skip max eccentricity below 0.90.
94  IF (emax.LT.0.0) go to 50
95 *
96 * Form minimum eccentricity (Douglas Heggie, Sept. 1996).
97  az = a**2 + z - 2.0
98  IF (az.GE.0.0) THEN
99  az1 = 1.0 + z - 4.0*a**2
100  emin2 = one6*(az1 - sqrt(az1**2 - 12.0*az))
101  ELSE
102  emin2 = 1.0 - 0.5*(a**2 + z)
103  END IF
104  emin2 = max(emin2,0.0001d0)
105  emin = sqrt(emin2)
106 *
107 * Estimate Kozai time-scale (N-body units).
108  tk = twopi*semi*sqrt(semi/body(i1))
109  tk1 = twopi*abs(semi1)*sqrt(abs(semi1)/body(i))
110  tkoz = tk1**2*body(i)*(1.0 - ecc1**2)**1.5/(body(i2)*tk)
111 *
112 * Evaluate numerical precession factor in Myr (elliptic integral).
113  const = pfac(a,z)
114 * Note published factor 2/(3*pi) (Kiseleva, Eggleton & Mikkola MN/98).
115  const = const*4.0/(1.5*twopi*sqrt(6.0))
116  tkoz = const*tkoz*tstar
117 *
118 * Obtain inclination (see routine INCLIN).
119  fac = min(cosj,1.0d0)
120  angle = acos(fac)
121  zi = angle*360.0/twopi
122 *
123 * Check optional output (every 100 orbits or each time for #42 > 1).
124  it = it + 1
125  IF ((mod(it,100).EQ.0.AND.emax.GT.0.99).OR.kz(42).GT.1) THEN
126  WRITE (42,30) time+toff, name(i1), zi, emin, emax, ecc,
127  & semi1/semi, tkoz, tmdis(im)
128  30 FORMAT (' KOZAI T NAM1 INC EM EX E A1/A TKOZ TMD ',
129  & f9.3,i7,f7.1,f6.2,2f8.4,f6.1,1p,2e9.1)
130  CALL flush(42)
131  IF (it.GT.2000000000) it = 0
132  END IF
133 *
134 * Search for circularization candidates with large EMAX & inclination.
135  IF (emax.GT.0.90.AND.zi.GT.40.0) THEN
136  CALL findj(i1,jg,im2)
137 * Include safety test to avoid negative JG or JG > N.
138  IF (jg.LT.0.OR.jg.GT.n) go to 50
139 * Consider termination condition for tidal circularization (#27 = 2).
140  IF (semi*(1.0 - emax).LT.2.0*(radius(i1) + radius(jg)).AND.
141  & kz(27).EQ.2) THEN
142 *
143 * Estimate circularization time in Myr (C. Tout, Cambody lecture 2006).
144  IF (radius(jg).GT.radius(i1)) THEN
145  j1 = jg
146  j2 = i1
147  q1 = cm(2,im)/cm(1,im)
148  ELSE
149  j1 = i1
150  j2 = jg
151  q1 = cm(1,im)/cm(2,im)
152  END IF
153 * SEMI0 = SEMI*(1.0 - ECC2)
154  semi0 = semi*(1.0 - emax**2)
155  tcirc = 2.0*q1**2/(1.0 + q1)*(semi0/radius(j1))**8
156  tcirc = 1.0d-06*tcirc
157 * WRITE (45,40) TIME+TOFF, NAME(I1), ECC, SEMI, SEMI0,
158 * & TCIRC
159 * 40 FORMAT (' CIRCULARIZED T NM E A0 A TC ',
160 * & F9.2,I6,F8.4,1P,4E10.2)
161 * Note possibility of switching if nothing happens.
162  rp = semi*(1.0 - emax)
163  tt = min(tkoz,tcirc)
164  WRITE (6,45) emax, semi, rp, radius(i1)+radius(jg), tt
165  45 FORMAT (' KOZAI TIDE EX A RP R1+R2 TT ',
166  & f10.6,1p,e12.4,3e10.2)
167 * ITERM = 1
168 * Perform circularization on small Kozai period or circ time.
169  IF (tkoz.LT.1.0.OR.tcirc.LT.1) iterm = 2
170  tev0(i) = time
171  tev(i) = time + 1.0
172  END IF
173  END IF
174 *
175  50 RETURN
176 *
177  END