Nbody6
histab.f
Go to the documentation of this file.
1  SUBROUTINE histab(IPAIR,J,PMIN,RSTAB)
2 *
3 *
4 * Hierarchical stability condition.
5 * ---------------------------------
6 *
7  include 'common6.h'
8  REAL*8 xx(3,3),vv(3,3)
9 *
10 *
11 * Define KS & c.m. indices and save IPAIR for RETURN.
12  i1 = 2*ipair - 1
13  i2 = i1 + 1
14  i = n + ipair
15  jpair = ipair
16 *
17 * Evaluate semi-major axis & eccentricity of inner binary.
18  semi = -0.5d0*body(i)/h(ipair)
19  ecc2 = (1.0d0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
20  ecc = sqrt(ecc2)
21 *
22 * Form stability ratio (Eggleton & Kiseleva, Ap.J. 455, 640, 1995).
23 * Q = BODY(I)/BODY(J)
24 * Q1 = MAX(BODY(I2)/BODY(I1),BODY(I1)/BODY(I2))
25 * Q3 = Q**0.33333
26 * Q13 = Q1**0.33333
27 * Note sign error in third term of Eq. (2).
28 * AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
29 *
30 * Save hierarchical stability separation in COMMON variable.
31 * PCRIT = AR*SEMI*(1.0D0 + ECC)
32 *
33 * Determine the outer eccentricity.
34  rij2 = 0.0
35  vij2 = 0.0
36  rdot = 0.0
37  DO 5 k = 1,3
38  rij2 = rij2 + (x(k,i) - x(k,j))**2
39  vij2 = vij2 + (xdot(k,i) - xdot(k,j))**2
40  rdot = rdot + (x(k,i) - x(k,j))*(xdot(k,i) - xdot(k,j))
41  5 CONTINUE
42  rij = sqrt(rij2)
43  a1 = 2.0/rij - vij2/(body(i) + body(j))
44 *
45 * Exit on hyperbolic orbit.
46  IF (a1.LT.0.0) THEN
47  pmin = 1.0
48  rstab = 0.0
49  go to 20
50  END IF
51 *
52  semi1 = 1.0/a1
53  ecc2 = (1.0 - rij/semi1)**2 + rdot**2/(semi1*(body(i) + body(j)))
54  ecc1 = sqrt(ecc2)
55  pmin = semi1*(1.0 - ecc1)
56 *
57 * Evaluate the basic stability condition without fudge factor.
58  q = body(j)/body(i)
59  IF (ecc1.LT.1.0) THEN
60  xfac = (1.0 + q)*(1.0 + ecc1)/sqrt(1.0 - ecc1)
61  ELSE
62  xfac = 40.0*(1.0 + q)
63  END IF
64  pcrit = 2.8*xfac**0.4*semi
65 *
66 * Exit if stability value falls outside practical limits.
67  IF ((pcrit.GT.1.5*pmin.OR.pcrit.LT.0.5*pmin).AND.j.LE.n) THEN
68  rstab = pcrit
69  go to 20
70  END IF
71 *
72 * Choose the most active triple in case of two binaries.
73  jj = j
74  IF (j.GT.n) THEN
75  semi2 = -0.5d0*body(j)/h(j-n)
76 * Adopt 10% fudge factor with linear dependence on smallest ratio.
77  yfac = 1.0 + 0.1*min(semi2/semi,semi/semi2)
78  IF (semi2.GT.semi) THEN
79  ecc2 = (1.0 - r(j-n)/semi2)**2 +
80  & tdot2(j-n)**2/(body(j)*semi2)
81  ecc = sqrt(ecc2)
82  semiz = semi2
83  semi2 = semi
84  semi = semiz
85  ipair = j - n
86  i1 = 2*ipair - 1
87  i2 = i1 + 1
88  jj = i
89  i = j
90  END IF
91  ELSE
92  yfac = 1.0
93  END IF
94 *
95 * Resolve weakly perturbed binary (prevent X(K,I1) = X(K,I2)).
96  IF (gamma(ipair).LT.gmin.OR.x(1,i1).EQ.x(1,i2)) THEN
97  CALL resolv(ipair,1)
98  END IF
99 *
100 * Determine inclination for triple configuration (NB! radians).
101  DO 10 k = 1,3
102  xx(k,1) = x(k,i1)
103  xx(k,2) = x(k,i2)
104  xx(k,3) = x(k,jj)
105  vv(k,1) = xdot(k,i1)
106  vv(k,2) = xdot(k,i2)
107  vv(k,3) = xdot(k,jj)
108  10 CONTINUE
109  CALL inclin(xx,vv,x(1,i),xdot(1,i),angle)
110 *
111 * Employ the improved stability criterion for doubtful cases.
112  rstab = stability(body(i1),body(i2),body(jj),ecc,ecc1,angle)*semi
113 * Note: the present stability routine includes inclination!
114 * RSTAB = YFAC*RSTAB
115  pcrit = rstab
116  ipair = jpair
117 *
118  20 RETURN
119 *
120  END