Nbody6
 All Files Functions Variables
nstab.f
Go to the documentation of this file.
1 c general three-body stability algorithm
2 c
3 c system is unstable if nstab=1 is returned
4 c system is stable if nstab=0 is returned
5 c
6 c Rosemary Mardling
7 c School of Mathematical Sciences, Monash University
8 c
9 c version as of 16-11-07
10 c email mardling@sci.monash.edu.au to be added to updates email list
11 c preprint on astro-ph by New Year :-)
12 c
13 c sigma=period ratio (outer/inner) **should be > 1**
14 c ai=inner semi-major axis
15 c a0=outer semi-major axis
16 c ei0=initial inner eccentricity
17 c eo=outer eccentricity
18 c relinc=relative inclination (radians)
19 c m1, m2, m3=masses (any units; m3=outer body)
20 c
21 c valid for all inclinations
22 c
23 c MASS RATIO CONDITIONS
24 c valid for systems with at least one of m_2/m_1>0.05 OR m_3/m_1>0.05
25 c (so that one could have, for example, m_2/m_1=0 and m_3/m_1=0.1)
26 c OR BOTH m2/m1>0.01 AND m3/m1>0.01
27 c **future version will include other resonances to cover smaller mass ratios
28 c
29 c assumes resonance angle phi=0 because resonance overlap criterion doesn't recognize
30 c instability outside separatrix.
31 c
32 c system is unstable if nstab=1 is returned
33 c system is stable if nstab=0 is returned
34 
35  integer function nstab(ai,a0,ei0,eo,relinc,m1,m2,m3)
36 
37  implicit real*8 (a-h,m,o-z)
38  common/params2/mm1,mm2,mm3
39  common/savepi/pi
40  save itime
41  data itime/0/
42 
43 c reject outer pericentre inside inner apocentre
44  if(a0*(1.-eo).lt.ai*(1.+ei0))then
45  nstab=1
46  return
47  endif
48 
49  if(itime.eq.0)then
50  pi=4.d0*datan(1.0d0)
51  itime=1
52  endif
53 
54  mm1=m1
55  mm2=m2
56  mm3=m3
57 
58  m12=m1+m2
59  m123=m12+m3
60 
61 c set period ratio (outer/inner)
62  a0ai=a0/ai
63  sigma=sqrt(a0ai**3*m12/m123)
64 c do not allow period ratio < 1
65  if(sigma.lt.1.0)then
66  nstab=1
67  return
68  endif
69 
70  mi2=m3/m123
71  mo2=(m1*m2/m12**2)*(m12/m123)**(2./3.)
72  mi3=(m3/m12)*(m12/m123)**(4./3.)*(m1-m2)/m12
73  mo3=(m1*m2/m12**2)*(m12/m123)*(m1-m2)/m12
74 
75  c22=3./8.
76  c20=0.25
77  c31=sqrt(3.)/4.
78  c33=-sqrt(5.)/4.
79 
80  e=eo
81 
82 c inclination coefficients
83 
84  win=0
85 
86  a=sqrt(1-ei0**2)*cos(relinc)
87  z=(1-ei0**2)*(1+sin(relinc)**2)+5*ei0**2*
88  & (sin(win)*sin(relinc))**2
89  del=z**2+25+16*a**4-10*z-20*a**2-8*a**2*z
90 
91  ek=sqrt(abs((z+1-4*a**2+sqrt(del))/6.))
92  cosik=a/sqrt(1-ek**2)
93  sinik=sqrt(1-cosik**2)
94 
95  gam222=0.25*(1+cosik)**2
96  gam22m2=0.25*(1-cosik)**2
97  gam220=0.5*sqrt(1.5)*sinik**2
98  gam200=0.5*(3*cosik**2-1)
99 
100 c induced inner eccentricity
101  ei=ein_induced(sigma,ei0,e,relinc)
102 
103 c octopole emax
104  if(m1.ne.m2)then
105  eoctmax=eoct(sigma,ei0,e)
106  ei=max(eoctmax,ei)
107  endif
108 
109  ei=max(ek,ei)
110  ei=min(ei,1.0d0)
111 
112  n=sigma
113  nstab=0
114 
115 c [n:1](222) resonance
116  s221=-3*ei+(13./8.)*ei**3+(5./192.)*ei**5
117 
118  f22n=flmn(2,2,n,e)/(1-e)**3
119 
120  an=abs(6*c22*s221*f22n*(mi2+mo2*sigma**0.666)*gam222)
121  phi=0
122  en=0.5*(sigma-n)**2-an*(1+cos(phi))
123 
124 c [n+1:1](222) resonance
125  f22n=flmn(2,2,n+1,e)/(1-e)**3
126 
127  an=abs(6*c22*s221*f22n*(mi2+mo2*sigma**0.666)*gam222)
128 
129  enp1=0.5*(sigma-(n+1))**2-an*(1+cos(phi))
130  if(en.lt.0.and.enp1.lt.0)nstab=1
131 
132 c [n:1](22-2) resonance
133  s22m1=-(ei**3*(4480 + 1880*ei**2 + 1091*ei**4))/15360.
134 
135  f22n=flmn(2,2,n,e)/(1-e)**3
136 
137  an=abs(6*c22*s22m1*f22n*(mi2+mo2*sigma**0.666)*gam22m2)
138 
139  phi=0
140  en=0.5*(sigma-n)**2-an*(1+cos(phi))
141 
142 c [n+1:1](22-2) resonance
143  f22n=flmn(2,2,n+1,e)/(1-e)**3
144 
145  an=abs(6*c22*s22m1*f22n*(mi2+mo2*sigma**0.666)*gam22m2)
146 
147  enp1=0.5*(sigma-(n+1))**2-an*(1+cos(phi))
148  if(en.lt.0.and.enp1.lt.0)nstab=1
149 
150 c [n:1](202) resonance
151  s201=(ei*(-9216 + 1152*ei**2 - 48*ei**4 + ei**6))/9216.
152 
153  f22n=flmn(2,2,n,e)/(1-e)**3
154 
155  an=abs(6*sqrt(c20*c22)*s201*f22n*(mi2+mo2*sigma**0.666)*gam220)
156 
157  phi=0
158  en=0.5*(sigma-n)**2-an*(1+cos(phi))
159 
160 c [n+1:1](202) resonance
161  f22n=flmn(2,2,n+1,e)/(1-e)**3
162 
163  an=abs(6*sqrt(c20*c22)*s201*f22n*(mi2+mo2*sigma**0.666)*gam220)
164 
165  enp1=0.5*(sigma-(n+1))**2-an*(1+cos(phi))
166  if(en.lt.0.and.enp1.lt.0)nstab=1
167 
168 c [n:1](002) resonance
169  s201=(ei*(-9216 + 1152*ei**2 - 48*ei**4 + ei**6))/9216.
170 
171  f20n=flmn(2,0,n,e)/(1-e)**3
172 
173  an=abs(3*c20*s201*f20n*(mi2+mo2*sigma**0.666)*gam200)
174 
175  phi=0
176  en=0.5*(sigma-n)**2-an*(1+cos(phi))
177 
178 c [n+1:1](002) resonance
179  f20n=flmn(2,0,n+1,e)/(1-e)**3
180 
181  an=abs(3*c20*s201*f20n*(mi2+mo2*sigma**0.666)*gam200)
182 
183  enp1=0.5*(sigma-(n+1))**2-an*(1+cos(phi))
184  if(en.lt.0.and.enp1.lt.0)nstab=1
185 
186  end
187 
188 c ---------------------------------------------------------------------------
189 c Asymptotic expression for f^(lm)_n(e) for all e<1 and n.
190 c
191  real*8 function flmn(l,m,n,e)
192 
193  implicit real*8 (a-h,o-z)
194  common/savepi/pi
195 
196  if(e.lt.5.e-3)then
197  if(m.eq.n)then
198  flmn=1
199  else
200  flmn=0
201  endif
202  return
203  endif
204 
205  rho=n*(1-e)**1.5
206 
207  xi=(acosh(1/e)-sqrt(1-e**2))/(1-e)**1.5
208 
209  flmn=(1/(2*pi*n))*2.0**m*(sqrt(2*pi)/facfac(l,m))*
210  . ((1+e)**(real(3*m-l-1)/4.)/e**m)*
211  . (rho**(real(l+m+1)/2.))*
212  . exp(-rho*xi)
213 
214  end
215 
216 
217 c ---------------------------------------------------------------------------
218  real*8 function ein_induced(sigma,ei0,e,relinc)
219 
220  implicit real*8 (a-h,m,o-z)
221  common/params2/m1,m2,m3
222  common/savepi/pi
223 
224  m123=m1+m2+m3
225  n=sigma
226 
227  gam222=0.25*(1+cos(relinc))**2
228  gam220=0.5*sqrt(1.5)*sin(relinc)**2
229  gam200=0.5*(3*cos(relinc)**2-1)
230 
231  f22n=flmn(2,2,n,e)/(1-e)**3
232  f20n=flmn(2,0,n,e)/(1-e)**3
233 
234  prod222=f22n*gam222
235  prod220=f22n*gam220
236  prod200=f20n*gam200
237 
238  prod=max(prod222,prod220,prod200)
239 
240  a=4.5*(m3/m123)*(2*pi*n)*prod/sigma**2
241 
242  ein_induced=sqrt(ei0**2+a**2)
243 
244  end
245 c ------------------------------------------------------------------------------
246 c eoct.f
247 c
248 c calculates maximum eccentricity for arbitrary coplanar system
249 c using Mardling (2007) MNRAS in press
250 c
251  real*8 function eoct(sigma,ei0,eo)
252  implicit real*8 (a-h,m,o-z)
253  common/params2/m1,m2,m3
254  common/savepi/pi
255 
256  m12=m1+m2
257  m123=m12+m3
258  aoai=((m123/m12)*sigma**2)**0.3333
259  al=1/aoai
260 
261  epso=sqrt(1-eo**2)
262 
263  eeq=1.25*al*eo/epso**2/abs(1-sqrt(al)*(m2/m3)/epso)
264 
265  aa=abs(1-ei0/eeq)
266 
267  if(aa.lt.1)then
268  eoct=(1+aa)*eeq
269  else
270  eoct=ei0+2*eeq
271  endif
272 
273  end
274 
275 c ---------------------------------------------------------------------------
276  real*8 function acosh(x)
277  real*8 x
278 
279  acosh=dlog(x+dsqrt(x**2-1.d0))
280 
281  end
282 c ---------------------------------------------------------------------------
283  real*8 function sgn(x)
284  real*8 x
285 
286  if(x.lt.0)then
287  sgn=-1
288  else
289  sgn=1
290  endif
291 
292  end
293 c ---------------------------------------------------------------------------
294  real*8 function facfac(l,m)
295  implicit real*8 (a-h,o-z)
296 
297  prod=1
298 
299  n=l+m-1
300 
301  do i=1,n,2
302  prod=prod*i
303  enddo
304 
305  facfac=prod
306 
307  end