Nbody6
stabl4.f
Go to the documentation of this file.
1  SUBROUTINE stabl4(ITERM)
2 *
3 *
4 * Stability test of four-body system.
5 * -----------------------------------
6 *
7  IMPLICIT REAL*8 (a-h,m,o-z)
8  LOGICAL switch,gtype,gtype0,aux
9  common/creg/ m(4),x(12),xd(12),p(12),q(12),time4,energy,epsr2,
10  & xr(9),w(9),r(6),ta(6),mij(6),cm(10),rmax4,tmax,
11  & ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
12  common/tpr/ switch,gtype,gtype0
13  common/ksave/ k1,k2
14  REAL*8 rc(3),vc(3),rc0(3),vc0(3)
15 *
16 *
17 * Transform to physical variables (retain value of SWITCH).
18  aux = switch
19  switch = .false.
20  CALL endreg
21  switch = aux
22 *
23 * Specify indices of two least dominant bodies (denoted K3 & K4).
24  k3 = 0
25  DO 1 l = 1,4
26  IF (l.EQ.k1.OR.l.EQ.k2) go to 1
27  IF (k3.EQ.0) THEN
28  k3 = l
29  ELSE
30  k4 = l
31  END IF
32  1 CONTINUE
33 *
34 * Initialize scalars for orbital elements.
35  rb0 = 0.0d0
36  rb = 0.0d0
37  rb1 = 0.0d0
38  rb2 = 0.0d0
39  rdot = 0.0d0
40  rdot1 = 0.0d0
41  rdot2 = 0.0d0
42  vrel2 = 0.0d0
43  vrel20 = 0.0d0
44  vrel21 = 0.0d0
45  vrel22 = 0.0d0
46 *
47 * Define binary masses of smallest & widest pair (K1 & K2 and K3 & K4).
48  mb0 = m(k1) + m(k2)
49  mb = m(k3) + m(k4)
50 *
51 * Form separations & velocities of MB0, MB and their relative orbit.
52  DO 10 k = 1,3
53  j1 = 3*(k1-1) + k
54  j2 = 3*(k2-1) + k
55  j3 = 3*(k3-1) + k
56  j4 = 3*(k4-1) + k
57  rc0(k) = (m(k1)*x(j1) + m(k2)*x(j2))/mb0
58  rc(k) = (m(k3)*x(j3) + m(k4)*x(j4))/mb
59  vc0(k) = (m(k1)*xd(j1) + m(k2)*xd(j2))/mb0
60  vc(k) = (m(k3)*xd(j3) + m(k4)*xd(j4))/mb
61  rb0 = rb0 + (x(j1) - x(j2))**2
62  rb = rb + (x(j3) - x(j4))**2
63  rb1 = rb1 + (rc(k) - rc0(k))**2
64  rb2 = rb2 + (rc0(k) - x(j3))**2
65  rdot = rdot + (x(j3) - x(j4))*(xd(j3) - xd(j4))
66  rdot1 = rdot1 + (rc(k) - rc0(k))*(vc(k) - vc0(k))
67  rdot2 = rdot2 + (rc0(k) - x(j3))*(vc0(k) - xd(j3))
68  vrel2 = vrel2 + (xd(j3) - xd(j4))**2
69  vrel20 = vrel20 + (xd(j1) - xd(j2))**2
70  vrel21 = vrel21 + (vc(k) - vc0(k))**2
71  vrel22 = vrel22 + (vc0(k) - xd(j3))**2
72  10 CONTINUE
73 *
74 * Determine semi-major axis of inner binary.
75  rb0 = sqrt(rb0)
76  semi0 = 2.0/rb0 - vrel20/mb0
77  semi0 = 1.0/semi0
78  e0 = 1.0 - rb0/semi0
79 *
80 * Form semi-major axis & eccentricity of outer pair.
81  rb = sqrt(rb)
82  semi = 2.0d0/rb - vrel2/mb
83  semi = 1.0/semi
84 * E = SQRT((1.0D0 - RB/SEMI)**2 + RDOT**2/(SEMI*MB))
85 *
86 * Evaluate orbital elements of relative c.m. motion.
87  rb1 = sqrt(rb1)
88  semi1 = 2.0d0/rb1 - vrel21/cm(7)
89  semi1 = 1.0/semi1
90  e1 = sqrt((1.0d0 - rb1/semi1)**2 + rdot1**2/(semi1*cm(7)))
91 *
92 * Consider the inner triple.
93  mb2 = mb0 + m(k3)
94  rb2 = sqrt(rb2)
95  semi2 = 2.0d0/rb2 - vrel22/mb2
96  semi2 = 1.0/semi2
97  e2 = sqrt((1.0d0 - rb2/semi2)**2 + rdot2**2/(semi2*mb2))
98 *
99 * Obtain standard stability ratio (outer pericentre / inner apocentre).
100  ratio = semi1*(1.0d0 - e1)/(semi0*(1.0d0 + e0))
101 *
102 * Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
103 * AM = (2.65 + E0)*(1.0 + MB0/MB)**0.3333
104 * FM = (2.0*MB0 - MB)/(3.0*MB)
105 *
106 * Expand natural logarithm for small arguments.
107 * IF (ABS(FM).LT.0.67) THEN
108 * BM = FM*(1.0 - (0.5 - 0.3333*FM)*FM)
109 * ELSE
110 * BM = LOG(1.0D0 + FM)
111 * END IF
112 *
113 * Adopt mass dependent criterion of Harrington (A.J. 80) & Bailyn.
114 * PCRIT = AM*(1.0 + 0.7*BM)*SEMI0
115 *
116 * Form hierarchical stability ratio (Kiseleva & Eggleton 1995).
117 * Q0 = MB/MB0
118 * Q1 = MAX(M(K2)/M(K1),M(K1)/M(K2))
119 * Q3 = Q0**0.33333
120 * Q13 = Q1**0.33333
121 * AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
122 * PCRIT = AR*SEMI0*(1.0D0 + E0)
123 *
124 * Check stability (AM 1997; inner triple or well separated quadruple).
125  iterm = 0
126  IF (rb1.GT.5.0*rb2.AND.e2.LT.1.0) THEN
127  q1 = m(k3)/mb0
128  xfac = (1.0 + q1)*(1.0 + e2)/sqrt(1.0 - e2)
129  pcrit = 2.8*xfac**0.4*semi0
130  pmin = semi2*(1.0 - e2)
131  IF (pcrit.LT.pmin) THEN
132  iterm = -1
133  ratio = semi2*(1.0d0 - e2)/(semi0*(1.0d0 + e0))
134  WRITE (6,15) semi0, semi2, e0, e2, ratio, rb0, rb2,
135  & pcrit, pmin
136  15 FORMAT (' STABT: A0 A2 E0 E2 RATIO R0 R2 PCR PM ',
137  & 1p,2e10.2,0p,2f7.3,f6.2,1p,4e9.1)
138  END IF
139  ELSE IF (rb1.GT.5.0*max(rb0,rb).AND.e1.LT.1.0.AND.
140  & min(semi0,semi).GT.0.0) THEN
141 * Choose smallest binary as third body and ignore fudge factor.
142  IF (semi.GT.semi0) THEN
143  q1 = mb0/mb
144  ain = semi
145  ELSE
146  q1 = mb/mb0
147  ain = semi0
148  END IF
149  xfac = (1.0 + q1)*(1.0 + e1)/sqrt(1.0 - e1)
150  pcrit = 2.8*xfac**0.4*ain
151  pmin = semi1*(1.0 - e1)
152  IF (pcrit.LT.pmin) THEN
153  iterm = -1
154  WRITE (6,20) ain, semi1, e0, e1, ratio, rb0, rb1,
155  & pcrit, pmin
156  20 FORMAT (' STABQ: AIN A1 E0 E1 RATIO R0 R1 PCR PM ',
157  & 1p,2e10.2,0p,2f7.3,f6.2,1p,4e9.1)
158  END IF
159  END IF
160 *
161  RETURN
162 *
163  END