Nbody6
 All Files Functions Variables
himax2.f
Go to the documentation of this file.
1  SUBROUTINE himax2(I1,ECC,SEMI,ECC1,SEMI1,EMAX,EMIN,ZI,TG,EDAV)
2 *
3 *
4 * Maximum eccentricity of inner 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 * Specify pair index of hierarchy and its c.m. index.
12  i2 = i1 + 1
13  ipair = kvec(i1)
14  icm = n + ipair
15 *
16 * Obtain global coordinates (routine KSINT uses local values).
17  CALL resolv(ipair,1)
18 *
19 * Form relative coordinates & velocities of current binary.
20  DO 2 k = 1,3
21  xrel(k) = x(k,i1) - x(k,i2)
22  vrel(k) = xdot(k,i1) - xdot(k,i2)
23  2 CONTINUE
24 *
25 * Obtain inner & outer angular momentum by cyclic notation.
26  a12 = 0.0
27  a22 = 0.0
28  a1a2 = 0.0
29  ri2 = 0.0
30  vi2 = 0.0
31  rvi = 0.0
32  DO 5 k = 1,3
33  k1 = k + 1
34  IF (k1.GT.3) k1 = 1
35  k2 = k1 + 1
36  IF (k2.GT.3) k2 = 1
37  a1(k) = xrel(k1)*vrel(k2) - xrel(k2)*vrel(k1)
38  a2(k) = (x(k1,jcomp)-x(k1,icm))*(xdot(k2,jcomp)-xdot(k2,icm))
39  & - (x(k2,jcomp)-x(k2,icm))*(xdot(k1,jcomp)-xdot(k1,icm))
40  a12 = a12 + a1(k)**2
41  a22 = a22 + a2(k)**2
42  a1a2 = a1a2 + a1(k)*a2(k)
43  ri2 = ri2 + xrel(k)**2
44  vi2 = vi2 + vrel(k)**2
45  rvi = rvi + xrel(k)*vrel(k)
46  5 CONTINUE
47 *
48 * Evaluate orbital parameters for inner outer orbit from KS elements.
49  zmb = body(i1) + body(i2)
50 * Determine inclination in radians.
51  fac = a1a2/sqrt(a12*a22)
52  zi = acos(fac)
53 *
54 * Construct the Runge-Lenz vector (Heggie & Rasio 1995, IAU174, Eq.5).
55  ei2 = 0.0
56  DO 10 k = 1,3
57  ei(k) = (vi2*xrel(k) - rvi*vrel(k))/body(icm) -
58  & xrel(k)/sqrt(ri2)
59  ei2 = ei2 + ei(k)**2
60  10 CONTINUE
61 *
62 * Define unit vectors for inner eccentricity and angular momenta.
63  cosj = 0.0
64  sjsg = 0.0
65  DO 15 k = 1,3
66  ei(k) = ei(k)/sqrt(ei2)
67  hi(k) = a1(k)/sqrt(a12)
68  ho(k) = a2(k)/sqrt(a22)
69  cosj = cosj + hi(k)*ho(k)
70  sjsg = sjsg + ei(k)*ho(k)
71  15 CONTINUE
72 *
73 * Form unit vector BHAT and scalars AH & BH (Douglas Heggie, 10/9/96).
74  ah = 0.0
75  bh = 0.0
76  DO 16 k = 1,3
77  k1 = k + 1
78  IF (k1.GT.3) k1 = 1
79  k2 = k1 + 1
80  IF (k2.GT.3) k2 = 1
81  bhat(k) = hi(k1)*ei(k2) - hi(k2)*ei(k1)
82  ah = ah + ei(k)*ho(k)
83  bh = bh + bhat(k)*ho(k)
84  16 CONTINUE
85 *
86 * Evaluate the expressions A & Z.
87  a = cosj*sqrt(1.0 - ei2)
88  z = (1.0 - ei2)*(2.0 - cosj**2) + 5.0*ei2*sjsg**2
89 *
90 * Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
91  z2 = z**2 + 25.0 + 16.0*a**4 - 10.0*z - 20.0*a**2 - 8.0*a**2*z
92  emax = one6*(z + 1.0 - 4.0*a**2 + sqrt(z2))
93  emax = sqrt(emax)
94 *
95 * Form minimum eccentricity (Douglas Heggie, Sept. 1996).
96  az = a**2 + z - 2.0
97  IF (az.GE.0.0) THEN
98  az1 = 1.0 + z - 4.0*a**2
99  emin2 = one6*(az1 - sqrt(az1**2 - 12.0*az))
100  ELSE
101  emin2 = 1.0 - 0.5*(a**2 + z)
102  END IF
103  emin2 = max(emin2,0.0d0)
104  emin = sqrt(emin2)
105 *
106 * Estimate eccentricity growth time-scale.
107  zmb2 = zmb + body(jcomp)
108  tk = twopi*semi*sqrt(semi/zmb)
109  tk1 = twopi*abs(semi1)*sqrt(abs(semi1)/zmb2)
110  tg = tk1**2*zmb2*(1.0 - ecc1**2)**1.5/(body(jcomp)*tk)
111 *
112 * Evaluate numerical precession factor (involves elliptic integral).
113  const = pfac(a,z)
114  const = const*4.0/(1.5*twopi*sqrt(6.0))
115 *
116 * Convert growth time to units of 10**6 yrs.
117  tg = const*tg*tstar
118 *
119 * Form doubly averaged eccentricity derivative (Douglas Heggie 9/96).
120  yfac = 15.0*body(jcomp)/(4.0*zmb2)*twopi*tk/tk1**2
121  yfac = yfac*ecc*sqrt(1.0 - ecc**2)/(1.0 - ecc1**2)**(1.5)
122  edav = yfac*ah*bh
123 *
124  RETURN
125 *
126  END