Nbody6
stabl3.f
Go to the documentation of this file.
1  SUBROUTINE stabl3(ITERM)
2 *
3 *
4 * Stability test of three-body system.
5 * ------------------------------------
6 *
7  IMPLICIT REAL*8 (a-h,m,o-z)
8  common/azreg/ time3,tmax,q(8),p(8),r1,r2,r3,energy,m(3),x3(3,3),
9  & xdot3(3,3),cm(10),c11,c12,c19,c20,c24,c25,
10  & nstep3,name3(3),kz15,kz27
11  common/azout/ mb,rb,ri,semi,semi1,e,e1
12 *
13 *
14 * Transform to physical variables.
15  CALL trans3(3)
16 *
17 * Identify index of second binary component and the third body.
18  im = 1
19  IF (r2.LT.r1) im = 2
20  i = 3 - im
21 *
22 * Form scalar product of distance & velocity for body #I.
23  ridot = x3(1,i)*xdot3(1,i) + x3(2,i)*xdot3(2,i) +
24  & x3(3,i)*xdot3(3,i)
25 *
26 * Set distance & radial velocity of body #I with respect to binary.
27  mb = m(im) + m(3)
28  fac = cm(7)/mb
29  ri = sqrt(x3(1,i)**2 + x3(2,i)**2 + x3(3,i)**2)
30  ridot = fac*ridot/ri
31  ri = fac*ri
32 *
33 * Evaluate orbital elements.
34  rdot = 0.0d0
35  vrel2 = 0.0d0
36  vi2 = 0.0d0
37  DO 10 k = 1,3
38  rdot = rdot +
39  & (x3(k,3) - x3(k,im))*(xdot3(k,3) - xdot3(k,im))
40  vrel2 = vrel2 + (xdot3(k,3) - xdot3(k,im))**2
41  vi2 = vi2 + xdot3(k,i)**2
42  10 CONTINUE
43 *
44 * Determine semi-major axis & eccentricity of inner binary.
45  rb = min(r1,r2)
46  semi = 2.0d0/rb - vrel2/mb
47  semi = 1.0/semi
48  e = sqrt((1.0d0 - rb/semi)**2 + rdot**2/(semi*mb))
49 *
50 * Form outer semi-major axis & eccentricity.
51  vi2 = vi2*fac**2
52  semi1 = 2.0d0/ri - vi2/cm(7)
53  semi1 = 1.0/semi1
54  e1 = sqrt((1.0d0 - ri/semi1)**2 + (ri*ridot)**2/(semi1*cm(7)))
55 *
56 * Obtain standard stability ratio (outer pericentre / inner apocentre).
57  ratio = semi1*(1.0d0 - e1)/(semi*(1.0d0 + e))
58 *
59 * Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
60 * AM = (2.65 + E)*(1.0 + M(I)/MB)**0.3333
61 * FM = (2.0*M(I) - MB)/(3.0*MB)
62 *
63 * Expand natural logarithm for small arguments.
64 * IF (ABS(FM).LT.0.67) THEN
65 * BM = FM*(1.0 - (0.5 - ONE3*FM)*FM)
66 * ELSE
67 * BM = LOG(1.0D0 + FM)
68 * END IF
69 * Define mass dependent criterion of Harrington (A.J. 80) & Bailyn.
70 * PCRIT = AM*(1.0 + 0.7*BM)*SEMI
71 *
72 * Form hierarchical stability ratio (Kiseleva & Eggleton 1995).
73 * Q0 = MB/M(I)
74 * Q1 = MAX(M(3)/M(IM),M(IM)/M(3))
75 * Q3 = Q0**0.33333
76 * Q13 = Q1**0.33333
77 * AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
78 * PCRIT = AR*SEMI*(1.0D0 + E)
79 *
80 * Adopt the semi-analytical stability criterion (MA 1997).
81  q1 = m(i)/mb
82  IF (e1.LT.1.0) THEN
83  xfac = (1.0 + q1)*(1.0 + e1)/sqrt(1.0 - e1)
84  ELSE
85  xfac = 40.0*(1.0 + q1)
86  END IF
87  pcrit = 2.8*xfac**0.4*semi
88  pmin = semi1*(1.0d0 - e1)
89 *
90 * Set negative termination index if system is stable and RB > SEMI.
91  IF (pcrit.LT.pmin.AND.e1.LT.1.0.AND.1./rb.LT.1./semi) THEN
92  iterm = -1
93 * Obtain Zare's stability parameter (valid for small inclinations).
94  m1 = m(im)
95  m2 = m(3)
96  m3 = m(i)
97  CALL stablz(m1,m2,m3,sp)
98  WRITE (6,20) semi, semi1, e, e1, ri, ratio, sp, pcrit, pmin
99  20 FORMAT (' STABT: A A1 E E1 RI RATIO SP PCR PM ',
100  & 1p,2e10.2,0p,2f6.2,f9.5,2f6.2,1p,2e9.1)
101 * Terminate if escaper is outside 3*SEMI (> 0).
102  ELSE IF (e1.GT.1.0.AND.1./ri.LT.1./(3.0*semi)) THEN
103  iterm = -1
104  ELSE
105  iterm = 0
106  END IF
107 *
108  RETURN
109 *
110  END