Nbody6
 All Files Functions Variables
himax.f
Go to the documentation of this file.
1  SUBROUTINE himax(I,IM,ECC,SEMI,EMAX,EMIN,ZI,TG,EDAV)
2 *
3 *
4 * Maximum eccentricity of hierarchical binary.
5 * --------------------------------------------
6 *
7  include 'common6.h'
8  common/binary/ cm(4,mmax),yrel(3,mmax),zrel(3,mmax),
9  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
10  & namem(mmax),nameg(mmax),kstarm(mmax),iflag(mmax)
11  REAL*8 a1(3),a2(3),xrel(3),vrel(3),ei(3),hi(3),ho(3),bhat(3),
12  & ui(4),v(4)
13 *
14 *
15 * Copy KS variables to local scalars.
16  DO 1 k = 1,4
17  ui(k) = um(k,im)
18  v(k) = umdot(k,im)
19  1 CONTINUE
20 *
21 * Transform to physical variables and multiply by 4 (momentum formula).
22  CALL ksphys(ui,v,xrel,vrel)
23  DO 2 k = 1,3
24  vrel(k) = 4.0*vrel(k)
25  2 CONTINUE
26 *
27 * Specify index for outer body (second KS component) & pair index.
28  jcomp = i + 1
29  ipair = kvec(i)
30 *
31 * Ensure that unperturbed binary is resolved.
32  IF (list(1,2*ipair-1).EQ.0.OR.x(1,jcomp).EQ.x(1,i)) THEN
33  CALL resolv(ipair,1)
34  END IF
35 *
36 * Obtain inner & outer angular momentum by cyclic notation.
37  a12 = 0.0
38  a22 = 0.0
39  a1a2 = 0.0
40  ri2 = 0.0
41  vi2 = 0.0
42  rvi = 0.0
43  DO 5 k = 1,3
44  k1 = k + 1
45  IF (k1.GT.3) k1 = 1
46  k2 = k1 + 1
47  IF (k2.GT.3) k2 = 1
48  a1(k) = xrel(k1)*vrel(k2) - xrel(k2)*vrel(k1)
49  a2(k) = (x(k1,jcomp) - x(k1,i))*(xdot(k2,jcomp) - xdot(k2,i))
50  & - (x(k2,jcomp) - x(k2,i))*(xdot(k1,jcomp) - xdot(k1,i))
51  a12 = a12 + a1(k)**2
52  a22 = a22 + a2(k)**2
53  a1a2 = a1a2 + a1(k)*a2(k)
54  ri2 = ri2 + xrel(k)**2
55  vi2 = vi2 + vrel(k)**2
56  rvi = rvi + xrel(k)*vrel(k)
57  5 CONTINUE
58 *
59 * Evaluate orbital parameters for outer orbit from KS elements.
60  zmb = body(i) + body(jcomp)
61  semi1 = -0.5*zmb/h(ipair)
62  ecc2 = (1.0 - r(ipair)/semi1)**2 + tdot2(ipair)**2/(semi1*zmb)
63  ecc1 = sqrt(ecc2)
64 * Determine inclination in radians.
65  fac = a1a2/sqrt(a12*a22)
66  zi = acos(fac)
67 *
68 * Construct the Runge-Lenz vector (Heggie & Rasio 1995, IAU174, Eq.5).
69  ei2 = 0.0
70  DO 10 k = 1,3
71  ei(k) = (vi2*xrel(k) - rvi*vrel(k))/body(i) -
72  & xrel(k)/sqrt(ri2)
73  ei2 = ei2 + ei(k)**2
74  10 CONTINUE
75 *
76 * Define unit vectors for inner eccentricity and angular momenta.
77  cosj = 0.0
78  sjsg = 0.0
79  DO 15 k = 1,3
80  ei(k) = ei(k)/sqrt(ei2)
81  hi(k) = a1(k)/sqrt(a12)
82  ho(k) = a2(k)/sqrt(a22)
83  cosj = cosj + hi(k)*ho(k)
84  sjsg = sjsg + ei(k)*ho(k)
85  15 CONTINUE
86 *
87 * Form unit vector BHAT and scalars AH & BH (Douglas Heggie, 10/9/96).
88  ah = 0.0
89  bh = 0.0
90  DO 16 k = 1,3
91  k1 = k + 1
92  IF (k1.GT.3) k1 = 1
93  k2 = k1 + 1
94  IF (k2.GT.3) k2 = 1
95  bhat(k) = hi(k1)*ei(k2) - hi(k2)*ei(k1)
96  ah = ah + ei(k)*ho(k)
97  bh = bh + bhat(k)*ho(k)
98  16 CONTINUE
99 *
100 * Evaluate the expressions A & Z.
101  a = cosj*sqrt(1.0 - ei2)
102  z = (1.0 - ei2)*(2.0 - cosj**2) + 5.0*ei2*sjsg**2
103 *
104 * Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
105  z2 = z**2 + 25.0 + 16.0*a**4 - 10.0*z - 20.0*a**2 - 8.0*a**2*z
106  emax = one6*(z + 1.0 - 4.0*a**2 + sqrt(z2))
107  emax = sqrt(emax)
108 *
109 * Form minimum eccentricity (Douglas Heggie, Sept. 1996).
110  az = a**2 + z - 2.0
111  IF (az.GE.0.0) THEN
112  az1 = 1.0 + z - 4.0*a**2
113  emin2 = one6*(az1 - sqrt(az1**2 - 12.0*az))
114  ELSE
115  emin2 = 1.0 - 0.5*(a**2 + z)
116  END IF
117  emin2 = max(emin2,0.0d0)
118  emin = sqrt(emin2)
119 *
120 * Estimate eccentricity growth time-scale.
121  tk = twopi*semi*sqrt(semi/body(i))
122  tk1 = twopi*abs(semi1)*sqrt(abs(semi1)/zmb)
123  tg = tk1**2*zmb*(1.0 - ecc1**2)**1.5/(body(jcomp)*tk)
124 *
125 * Evaluate numerical precession factor (involves elliptic integral).
126  const = pfac(a,z)
127  const = const*4.0/(1.5*twopi*sqrt(6.0))
128 *
129 * Convert growth time to units of 10**6 yrs.
130  tg = const*tg*tstar
131 *
132 * Form doubly averaged eccentricity derivative (Douglas Heggie 9/96).
133  yfac = 15.0*body(jcomp)/(4.0*zmb)*twopi*tk/tk1**2
134  yfac = yfac*ecc*sqrt(1.0 - ecc**2)/(1.0 - ecc1**2)**(1.5)
135  edav = yfac*ah*bh
136 *
137  RETURN
138 *
139  END