PeTar
N-body code for collisional gravitational systems
stability.hpp
Go to the documentation of this file.
1 #pragma once
2 #include "Common/binary_tree.h"
3 
5 template <class Tptcl>
6 class Stability {
7 public:
8  typedef COMM::BinaryTree<Tptcl,COMM::Binary> BinTree;
9  typedef COMM::Binary Bin;
10  PS::F64 t_crit; // period criterion to determine stability
11  PS::ReallocatableArray<BinTree*> stable_binary_tree; // to store the address of sub-branch of stable systems from a binary tree
12 
13  Stability(): t_crit(-1.0), stable_binary_tree() {}
14 
15  void clear() {
16  t_crit=-1.0;
17  stable_binary_tree.clear();
18  }
19 
21  /* Use Myllaeri et al. (2018, MNRAS, 476, 830) stability criterion to check whether the system is stable.
22  @param[in] _bin_in: inner binary
23  @param[in] _bin_out: outer binary (m1 is the inner c.m. mass)
24  @param[in] _incline: inclination angle between inner and outer orbit (radians)
25  @param[in] _t_crit: time interval for stable check
26  @param[in] _is_cm_in_first: true: _bin_out m1 is the inner c.m. mass; otherwise m2
27  \return stability factor <1 stable; >1 unstable
28  */
29  static PS::F64 stable3body(const Bin& _bin_in,
30  const Bin& _bin_out,
31  const PS::F64 _incline,
32  const PS::F64 _t_crit,
33  const bool is_cm_in_first) {
34 
35  PS::F64 fac = 1.0 - 2.0*_bin_in.ecc/3.0 * (1.0 - 0.5*_bin_in.ecc*_bin_in.ecc)
36  - 0.3*std::cos(_incline)*(1.0 - 0.5*_bin_in.ecc + 2.0*std::cos(_incline)*(1.0 - 2.5*std::pow(_bin_in.ecc,1.5) - std::cos(_incline)));
37 
38  PS::F64 mout, min;
39  if (is_cm_in_first) {
40  min = _bin_out.m1;
41  mout = _bin_out.m2;
42  }
43  else {
44  min = _bin_out.m2;
45  mout = _bin_out.m1;
46  }
47  PS::F64 g = std::sqrt(std::max(_bin_in.m1,_bin_in.m2) /min)*(1.0 + mout/min);
48 
49  //Adopt at least 10,000 outer orbits for random walk time-scale.
50  PS::F64 q = 1.52*std::pow(std::sqrt(_t_crit/_bin_out.period)/(1.0 - _bin_out.ecc),1.0/6.0)*std::pow(fac*g,1.0/3.0);
51 
52  PS::F64 peri_out = _bin_out.semi * (_bin_out.ecc + 1.0);
53  PS::F64 rp = peri_out/_bin_in.semi;
54 
55  PS::F64 stab = q/rp;
56  return stab;
57  }
58 
59 
61  /* return case
62  Unstable: false hyperbolic orbit
63  stable: true apo-center < _r_crit && period < _t_crit
64  false others
65 
66  @param[in,out] _bin: binary information
67  @param[in] _r_crit: distance criterion to select stable binary
68  @param[in] _t_crit: period criterion to select stable binary
69  */
70  static bool stable2check(const Bin& _bin, const PS::F64 _r_crit, const PS::F64 _t_crit) {
71 
72 #ifdef STABLE_CHECK_DEBUG_PRINT
73  std::cerr<<"STAB2: semi="<<_bin.semi
74  <<" ecc="<<_bin.ecc
75  <<" m1="<<_bin.m1
76  <<" m2="<<_bin.m2
77  <<" period="<<_bin.period
78  <<" apo="<<_bin.semi*(1.0+_bin.ecc)
79  <<" pec="<<_bin.semi*(1.0-_bin.ecc)
80  <<std::endl;
81 #endif
82  // hyperbolic case
83  PS::F64 semi = _bin.semi;
84  if(semi<0) {
85 #ifdef STABLE_CHECK_DEBUG_PRINT
86  std::cerr<<"STAB2 reject: Hyperbolic"<<std::endl;
87 #endif
88  return false;
89  }
90 
91  // binary case
92  PS::F64 apo = _bin.semi*(1.0+_bin.ecc);
93 
94  if (apo<_r_crit&&_bin.period<_t_crit) {
95 #ifdef STABLE_CHECK_DEBUG_PRINT
96  std::cerr<<"STAB2 accept: Stable"<<std::endl;
97 #endif
98  return true;
99  }
100  else {
101 #ifdef STABLE_CHECK_DEBUG_PRINT
102  std::cerr<<"STAB2 reject: Stable but too large orbit, apo: "<<apo<<" ecc: "<<_bin.ecc<<" r_crit: "<<_r_crit
103  <<" period: "<<_bin.period<<" t_crit: "<<_t_crit<<std::endl;
104 #endif
105  return false;
106  }
107  }
108 
110  /* return case
111  false hyperbolic outer orbit
112  false apo-center outer > r_crit || period outer > t_crit
113  false stab3 >1
114  true stab3 <1
115  @param[in] _bin: inner orbit parameter
116  @param[in] _bout: outer orbit parameter
117  @param[in] _r_crit: distance criterion to select stable binary
118  @param[in] _t_crit: period criterion to select stable binary
119  @param[in] _is_cm_in_first: true: _bin_out m1 is the inner c.m. mass; otherwise m2
120  */
121  static bool stable3check(const Bin& _bin,
122  const Bin& _bout,
123  const PS::F64 _r_crit,
124  const PS::F64 _t_crit,
125  const PS::F64 _is_cm_in_first) {
126 #ifdef STABLE_CHECK_DEBUG_PRINT
127  std::cerr<<"STAB3 bout semi="<<_bout.semi
128  <<" ecc="<<_bout.ecc
129  <<" m1="<<_bout.m1
130  <<" m2="<<_bout.m2
131  <<" period="<<_bout.period
132  <<" apo="<<_bout.semi*(1.0+_bout.ecc)
133  <<" pec="<<_bout.semi*(1.0-_bout.ecc)
134  <<" bin semi="<<_bin.semi
135  <<" ecc="<<_bin.ecc
136  <<" m1="<<_bin.m1
137  <<" m2="<<_bin.m2
138  <<" period="<<_bin.period
139  <<" apo="<<_bin.semi*(1.0+_bin.ecc)
140  <<" pec="<<_bin.semi*(1.0-_bin.ecc)
141  <<std::endl;
142 #endif
143  // hyperbolic okuter orbit
144  if(_bout.semi<0) {
145 #ifdef STABLE_CHECK_DEBUG_PRINT
146  std::cerr<<"STAB3 reject: Unstable, Outer body hyperbolic, semi_out: "<<_bout.semi<<std::endl;
147 #endif
148  return false;
149  }
150  PS::F64 apo_out=_bout.semi*(1.0+_bout.ecc);
151 
152  // too large orbit
153  if(apo_out>_r_crit) {
154 #ifdef STABLE_CHECK_DEBUG_PRINT
155  std::cerr<<"STAB3 reject: Too large outer orbit, apo_out: "<<apo_out<<" r_crit: "<<_r_crit<<std::endl;
156 #endif
157  return false;
158  }
159 
160  // too large period
161  if(_bout.period>_t_crit) {
162 #ifdef STABLE_CHECK_DEBUG_PRINT
163  std::cerr<<"STAB3 reject: Too large outer period, period_out: "<<_bout.period<<" t_crit: "<<_t_crit<<std::endl;
164 #endif
165  return false;
166  }
167 
168  // stability check
169  // inclination between inner and outer orbit
170  PS::F64 incline=std::acos(std::min(1.0, _bout.am*_bin.am/std::sqrt((_bout.am*_bout.am)*(_bin.am*_bin.am))));
171  PS::F64 stab3 = stable3body(_bin, _bout, incline, _t_crit, _is_cm_in_first);
172  if(stab3>1) {
173  // Unstable case
174 #ifdef STABLE_CHECK_DEBUG_PRINT
175  std::cerr<<"STAB3 reject: Unstable, stab3: "<<stab3<<std::endl;
176 #endif
177  return false;
178  }
179  else {
180  // stable case
181 #ifdef STABLE_CHECK_DEBUG_PRINT
182  std::cerr<<"STAB3 accept: Stable, stab3: "<<stab3<<" period_in: "<<_bin.period<<" period_out: "<<_bout.period<<std::endl;
183 #endif
184  return true;
185  }
186  }
187 
189  /* return tstep stable_factor case
190  Unstable: false -1 -1 hyperbolic outer orbit
191  false -1 -1 apo-center outer > _r_crit
192  false -1 -1 period outer > 0.25 * dt_tree
193  true inner -stab3_max stab3_1 >0.8 || stab3_2 > 0.8 & apo_out <= r_out
194  false -1 -1 & apo_out > r_out
195  stable: true inner -1 stab3_1 <=0.8 & stab3_2 <=0.8 & unpert & outer period > 1/8 dt_tree
196  false -1 -1 -- & acceleration ratio(out/in) <1e-6 and outer period >1e-4 dt_tree
197  true inner fpert/(m_out/apo_out^2) -- & other cases
198  @param[in] _bout: outer orbit parameter
199  @param[in] _bin1: first inner orbit parameter
200  @param[in] _bin2: second inner orbit parameter
201  @param[in] _r_crit: distance criterion to select stable binary
202  @param[in] _t_crit: period criterion to select stable binary
203  */
204  static bool stable4check(const Bin& _bin1,
205  const Bin& _bin2,
206  const Bin& _bout,
207  const PS::F64 _r_crit,
208  const PS::F64 _t_crit) {
209 #ifdef STABLE_CHECK_DEBUG_PRINT
210  std::cerr<<"STAB4 bout semi="<<_bout.semi
211  <<" ecc="<<_bout.ecc
212  <<" m1="<<_bout.m1
213  <<" m2="<<_bout.m2
214  <<" period="<<_bout.period
215  <<" apo="<<_bout.semi*(1.0+_bout.ecc)
216  <<" pec="<<_bout.semi*(1.0-_bout.ecc)
217  <<" bin1 semi="<<_bin1.semi
218  <<" ecc="<<_bin1.ecc
219  <<" m1="<<_bin1.m1
220  <<" m2="<<_bin1.m2
221  <<" period="<<_bin1.period
222  <<" apo="<<_bin1.semi*(1.0+_bin1.ecc)
223  <<" pec="<<_bin1.semi*(1.0-_bin1.ecc)
224  <<" bin2 semi="<<_bin2.semi
225  <<" ecc="<<_bin2.ecc
226  <<" m1="<<_bin2.m1
227  <<" m2="<<_bin2.m2
228  <<" period="<<_bin2.period
229  <<" apo="<<_bin2.semi*(1.0+_bin2.ecc)
230  <<" pec="<<_bin2.semi*(1.0-_bin2.ecc)
231  <<std::endl;
232 #endif
233  // hyperbolic outer orbit
234  if(_bout.semi<0) {
235 #ifdef STABLE_CHECK_DEBUG_PRINT
236  std::cerr<<"STAB4 reject: Unstable, Outer body hyperbolic, semi_out: "<<_bout.semi<<std::endl;
237 #endif
238  return false;
239  }
240 
241  PS::F64 apo_out=_bout.semi*(1.0+_bout.ecc);
242  //PS::F64 pec_out=_bout.semi*(1.0-_bout.ecc);
243  //PS::F64 apo_in1=_bin1.semi*(1.0+_bin1.ecc);
244  //PS::F64 apo_in2=_bin2.semi*(1.0+_bin2.ecc);
245 
246 
247  // too large orbit
248  if(apo_out>_r_crit) {
249 #ifdef STABLE_CHECK_DEBUG_PRINT
250  std::cerr<<"STAB4 reject: Unstable, Too large outer orbit, apo_out: "<<apo_out<<" r_crit: "<<_r_crit<<std::endl;
251 #endif
252  return false;
253  }
254 
255  // too large period
256  if(_bout.period>_t_crit) {
257 #ifdef STABLE_CHECK_DEBUG_PRINT
258  std::cerr<<"STAB4 reject: Too large outer period, period_out: "<<_bout.period<<" t_crit: "<<_t_crit<<std::endl;
259 #endif
260  return false;
261  }
262 
263  // stability check
264  // inclination between inner and outer orbit
265  PS::F64 incl1=std::acos(std::min(1.0, _bout.am*_bin1.am/std::sqrt((_bout.am*_bout.am)*(_bin1.am*_bin1.am))));
266  PS::F64 incl2=std::acos(std::min(1.0, _bout.am*_bin2.am/std::sqrt((_bout.am*_bout.am)*(_bin2.am*_bin2.am))));
267  PS::F64 stab3_1 = stable3body(_bin1, _bout, incl1, _t_crit, true);
268  PS::F64 stab3_2 = stable3body(_bin2, _bout, incl2, _t_crit, false);
269  if(stab3_1>0.8||stab3_2>0.8) {
270  // Unstable case
271 #ifdef STABLE_CHECK_DEBUG_PRINT
272  std::cerr<<"STAB4 reject: Unstable, stab3_1: "<<stab3_1<<" stab3_2: "<<stab3_2<<std::endl;
273 #endif
274  return false;
275  }
276  else {
277  // stable case
278 #ifdef STABLE_CHECK_DEBUG_PRINT
279  std::cerr<<"STAB4 accept: Stable, stab3_1: "<<stab3_1<<" stab3_2: "<<stab3_2<<std::endl;
280 #endif
281  return true;
282  }
283  }
284 
286 
291  static PS::S32 closedSystemCheckIter(Stability& _stab, const PS::S32& _stab_res1, const PS::S32& _stab_res2, BinTree& _bin) {
292  // unstable case
293  if (_stab_res1==0||_stab_res2==0) {
294  // if any branch is stable, save data
295  if (_stab_res1==1) _stab.stable_binary_tree.push_back((BinTree*)_bin.getLeftMember());
296  if (_stab_res2==1) _stab.stable_binary_tree.push_back((BinTree*)_bin.getRightMember());
297  return 0;
298  }
299 
300  PS::F64 apo = _bin.semi*(1+_bin.ecc);
301  bool stab_flag = (apo>0.0 && apo<=_bin.getRGroup());
302  // B-B system
303  if (_stab_res1==1&&_stab_res2==1) {
304  if (stab_flag) return 1;
305  else {
306  _stab.stable_binary_tree.push_back((BinTree*)_bin.getLeftMember());
307  _stab.stable_binary_tree.push_back((BinTree*)_bin.getRightMember());
308  return 0;
309  }
310  }
311  // triple case
312  if (_stab_res1==-1&&_stab_res2==1) {
313  if (stab_flag) return 1;
314  else {
315  _stab.stable_binary_tree.push_back((BinTree*)_bin.getRightMember());
316  return 0;
317  }
318  }
319  if (_stab_res2==-1&&_stab_res1==1) {
320  if (stab_flag) return 1;
321  else {
322  _stab.stable_binary_tree.push_back((BinTree*)_bin.getLeftMember());
323  return 0;
324  }
325  }
326  // binary case
327 #ifdef STABLE_CHECK_DEBUG
328  assert(_stab_res1==-1&&_stab_res2==-1);
329 #endif
330  if (stab_flag) return 1;
331  else return 0;
332  }
333 
335 
339  void findClosedTree(BinTree& _bin) {
340  stable_binary_tree.resizeNoInitialize(0);
341  bool is_stable= _bin.processTreeIter(*this, -1, -1, closedSystemCheckIter);
342  if (is_stable) stable_binary_tree.push_back(&_bin);
343  }
344 
346 
351  static PS::S32 stabilityCheckIter(Stability& _stab, const PS::S32& _stab_res1, const PS::S32& _stab_res2, BinTree& _bin) {
352  // B-B system
353  if (_stab_res1==1&&_stab_res2==1) {
354  bool stab4=stable4check(*(BinTree*)_bin.getLeftMember(), *(BinTree*)_bin.getRightMember(), _bin, _bin.getRGroup(), _stab.t_crit);
355  if (stab4) return 1;
356  else {
357  _stab.stable_binary_tree.push_back((BinTree*)_bin.getLeftMember());
358  _stab.stable_binary_tree.push_back((BinTree*)_bin.getRightMember());
359  return 0;
360  }
361  }
362  // unstable case
363  if (_stab_res1==0||_stab_res2==0) {
364  // if any branch is stable, save data
365  if (_stab_res1==1) _stab.stable_binary_tree.push_back((BinTree*)_bin.getLeftMember());
366  if (_stab_res2==1) _stab.stable_binary_tree.push_back((BinTree*)_bin.getRightMember());
367  return 0;
368  }
369  // triple case
370  if (_stab_res1==-1&&_stab_res2==1) {
371  bool stab3=stable3check(*(BinTree*)_bin.getRightMember(), _bin, _bin.getRGroup(), _stab.t_crit, false);
372  if (stab3) return 1;
373  else {
374  _stab.stable_binary_tree.push_back((BinTree*)_bin.getRightMember());
375  return 0;
376  }
377  }
378  if (_stab_res2==-1&&_stab_res1==1) {
379  bool stab3=stable3check(*(BinTree*)_bin.getLeftMember(), _bin, _bin.getRGroup(), _stab.t_crit, true);
380  if (stab3) return 1;
381  else {
382  _stab.stable_binary_tree.push_back((BinTree*)_bin.getLeftMember());
383  return 0;
384  }
385  }
386  // binary case
387 #ifdef STABLE_CHECK_DEBUG
388  assert(_stab_res1==-1&&_stab_res2==-1);
389 #endif
390  bool stab2 = stable2check(_bin, _bin.getRGroup(), _stab.t_crit);
391  if (stab2) return 1;
392  else return 0;
393  }
394 
396 
399  void findStableTree(BinTree& _bin) {
400  stable_binary_tree.resizeNoInitialize(0);
401  bool is_stable= _bin.processTreeIter(*this, -1, -1, stabilityCheckIter);
402  if (is_stable) stable_binary_tree.push_back(&_bin);
403  }
404 
405 };
Stability
stability checker for binary tree
Definition: stability.hpp:6
PIKG::S32
int32_t S32
Definition: pikg_vector.hpp:24
Stability::findClosedTree
void findClosedTree(BinTree &_bin)
Check binary tree and collect closed systems subtree root.
Definition: stability.hpp:339
Stability::closedSystemCheckIter
static PS::S32 closedSystemCheckIter(Stability &_stab, const PS::S32 &_stab_res1, const PS::S32 &_stab_res2, BinTree &_bin)
Closed system check iteration function for BinaryTree.
Definition: stability.hpp:291
Stability::Bin
COMM::Binary Bin
Definition: stability.hpp:9
Stability::stable4check
static bool stable4check(const Bin &_bin1, const Bin &_bin2, const Bin &_bout, const PS::F64 _r_crit, const PS::F64 _t_crit)
Four-body (B-B) stability check.
Definition: stability.hpp:204
Stability::stable2check
static bool stable2check(const Bin &_bin, const PS::F64 _r_crit, const PS::F64 _t_crit)
check two-body stability
Definition: stability.hpp:70
PIKG::F64
double F64
Definition: pikg_vector.hpp:17
Stability::Stability
Stability()
Definition: stability.hpp:13
PIKG::min
T min(const Vector3< T > &v)
Definition: pikg_vector.hpp:150
Stability::t_crit
PS::F64 t_crit
Definition: stability.hpp:10
Stability::stable3body
static PS::F64 stable3body(const Bin &_bin_in, const Bin &_bin_out, const PS::F64 _incline, const PS::F64 _t_crit, const bool is_cm_in_first)
Three-body stability function.
Definition: stability.hpp:29
Stability::stabilityCheckIter
static PS::S32 stabilityCheckIter(Stability &_stab, const PS::S32 &_stab_res1, const PS::S32 &_stab_res2, BinTree &_bin)
Stability check iteraction function for BinaryTree.
Definition: stability.hpp:351
Stability::clear
void clear()
Definition: stability.hpp:15
Stability::stable3check
static bool stable3check(const Bin &_bin, const Bin &_bout, const PS::F64 _r_crit, const PS::F64 _t_crit, const PS::F64 _is_cm_in_first)
Three-body (B-S) stability check.
Definition: stability.hpp:121
PIKG::max
T max(const Vector3< T > &v)
Definition: pikg_vector.hpp:143
Stability::stable_binary_tree
PS::ReallocatableArray< BinTree * > stable_binary_tree
Definition: stability.hpp:11
Stability::BinTree
COMM::BinaryTree< Tptcl, COMM::Binary > BinTree
Definition: stability.hpp:8
Stability::findStableTree
void findStableTree(BinTree &_bin)
Check binary tree and collect stable subtree.
Definition: stability.hpp:399