SlowDown Algorithmic Regularization (SDAR)
Algorithmic Regularization with slowdown method for integrating few-body motions
binary_tree.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include<cmath>
4 #include<algorithm>
5 #include"Common/Float.h"
6 #include"Common/matrix.h"
7 
8 namespace COMM{
9  const Float PI = 4.0*atan(1.0);
10 
12  class Binary{
13  public:
14  Float semi; // semi: semi-major axis
15  Float ecc; // ecc: eccentricity
16  Float incline; // incline: inclination angle, Eular angle beta
17  Float rot_horizon; // rotational angle in horizon, Eular angle alpha
18  Float rot_self; // rotational angle in orbital plane, Eular angle gamma
19  Float t_peri; // t_peri: time to peri-center
20  Float period; // period: period
21  Float ecca; // ecca: eccentricty anomaly (-pi, pi)
22  Float m1; // m1: mass 1
23  Float m2; // m2: mass 2
24  Float r; // distance between too members
25  Vector3<Float> am; // specific angular momentum
26  Float stab; // stability factor
27 
29 
36  template <class Tptcl>
37  static void orbitToParticle(Tptcl& _p1, Tptcl& _p2, const Binary& _bin, const Float& _ecca, const Float _G) {
38  Float m_tot = _bin.m1 + _bin.m2;
39  Vector3<Float> pos_star, vel_star;
40  // hyper bolic orbit
41  if (_bin.semi<0) {
42  Float coshu = cosh(_ecca);
43  Float sinhu = sinh(_ecca);
44  Float eccfac = sqrt(_bin.ecc*_bin.ecc - 1.0);
45  pos_star.x = -_bin.semi*(_bin.ecc - coshu); // -a (e - cosh(E))
46  pos_star.y = -_bin.semi*eccfac*sinhu; // -a sqrt(e^2-1) sinh(E)
47  pos_star.z = 0.0;
48  Float afac = sqrt(_G*m_tot / -_bin.semi); // mean_motion * a
49  Float fac = afac/(_bin.ecc*coshu-1.0); // a/r * sqrt(GM/a)
50  vel_star.x = -fac*sinhu; // - sqrt(GM/-a) sinh(E)/(e cosh(E)-1)
51  vel_star.y = fac*eccfac*coshu; // sqrt(GM(e^2-1)/-a) cosh(E)/(e cosh(E)-1)
52  vel_star.z = 0.0;
53  }
54  else{ // ellipse orbit
55  Float cosu = cos(_ecca);
56  Float sinu = sin(_ecca);
57  Float eccfac = sqrt(1.0 - _bin.ecc*_bin.ecc);
58  pos_star.x = _bin.semi*(cosu - _bin.ecc); // a (cos(E) - e)
59  pos_star.y = _bin.semi*eccfac*sinu; // a sqrt(1-e^2) sin(E)
60  pos_star.z = 0.0;
61  Float afac = sqrt(_G*m_tot / _bin.semi);
62  Float fac = afac/(1.0-_bin.ecc*cosu); // a/r * sqrt(GM/a)
63  vel_star.x = -fac*sinu; // - sqrt(GM/a) sin(E)/(1 - e cos(E))
64  vel_star.y = fac*eccfac*cosu; // sqrt(GM(1-e^2)/a) cos(E)/(1 - e cos(E))
65  vel_star.z = 0.0;
66  }
67  Matrix3<Float> rot;
68  rot.rotation(_bin.incline, _bin.rot_horizon, _bin.rot_self);
69  Vector3<Float> pos_red = rot*pos_star;
70  Vector3<Float> vel_red = rot*vel_star;
71  _p1.mass = _bin.m1;
72  _p2.mass = _bin.m2;
73 
74  Float m2_mt = _p2.mass / m_tot;
75  Float m1_mt = _p1.mass / m_tot;
76  _p1.pos[0] = -m2_mt * pos_red.x;
77  _p1.pos[1] = -m2_mt * pos_red.y;
78  _p1.pos[2] = -m2_mt * pos_red.z;
79 
80  _p2.pos[0] = m1_mt * pos_red.x;
81  _p2.pos[1] = m1_mt * pos_red.y;
82  _p2.pos[2] = m1_mt * pos_red.z;
83 
84  _p1.vel[0] = -m2_mt * vel_red.x;
85  _p1.vel[1] = -m2_mt * vel_red.y;
86  _p1.vel[2] = -m2_mt * vel_red.z;
87 
88  _p2.vel[0] = m1_mt * vel_red.x;
89  _p2.vel[1] = m1_mt * vel_red.y;
90  _p2.vel[2] = m1_mt * vel_red.z;
91  }
92 
94  /* refer to the P3T code developed by Iwasawa M.
95  @param[out] _bin: binary parameter
96  @param[in] _p1: particle 1
97  @param[in] _p2: particle 2
98  @param[in] _G: gravitational constant
99  */
100  template <class Tptcl>
101  static void particleToOrbit(Binary& _bin, const Tptcl& _p1, const Tptcl& _p2, const Float _G){
102  _bin.m1 = _p1.mass;
103  _bin.m2 = _p2.mass;
104  Float m_tot = _p1.mass + _p2.mass;
105  Float Gm_tot = _G*m_tot;
106  Vector3<Float> pos_red(_p2.pos[0] - _p1.pos[0], _p2.pos[1] - _p1.pos[1], _p2.pos[2] - _p1.pos[2]);
107  Vector3<Float> vel_red(_p2.vel[0] - _p1.vel[0], _p2.vel[1] - _p1.vel[1], _p2.vel[2] - _p1.vel[2]);
108  Float r_sq = pos_red * pos_red;
109  _bin.r = sqrt(r_sq);
110  Float inv_dr = 1.0 / _bin.r;
111  Float v_sq = vel_red * vel_red;
112  _bin.semi = 1.0 / (2.0*inv_dr - v_sq / Gm_tot);
113  _bin.am = pos_red ^ vel_red;
114  _bin.incline = atan2( sqrt(_bin.am.x*_bin.am.x+_bin.am.y*_bin.am.y), _bin.am.z);
115  _bin.rot_horizon = _bin.am.y!=0.0? atan2(_bin.am.x, -_bin.am.y) : 0.0;
116 
117  // position and velocity in rest frame
118  Vector3<Float> pos_bar, vel_bar;
119  // OMG: rot_horizon; omg: rot_self; inc: incline
120  Float cosOMG = cos(_bin.rot_horizon);
121  Float sinOMG = sin(_bin.rot_horizon);
122  Float cosinc = cos(_bin.incline);
123  Float sininc = sin(_bin.incline);
124  pos_bar.x = pos_red.x*cosOMG + pos_red.y*sinOMG;
125  pos_bar.y = (-pos_red.x*sinOMG + pos_red.y*cosOMG)*cosinc + pos_red.z*sininc;
126  pos_bar.z = 0.0;
127  vel_bar.x = vel_red.x*cosOMG + vel_red.y*sinOMG;
128  vel_bar.y = (-vel_red.x*sinOMG + vel_red.y*cosOMG)*cosinc + vel_red.z*sininc;
129  vel_bar.z = 0.0;
130  // h: specific angular momentum value (|r x v|)
131  Float h = sqrt(_bin.am*_bin.am);
132  // eccentricity vector: (v x h)/GM - r/|r|
133  // angle between AM and vel is 90 degree
134  Float ecccosomg = h/Gm_tot*vel_bar.y - pos_bar.x*inv_dr;
135  Float eccsinomg = -h/Gm_tot*vel_bar.x - pos_bar.y*inv_dr;
136  _bin.ecc = sqrt( ecccosomg*ecccosomg + eccsinomg*eccsinomg );
137  _bin.rot_self = atan2(eccsinomg, ecccosomg);
138  // angle of position referring to the rest frame: f + omg (f: true anomaly)
139  Float phi = atan2(pos_bar.y, pos_bar.x);
140  Float true_anomaly = phi - _bin.rot_self;
141  // eccentric anomaly
142  //Float sin_ecca = _bin.r*sin(true_anomaly) / (_bin.semi*sqrt(1.0 - _bin.ecc*_bin.ecc));
143  //Float cos_ecca = (_bin.r*cos(true_anomaly) / _bin.semi) + _bin.ecc;
144  if (_bin.semi>0) {
145  _bin.ecca = atan2(sin(true_anomaly)*sqrt(1.0 - _bin.ecc*_bin.ecc), _bin.ecc+ cos(true_anomaly));
146  Float mean_motion = sqrt(Gm_tot/(_bin.semi*_bin.semi*_bin.semi));
147  _bin.period = 2*PI/mean_motion;
148  Float mean_anomaly = _bin.ecca - _bin.ecc*sin(_bin.ecca);
149  _bin.t_peri = mean_anomaly / mean_motion;
150  }
151  else {
152  _bin.ecca = atanh(sin(true_anomaly)*sqrt(_bin.ecc*_bin.ecc-1)/(_bin.ecc+ cos(true_anomaly)));
153  Float mean_motion = sqrt(Gm_tot/(-_bin.semi*_bin.semi*_bin.semi));
154  _bin.period = 2*PI/mean_motion;
155  Float mean_anomaly = _bin.ecc*sinh(_bin.ecca) - _bin.ecca;
156  _bin.t_peri = mean_anomaly / mean_motion;
157  }
158  }
159 
161 
166  static Float periodToSemi(const Float& _period, const Float& _mtot, const Float& _G) {
167  return pow(_period*_period*_G*_mtot/(4*PI*PI),1.0/3.0);
168  }
169 
171 
176  static Float semiToPeriod(const Float& _semi, const Float& _mtot, const Float& _G) {
177  Float mean_motion = sqrt(_G*_mtot/(fabs(_semi*_semi*_semi)));
178  return 2*PI/mean_motion;
179  }
180 
181 
183 
191  template <class Tpi, class Tpj>
192  static void particleToSemiEcc(Float& _semi, Float& _ecc, Float& _r, Float& _rv, const Tpi& _p1, const Tpj& _p2, const Float _G){
193  Float m_tot = _p1.mass + _p2.mass;
194  Float Gm_tot = _G*m_tot;
195  Vector3<Float> pos_red(_p2.pos[0] - _p1.pos[0], _p2.pos[1] - _p1.pos[1], _p2.pos[2] - _p1.pos[2]);
196  Vector3<Float> vel_red(_p2.vel[0] - _p1.vel[0], _p2.vel[1] - _p1.vel[1], _p2.vel[2] - _p1.vel[2]);
197  Float r_sq = pos_red * pos_red;
198  _r = sqrt(r_sq);
199  Float v_sq = vel_red * vel_red;
200  _rv = pos_red * vel_red;
201  _semi = 1.0 / (2.0 / _r - v_sq / Gm_tot);
202  Float p = 1.0 - _r/_semi;
203  _ecc = sqrt(p*p + _rv*_rv/_semi/Gm_tot);
204  }
205 
207 
212  template <class Tpi, class Tpj>
213  void particleToSemiEccPeriod(const Tpi& _p1, const Tpj& _p2, const Float _G) {
214  Float r,rv;
215  particleToSemiEcc(semi, ecc, r, rv, _p1, _p2, _G);
216  Float mtot=m1+m2;
217  period = semiToPeriod(semi, mtot, _G);
218  }
219 
221 
225  Float calcEccAnomaly(const Float _r) {
226  if (semi>0) {
227  Float cos_ecca = (1.0-_r/semi)/ecc;
228  return acos(cos_ecca);
229  }
230  else {
231  Float cosh_ecca = (1.0-_r/semi)/ecc;
232  return acosh(cosh_ecca);
233  }
234  }
235 
237 
242  static Float calcMeanAnomaly(const Float _ecca, const Float _ecc) {
243  if (_ecc<1.0) return _ecca - _ecc*sin(_ecca);
244  else return _ecc*sinh(_ecca) - _ecca;
245  }
246 
248  /* refer to the P3T code developed by Iwasawa M.
249  @param[in] _mean_anomaly: mean_anomaly
250  @param[in] _ecc: eccentricity
251  \return eccentric anomaly
252  */
253  static Float calcEccAnomaly(const Float _mean_anomaly,
254  const Float _ecc){
255  // a: semi-major axis
256  // l: mean anomaly
257  // e: eccentricity
258  // u: eccentric anomaly
259  // n: mean mortion
260  Float u0 = _mean_anomaly;
261  Float u1;
262  int loop = 0;
263  while(1){
264  loop++;
265  Float su0 = sin(u0);
266  Float cu0 = sqrt(1.0 - su0*su0);
267  //u1 = u0 - keplereq(l, e, u0)/keplereq_dot(e, u0);
268  u1 = u0 - ((u0- _ecc*su0-_mean_anomaly)/(1.0 - _ecc*cu0));
269  //if( fabs(u1-u0) < 1e-13 ){ return u1; }
270  if( fabs(u1-u0) < 1e-15 ){ return u1; }
271  else{ u0 = u1; }
272  if (loop>1e5) {
273  std::cerr<<"Error: kepler solver cannot converge to find correct eccentricity anomaly!\n";
274  abort();
275  }
276  }
277  }
278 
280 
284  static void solveKepler(Binary& _bin,
285  const Float _dt) {
286  Float freq = sqrt( (_bin.m1+_bin.m2) / (_bin.semi*_bin.semi*_bin.semi) );
287  Float mean_anomaly_old = _bin.ecca - _bin.ecc * sin(_bin.ecca);
288  Float dt_tmp = _dt - to_int(_dt/_bin.period)*_bin.period;
289  Float mean_anomaly_new = freq * dt_tmp + mean_anomaly_old; // mean anomaly
290  _bin.ecca = calcEccAnomaly(mean_anomaly_new, _bin.ecc); // eccentric anomaly
291  }
292 
294 
299  template <class Tptcl>
300  void calcOrbit(const Tptcl& _p1, const Tptcl& _p2, const Float _G) {
301  particleToOrbit(*this, _p1, _p2, _G);
302  }
303 
305 
310  template <class Tptcl>
311  void calcParticles(Tptcl& _p1, Tptcl& _p2, const Float& _G) {
312  orbitToParticle(_p1, _p2, *this, this->ecca, _G);
313  }
314 
316 
319  void calcSemiFromPeriod(const Float& _G) {
320  Float mtot = m1+m2;
321  semi = periodToSemi(period, mtot, _G);
322  }
323 
325 
331  template <class Tptcl>
332  void calcParticlesEcca(Tptcl& _p1, Tptcl& _p2, const Float _ecca, const Float _G) const {
333  orbitToParticle(_p1, _p2, *this, _ecca, _G);
334  }
335 
337 
340  void rotateToOriginalFrame(Float* _vec ) {
341  Matrix3<Float> rot;
342  rot.rotation(incline, rot_horizon, rot_self);
343  Vector3<Float> vec (_vec[0], _vec[1], _vec[2]);
344  vec = rot*vec;
345  _vec[0] = vec.x;
346  _vec[1] = vec.y;
347  _vec[2] = vec.z;
348  }
349 
351 
354  void evolve(const Float _dt) {
355  solveKepler(*this, _dt);
356  }
357 
358  // IO functions
359  void print(std::ostream& _os) const{
360  _os<<"semi: semi-major axis: "<<semi<<std::endl
361  <<"ecc: eccentricity: "<<ecc<<std::endl
362  <<"incline: inclination angle: "<<incline<<std::endl
363  <<"rot_horizon: rotational angle in horizon: "<<rot_horizon<<std::endl
364  <<"rot_self: rotational angle in orbital plane: "<<rot_self<<std::endl
365  <<"t_peri: phase to peri-center: "<<t_peri<<std::endl
366  <<"period: period: "<<period<<std::endl
367  <<"ecca: eccentricty anomaly: "<<ecca<<std::endl
368  <<"m1: mass 1: "<<m1<<std::endl
369  <<"m2: mass 2: "<<m2<<std::endl
370  <<"r: separation: "<<r<<std::endl
371  <<"L: Angular momentum: "<<am.x<<" "<<am.y<<" "<<am.z<<std::endl
372  <<"stab: stability factor: "<<stab<<std::endl;
373  }
374 
376 
380  static void printColumnTitle(std::ostream & _fout, const int _width=20) {
381  _fout<<std::setw(_width)<<"semi"
382  <<std::setw(_width)<<"ecc"
383  <<std::setw(_width)<<"incline"
384  <<std::setw(_width)<<"rot_horizon"
385  <<std::setw(_width)<<"rot_self"
386  <<std::setw(_width)<<"t_peri"
387  <<std::setw(_width)<<"period"
388  <<std::setw(_width)<<"ecca"
389  <<std::setw(_width)<<"m1"
390  <<std::setw(_width)<<"m2"
391  <<std::setw(_width)<<"r"
392  <<std::setw(_width)<<"L.x"
393  <<std::setw(_width)<<"L.y"
394  <<std::setw(_width)<<"L.z"
395  <<std::setw(_width)<<"stab";
396  }
397 
399 
403  void printColumn(std::ostream & _fout, const int _width=20){
404  _fout<<std::setw(_width)<<semi
405  <<std::setw(_width)<<ecc
406  <<std::setw(_width)<<incline
407  <<std::setw(_width)<<rot_horizon
408  <<std::setw(_width)<<rot_self
409  <<std::setw(_width)<<t_peri
410  <<std::setw(_width)<<period
411  <<std::setw(_width)<<ecca
412  <<std::setw(_width)<<m1
413  <<std::setw(_width)<<m2
414  <<std::setw(_width)<<r
415  <<std::setw(_width)<<am.x
416  <<std::setw(_width)<<am.y
417  <<std::setw(_width)<<am.z
418  <<std::setw(_width)<<stab;
419  }
420 
421 
423 
425  void writeBinary(FILE *_fp) const {
426  fwrite(this, sizeof(*this),1,_fp);
427  }
428 
430 
432  void readBinary(FILE *_fin) {
433  size_t rcount = fread(this, sizeof(*this),1,_fin);
434  if (rcount<1) {
435  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
436  abort();
437  }
438  }
439 
441 
443  void writeAscii(std::ostream& _fout) const {
444  _fout<<semi<<" "
445  <<ecc<<" "
446  <<incline<<" "
447  <<rot_horizon<<" "
448  <<rot_self<<" "
449  <<t_peri<<" "
450  <<period<<" "
451  <<ecca<<" "
452  <<m1<<" "
453  <<m2<<" "
454  <<r<<" "
455  <<am.x<<" "
456  <<am.y<<" "
457  <<am.z<<" "
458  <<stab<<" ";
459  }
460 
462 
464  void readAscii(std::istream& _fin) {
465  _fin>>semi>>ecc>>incline>>rot_horizon>>rot_self>>t_peri>>period>>ecca>>m1>>m2>>r>>am.x>>am.y>>am.z>>stab;
466  }
467 
468  //void PosVel2SemiEcc(Float & semi,
469  // Float & ecc,
470  // const Vector3<Float> & pos0, const Vector3<Float> & pos1,
471  // const Vector3<Float> & vel0, const Vector3<Float> & vel1,
472  // const Float mass0, const Float mass1) {
473  // Float m_tot = mass0 + mass1;
474  // Vector3<Float> pos_red = pos1 - pos0;
475  // Vector3<Float> vel_red = vel1 - vel0;
476  // Float r_sq = pos_red * pos_red;
477  // Float r = sqrt(r_sq);
478  // Float inv_dr = 1.0 / r;
479  // Float v_sq = vel_red * vel_red;
480  // semi = 1.0 / (2.0*inv_dr - v_sq / m_tot);
481  //
482  // Float rs = 1.0 - r/semi;
483  // Float rv = pos_red * vel_red;
484  // ecc = std::sqrt(rs*rs + rv*rv/(m_tot*semi));
485  //}
486  };
487 
489 
491  template <class Tptcl, class Tbinary>
492  class BinaryTree: public Tptcl, public Tbinary {
493  private:
494  typedef std::pair<Float, int> R2Index; // First is distance square, second is index
495  typedef BinaryTree<Tptcl,Tbinary> BinaryTreeLocal; // local defined binarytree
496  /* Example for level & branch
497  -------------------------------------------------
498  level branch
499  0 0
500  / \
501  1 0 1
502  / \ / \
503  2 0 1 2 3
504  -------------------------------------------------*/
505  int n_members;
506  int member_index[2];
507  Tptcl* member[2];
508 
509  public:
510  int level;
511  int branch;
512 
513  private:
514  static bool pairLess(const R2Index& a, const R2Index & b) {
515  return a.first < b.first;
516  }
517 
519  /* @param[in,out] _member_list: group particle member index, will be reordered by the minimum distance chain.
520  @param[in,out] _r2_list: a pair. first is the square distance of closes neighbor (next particle i+1 in _member_list); second is equal to particle current index i
521  @param[in] _n_members: number of members
522  @param[in] _ptcl_org: original particle data
523  */
524  static void calcMinDisList(int _ptcl_list[],
525  R2Index _r2_list[],
526  const int _n_members,
527  const Tptcl* _ptcl) {
528  for (int i=0; i<_n_members-1; i++) {
529  int k = _ptcl_list[i];
530  int jc=-1;
531  Float r2min = NUMERIC_FLOAT_MAX;
532  for(int j=i+1; j<_n_members; j++) {
533  const int kj =_ptcl_list[j];
534  Vector3<Float> dr = {_ptcl[k].pos[0] - _ptcl[kj].pos[0],
535  _ptcl[k].pos[1] - _ptcl[kj].pos[1],
536  _ptcl[k].pos[2] - _ptcl[kj].pos[2]};
537  Float r2 = dr*dr;
538  if(r2<r2min) {
539  r2min = r2;
540  jc = j;
541  }
542  }
543 #ifdef BINARY_DEBUG
544  ASSERT(jc>=0);
545 #endif
546  if (jc!=i+1) {
547  int jtmp = _ptcl_list[i+1];
548  _ptcl_list[i+1] = _ptcl_list[jc];
549  _ptcl_list[jc] = jtmp;
550  }
551  _r2_list[i].first = r2min;
552  _r2_list[i].second = i;
553  }
554  }
555 
556 
557  public:
558  BinaryTree(): Tptcl(), Tbinary(), n_members(-1), member_index{-1,-1}, member{NULL,NULL}, level(-1), branch(-1) {}
559 
560 
561  // copy function
562  //BinaryTreeLocal & operator = (const BinaryTreeLocal& _bin) {
563  // std::memcpy(this, &_bin, sizeof(Tbinary)+sizeof(Tptcl)+sizeof(int));
564  // return *this;
565  //}
566 
568  /* Use Myllaeri et al. (2018, MNRAS, 476, 830) 3-body stability criterion to check whether any B-S/B-B sub-system is stable, return the maximim stabliity factor
569  @param[in] _bin: binary tree to check
570  @param[in] _t_crit: time interval for stable check
571  \return maximum stability factor <1 stable; >1 unstable; =0 binary system
572  */
573  static Float& stableCheckIter(BinaryTreeLocal& _bin, const Float _t_crit) {
574  if (_bin.semi<0) _bin.stab=NUMERIC_FLOAT_MAX;
575  else _bin.stab = 0.0;
576  Float mout = _bin.mass;
577 
578  for (int k=0; k<2; k++) {
579  if (_bin.isMemberTree(k)) {
580  BinaryTreeLocal* bin_in = _bin.getMemberAsTree(k);
581  _bin.stab = std::max(_bin.stab,stableCheckIter(*bin_in, _t_crit));
582  if (_bin.stab<1.0) {
583  // if zero mass exist, set unstable
584  if (_bin.getMember(1-k)->mass ==0.0) {
585  _bin.stab = std::max(_bin.stab, Float(1000.0));
586  }
587  else {
588 
589  Float incline=acos(std::min(Float(1.0), _bin.am*bin_in->am/sqrt((_bin.am*_bin.am)*(bin_in->am*bin_in->am))));
590 
591  Float cosinc = cos(incline);
592  Float fac = 1.0 - 2.0*bin_in->ecc/3.0 * (1.0 - 0.5*bin_in->ecc*bin_in->ecc)
593  - 0.3*cosinc*(1.0 - 0.5*bin_in->ecc + 2.0*cosinc*(1.0 - 2.5*pow(bin_in->ecc,1.5) - cosinc));
594 
595  Float min = bin_in->mass;
596  Float g = sqrt(std::max(bin_in->m1,bin_in->m2) /min)*(1.0 + mout/min);
597 
598  Float qst = 1.52*pow(sqrt(_t_crit/_bin.period)/(1.0 - _bin.ecc),1.0/6.0)*pow(fac*g,1.0/3.0);
599 
600  Float peri_out = _bin.semi * (_bin.ecc + 1.0);
601  Float q = peri_out/bin_in->semi;
602 
603  bin_in->stab = std::max(bin_in->stab, qst/q);
604  _bin.stab = std::max(_bin.stab, qst/q);
605  }
606  }
607  }
608  }
609  return _bin.stab;
610  }
611 
613 
616  static void generateLevelBranch(BinaryTreeLocal& _bin, const bool _reset) {
617  if (_reset) {
618  _bin.level=0;
619  _bin.branch=0;
620  }
621  else {
622  ASSERT(_bin.level>=0&&_bin.branch>=0);
623  }
624  for (int i=0; i<2; i++) {
625  if(_bin.isMemberTree(i)) {
626  auto* bini = _bin.getMemberAsTree(i);
627  bini->level = _bin.level+1;
628  bini->branch = 2*_bin.branch + i;
629  generateLevelBranch(*bini, false);
630  }
631  }
632  };
633 
635 
638  int isSameBranch(const BinaryTreeLocal& _bin) const {
639  ASSERT(level>=0&&branch>=0&&_bin.level>=0&&_bin.branch>=0);
640  int dlevel = _bin.level - level;
641  if (dlevel>0) { // check leaf
642  int left_edge = (branch<<dlevel);
643  int right_edge = ((branch+1)<<dlevel);
644  if (_bin.branch>=left_edge&&_bin.branch<right_edge) return 2;
645  else return 0;
646  }
647  else if (dlevel<0) { // check mother
648  int left_edge = (_bin.branch<<-dlevel);
649  int right_edge = ((_bin.branch+1)<<-dlevel);
650  if (branch>=left_edge&&branch<right_edge) return -2;
651  else return 0;
652  }
653  else if(_bin.branch==branch) return 1;
654  else return 0;
655  }
656 
658  /*
659  @param[in,out] _bins: binary tree, size of n_members
660  @param[in,out] _ptcl_list: group particle member index in _ptcl, will be reordered by the minimum distance chain.
661  @param[in] _n: number of particles in _ptcl_list
662  @param[in] _ptcl: particle data array
663  @param[in] _G: gravtational constant
664  \return binary tree number
665  */
666  static void generateBinaryTree(BinaryTreeLocal _bins[], // make sure bins.size = n_members-1!
667  int _ptcl_list[], // make sure list.size = n_members!
668  const int _n,
669  Tptcl* _ptcl,
670  const Float _G) {
671  ASSERT(_n>=2);
672 
673  R2Index r2_list[_n];
674  // reorder _ptcl_list by minimum distance of each particles, and save square minimum distance r2min and index i in _ptcl_list (not particle index in _ptcl) in r2_list
675  calcMinDisList(_ptcl_list, r2_list, _n, _ptcl);
676  // sort r2_list by r2min
677  if(_n>2) std::sort(r2_list, r2_list+_n-1, pairLess);
678 
679  // tree root for each binary pair
680  BinaryTreeLocal* bin_host[_n];
681  for(auto &p : bin_host) p=NULL;
682 
683  // check binary from the closest pair to longest pair
684  for(int i=0; i<_n-1; i++) {
685  // get the pair index k in _ptcl_list with sorted r2_list, i represent the sorted r2min order.
686  int k = r2_list[i].second;
687  Tptcl* p[2];
688  int pindex[2]={-1,-1};
689  // if no tree root assign, set member 1 to particle and their host to current bins i
690  if(bin_host[k]==NULL) {
691 #ifdef BINARY_DEBUG
692  ASSERT(_ptcl_list[k]>=0);
693 #endif
694  p[0] = &_ptcl[_ptcl_list[k]];
695  pindex[0] = _ptcl_list[k];
696  bin_host[k] = &_bins[i];
697  }
698 
699  // if tree root already exist, member 1 assigned to tree root bin_host[k], tree members' bin_host -> current bins i
700  else {
701  p[0] = bin_host[k];
702  int ki = k;
703  while(bin_host[ki]==p[0]&&ki>=0) bin_host[ki--] = &_bins[i];
704  }
705 
706  // if no tree root assign, set member 2 to particle and their host to current bins i
707  if(bin_host[k+1]==NULL) {
708 #ifdef BINARY_DEBUG
709  ASSERT(_ptcl_list[k+1]>0);
710 #endif
711  p[1] = &_ptcl[_ptcl_list[k+1]];
712  pindex[1] = _ptcl_list[k+1];
713  bin_host[k+1] = &_bins[i];
714  }
715 
716  // if tree root already exist, member 2 assigned to tree root bin_host[k], tree members' bin_host -> current bins i
717  else {
718  p[1] = bin_host[k+1];
719  int ki = k+1;
720  while(bin_host[ki]==p[1]&&ki>=0) bin_host[ki++] = &_bins[i];
721  }
722 
723  // calculate binary parameter
724  _bins[i].setMembers(p[0], p[1], pindex[0], pindex[1]);
725  // calculate kepler orbit
726  _bins[i].calcOrbit(_G);
727  // calculate center-of-mass
728  _bins[i].calcCenterOfMass();
729  }
730 
731  // generate level branch
732  generateLevelBranch(_bins[_n-2],true);
733 
734 #ifdef BINARY_DEBUG
735  for(int i=0; i<_n; i++) ASSERT(bin_host[i]==&_bins[_n-2]); // check whether all bin_host point to the last of bins
736  ASSERT(_bins[_n-2].getMemberN() == _n); // makesure the particle number match the binarytree recored member
737 #endif
738  }
739 
741 
744  void getherBinaryTreeIter(BinaryTreeLocal _bin[]) {
745  copyDataIter(_bin, n_members-2);
746  }
747 
749  void calcCenterOfMass() {
750 #ifdef BINARY_DEBUG
751  ASSERT(member[0]!=NULL);
752  ASSERT(member[1]!=NULL);
753 #endif
754  Float m1 = member[0]->mass;
755  Float m2 = member[1]->mass;
756  this->pos[0] = m1*member[0]->pos[0] + m2*member[1]->pos[0];
757  this->pos[1] = m1*member[0]->pos[1] + m2*member[1]->pos[1];
758  this->pos[2] = m1*member[0]->pos[2] + m2*member[1]->pos[2];
759  this->vel[0] = m1*member[0]->vel[0] + m2*member[1]->vel[0];
760  this->vel[1] = m1*member[0]->vel[1] + m2*member[1]->vel[1];
761  this->vel[2] = m1*member[0]->vel[2] + m2*member[1]->vel[2];
762  this->mass = m1+m2;
763  this->pos[0] /=this->mass;
764  this->pos[1] /=this->mass;
765  this->pos[2] /=this->mass;
766  this->vel[0] /=this->mass;
767  this->vel[1] /=this->mass;
768  this->vel[2] /=this->mass;
769  }
770 
772  void calcSemiEccPeriod(const Float& _G) {
773  Tbinary::particleToSemiEccPeriod(*member[0], *member[1], _G);
774  }
775 
777  void calcOrbit(const Float _G) {
778  Tbinary::calcOrbit(*member[0], *member[1], _G);
779  }
780 
782  void calcParticles(const Float _G) {
783  Tbinary::calcParticles(*member[0], *member[1], _G);
784  }
785 
787  int calcMemberN() {
788 #ifdef BINARY_DEBUG
789  ASSERT(member[0]!=NULL);
790  ASSERT(member[1]!=NULL);
791 #endif
792  n_members = 0;
793  if (member_index[0]<0) n_members += ((BinaryTreeLocal*)member[0])->getMemberN();
794  else n_members++;
795  if (member_index[1]<0) n_members += ((BinaryTreeLocal*)member[1])->getMemberN();
796  else n_members++;
797  return n_members;
798  }
799 
801  int getMemberN() const {
802  return n_members;
803  }
804 
806  int calcMemberNIter() {
807 #ifdef BINARY_DEBUG
808  ASSERT(member[0]!=NULL);
809  ASSERT(member[1]!=NULL);
810 #endif
811  n_members = 0;
812  if (member_index[0]<0) n_members += (BinaryTreeLocal*)member[0]->calcMemberNIter();
813  else n_members++;
814  if (member_index[1]<0) n_members += (BinaryTreeLocal*)member[1]->calcMemberNIter();
815  else n_members++;
816  return n_members;
817  }
818 
820  template <class T>
821  using ProcessFunctionLeaf = void (*) (T&, Tptcl* &);
822 
824 
829  template <class T>
830  void processLeafIter(T& _dat, ProcessFunctionLeaf<T> _f) {
831 #ifdef BINARY_DEBUG
832  ASSERT(member[0]!=NULL);
833  ASSERT(member[1]!=NULL);
834 #endif
835  if (member_index[0]<0) ((BinaryTreeLocal*)member[0])->processLeafIter(_dat, _f);
836  else _f(_dat, member[0]);
837  if (member_index[1]<0) ((BinaryTreeLocal*)member[1])->processLeafIter(_dat, _f);
838  else _f(_dat, member[1]);
839  }
840 
842  template <class T>
843  using ProcessFunctionRoot = T (*) (T&, BinaryTreeLocal& );
844 
846 
851  template <class T>
852  T processRootIter(T& _dat, ProcessFunctionRoot<T> _f){
853  ASSERT(member[0]!=NULL);
854  ASSERT(member[1]!=NULL);
855  T dat_new = _f(_dat, *this);
856  for (int k=0; k<2; k++)
857  if (member_index[k]<0) dat_new = ((BinaryTreeLocal*)member[k])->processRootIter(dat_new, _f);
858  return dat_new;
859  }
860 
862 
869  template <class Troot, class Tleaf>
870  Troot processRootLeafIter(Troot& _dat_root, ProcessFunctionRoot<Troot> _f_root, Tleaf& _dat_leaf, ProcessFunctionLeaf<Tleaf> _f_leaf){
871  ASSERT(member[0]!=NULL);
872  ASSERT(member[1]!=NULL);
873  Troot dat_new = _f_root(_dat_root, *this);
874  for (int k=0; k<2; k++) {
875  if (member_index[k]<0) dat_new = ((BinaryTreeLocal*)member[k])->processRootLeafIter(dat_new, _f_root, _dat_leaf, _f_leaf);
876  else _f_leaf(_dat_leaf, member[k]);
877  }
878  return dat_new;
879  }
880 
882  template <class T, class Tr>
883  using ProcessFunctionTree = Tr (*) (T&, const Tr&, const Tr&, BinaryTreeLocal& );
884 
886 
893  template <class T, class Tr>
894  Tr processTreeIter(T& _dat, const Tr& _res1, const Tr& _res2, ProcessFunctionTree<T,Tr> _f){
895  ASSERT(member[0]!=NULL);
896  ASSERT(member[1]!=NULL);
897  Tr res[2] = {_res1, _res2};
898  for (int k=0; k<2; k++)
899  if (member_index[k]<0) res[k] = ((BinaryTreeLocal*)member[k])->processTreeIter(_dat, _res1, _res2, _f);
900  return _f(_dat, res[0], res[1], *this);
901  }
902 
903 
905 
909  int copyDataIter(BinaryTreeLocal _bin_out[], const int _i) {
910  int n_iter = 1;
911  int inow = _i; // array index to fill (backwards moving)
912  _bin_out[_i] = *this;
913  for (int k=0; k<2; k++) {
914  if (member_index[k]<0) {
915  inow--;
916  ASSERT(inow>=0);
917  _bin_out[_i].member[k] = &_bin_out[inow];
918  n_iter += ((BinaryTreeLocal*)member[k])->copyDataIter(_bin_out, inow);
919  }
920  }
921  return n_iter;
922  }
923 
925 
931  void setMembers(Tptcl* _p1, Tptcl* _p2, const int _i1, const int _i2) {
932  member[0] = _p1;
933  member[1] = _p2;
934  member_index[0] = _i1;
935  member_index[1] = _i2;
936  calcMemberN();
937  }
938 
940  Tptcl* getMember(const size_t i) const {
941  ASSERT(i<2);
942  return member[i];
943  }
944 
946  BinaryTreeLocal* getMemberAsTree(const size_t i) const {
947  ASSERT(i<2);
948  return (BinaryTreeLocal*)member[i];
949  }
950 
952  int getMemberIndex(const size_t i) const {
953  ASSERT(i<2);
954  return member_index[i];
955  }
956 
957  bool isMemberTree(const size_t i) const {
958  ASSERT(i<2);
959  return (member_index[i]<0);
960  }
961 
963  Tptcl* getLeftMember() const {
964  return member[0];
965  }
966 
968  BinaryTreeLocal* getLeftMemberAsTree() const {
969  return (BinaryTreeLocal*)member[0];
970  }
971 
973  Tptcl* getRightMember() const {
974  return member[1];
975  }
976 
978  BinaryTreeLocal* getRightMemberAsTree() const {
979  return (BinaryTreeLocal*)member[1];
980  }
981 
984  // */
985  //BinaryTreeLocal & operator = (const BinaryTreeLocal& _bin) {
986  // // copy everything first
987  // std::memcpy(this, &_bin, sizeof(BinaryTreeLocal));
988  // // get Different of address
989  // const BinaryTreeLocal* adr_diff = (const BinaryTreeLocal*)this - &_bin;
990  // // correct member address
991  // for (int k=0; k<2; k++) {
992  // if (_bin.member[k]!=NULL) {
993  // // if a member is another binarytree
994  // if (_bin.member_index[k]<0) {
995  // // Add the address difference to make the new member address consistent in the new binary tree array
996  // member[k] += adr_diff;
997  // }
998  // }
999  // }
1000  // return *this;
1001  //}
1002 
1004 
1008  static void printColumnTitle(std::ostream & _fout, const int _width=20) {
1009  Tptcl::printColumnTitle(_fout, _width);
1010  Tbinary::printColumnTitle(_fout, _width);
1011  }
1012 
1014 
1018  void printColumn(std::ostream & _fout, const int _width=20){
1019  Tptcl::printColumn(_fout, _width);
1020  Tbinary::printColumn(_fout, _width);
1021  }
1022 
1024 
1026  void writeAscii(std::ostream& _fout) const {
1027  Tptcl::writeAscii(_fout);
1028  Tbinary::writeAscii(_fout);
1029  }
1030 
1032 
1034  void readAscii(std::istream& _fin) {
1035  Tptcl::readAscii(_fin);
1036  Tbinary::readAscii(_fin);
1037  }
1038 
1040  void printBinaryTreeIter(std::ostream & _fout, const int _width=20){
1041  Tbinary::printColumn(_fout, _width);
1042  for (int k=0; k<2; k++) member[k]->printColumn(_fout,_width);
1043  for (int k=0; k<2; k++)
1044  if (member_index[k]<0) ((BinaryTreeLocal*)member[k])->printBinaryTreeIter(_fout, _width);
1045  }
1046 
1048  void printBinaryIter(std::ostream & _fout, const int _width=20){
1049  printColumn(_fout, _width);
1050  for (int k=0; k<2; k++) {
1051  if (member_index[k]<0) ((BinaryTreeLocal*)member[k])->printBinaryIter(_fout, _width);
1052  }
1053  }
1054  };
1055 }
COMM::Binary::readBinary
void readBinary(FILE *_fin)
read class data to file with binary format
Definition: binary_tree.h:432
COMM
Definition: binary_tree.h:8
COMM::BinaryTree::calcMemberN
int calcMemberN()
calc total number of members
Definition: binary_tree.h:787
COMM::Binary::rot_horizon
Float rot_horizon
Definition: binary_tree.h:17
COMM::Vector3< Float >
COMM::BinaryTree::level
int level
member pointer
Definition: binary_tree.h:510
COMM::BinaryTree::readAscii
void readAscii(std::istream &_fin)
read class data to file with ASCII format
Definition: binary_tree.h:1034
COMM::Binary::t_peri
Float t_peri
Definition: binary_tree.h:19
COMM::Matrix3
Definition: matrix.h:29
COMM::Binary::stab
Float stab
Definition: binary_tree.h:26
COMM::Binary::semiToPeriod
static Float semiToPeriod(const Float &_semi, const Float &_mtot, const Float &_G)
from semi-major axis calculate period;
Definition: binary_tree.h:176
COMM::BinaryTree::calcMemberNIter
int calcMemberNIter()
calc total number of members iteratively
Definition: binary_tree.h:806
COMM::BinaryTree< Tparticle >::ProcessFunctionTree
Tr(*)(T &, const Tr &, const Tr &, BinaryTreeLocal &) ProcessFunctionTree
iteratively tree processing function template
Definition: binary_tree.h:883
COMM::BinaryTree::BinaryTree
BinaryTree()
Definition: binary_tree.h:558
COMM::Binary::period
Float period
Definition: binary_tree.h:20
COMM::BinaryTree::isSameBranch
int isSameBranch(const BinaryTreeLocal &_bin) const
check whether a given binarytree data is in the same branch
Definition: binary_tree.h:638
COMM::Binary::ecca
Float ecca
Definition: binary_tree.h:21
Float.h
COMM::Binary::incline
Float incline
Definition: binary_tree.h:16
COMM::Binary::solveKepler
static void solveKepler(Binary &_bin, const Float _dt)
solve kepler orbit after dt
Definition: binary_tree.h:284
COMM::BinaryTree::getRightMemberAsTree
BinaryTreeLocal * getRightMemberAsTree() const
get right member as tree
Definition: binary_tree.h:978
COMM::Binary::calcParticlesEcca
void calcParticlesEcca(Tptcl &_p1, Tptcl &_p2, const Float _ecca, const Float _G) const
calculate two components from kepler Orbit with input eccentricity anomaly
Definition: binary_tree.h:332
NUMERIC_FLOAT_MAX
const Float NUMERIC_FLOAT_MAX
Definition: Float.h:29
COMM::Binary
Binary parameter class.
Definition: binary_tree.h:12
COMM::BinaryTree::printBinaryIter
void printBinaryIter(std::ostream &_fout, const int _width=20)
print binary function
Definition: binary_tree.h:1048
COMM::Binary::evolve
void evolve(const Float _dt)
Solve kepler motion for dt.
Definition: binary_tree.h:354
COMM::Matrix3::rotation
void rotation(const T I, const T OMEGA, const T omega)
Definition: matrix.h:51
COMM::BinaryTree::processRootIter
T processRootIter(T &_dat, ProcessFunctionRoot< T > _f)
Process root data and return result iteratively.
Definition: binary_tree.h:852
COMM::BinaryTree::branch
int branch
binary tree level, start from zero for root tree. Then each subtree increase level by one.
Definition: binary_tree.h:511
COMM::Binary::calcOrbit
void calcOrbit(const Tptcl &_p1, const Tptcl &_p2, const Float _G)
calculate kepler Orbit from particles
Definition: binary_tree.h:300
COMM::BinaryTree::getLeftMember
Tptcl * getLeftMember() const
get left member
Definition: binary_tree.h:963
COMM::BinaryTree::writeAscii
void writeAscii(std::ostream &_fout) const
write class data to file with ASCII format
Definition: binary_tree.h:1026
COMM::Binary::rot_self
Float rot_self
Definition: binary_tree.h:18
COMM::BinaryTree::copyDataIter
int copyDataIter(BinaryTreeLocal _bin_out[], const int _i)
copy a binary tree to a binaryTree array iteratively
Definition: binary_tree.h:909
matrix.h
COMM::Binary::semi
Float semi
Definition: binary_tree.h:14
COMM::BinaryTree::getherBinaryTreeIter
void getherBinaryTreeIter(BinaryTreeLocal _bin[])
collect BinaryTree data iteratively and save to a BinaryTree array
Definition: binary_tree.h:744
COMM::BinaryTree::printColumnTitle
static void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: binary_tree.h:1008
COMM::BinaryTree::generateBinaryTree
static void generateBinaryTree(BinaryTreeLocal _bins[], int _ptcl_list[], const int _n, Tptcl *_ptcl, const Float _G)
Generate kepler binary tree for a group of particles.
Definition: binary_tree.h:666
COMM::Binary::particleToSemiEcc
static void particleToSemiEcc(Float &_semi, Float &_ecc, Float &_r, Float &_rv, const Tpi &_p1, const Tpj &_p2, const Float _G)
position velocity to orbit semi-major axis and eccentricity
Definition: binary_tree.h:192
COMM::Binary::rotateToOriginalFrame
void rotateToOriginalFrame(Float *_vec)
rotate position vector from binary rest-frame to original frame based on three angles
Definition: binary_tree.h:340
COMM::BinaryTree::calcParticles
void calcParticles(const Float _G)
calculate particles from orbit
Definition: binary_tree.h:782
COMM::BinaryTree::getMemberN
int getMemberN() const
get total number of members
Definition: binary_tree.h:801
COMM::Vector3::y
T y
Definition: matrix.h:10
COMM::Binary::calcEccAnomaly
Float calcEccAnomaly(const Float _r)
calculate eccentric anomaly from separation (0-pi)
Definition: binary_tree.h:225
COMM::BinaryTree::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: binary_tree.h:1018
COMM::Binary::calcParticles
void calcParticles(Tptcl &_p1, Tptcl &_p2, const Float &_G)
calculate two components from kepler Orbit
Definition: binary_tree.h:311
COMM::Binary::m2
Float m2
Definition: binary_tree.h:23
Float
double Float
Definition: Float.h:25
COMM::Binary::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: binary_tree.h:403
COMM::Binary::periodToSemi
static Float periodToSemi(const Float &_period, const Float &_mtot, const Float &_G)
from period calculate semi-major axis
Definition: binary_tree.h:166
COMM::BinaryTree::calcSemiEccPeriod
void calcSemiEccPeriod(const Float &_G)
calculate Kepler orbit semi, ecc and period only
Definition: binary_tree.h:772
COMM::BinaryTree::calcOrbit
void calcOrbit(const Float _G)
calculate Kepler orbit from members
Definition: binary_tree.h:777
COMM::BinaryTree::isMemberTree
bool isMemberTree(const size_t i) const
Definition: binary_tree.h:957
COMM::Binary::m1
Float m1
Definition: binary_tree.h:22
COMM::BinaryTree::getMemberAsTree
BinaryTreeLocal * getMemberAsTree(const size_t i) const
get member as tree
Definition: binary_tree.h:946
COMM::BinaryTree::processRootLeafIter
Troot processRootLeafIter(Troot &_dat_root, ProcessFunctionRoot< Troot > _f_root, Tleaf &_dat_leaf, ProcessFunctionLeaf< Tleaf > _f_leaf)
Process root and leaf data and return result iteratively.
Definition: binary_tree.h:870
COMM::Binary::printColumnTitle
static void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: binary_tree.h:380
COMM::BinaryTree::processTreeIter
Tr processTreeIter(T &_dat, const Tr &_res1, const Tr &_res2, ProcessFunctionTree< T, Tr > _f)
Process tree data with extra dat iteratively (from bottom to top)
Definition: binary_tree.h:894
COMM::BinaryTree::printBinaryTreeIter
void printBinaryTreeIter(std::ostream &_fout, const int _width=20)
print binary and member information
Definition: binary_tree.h:1040
COMM::BinaryTree::getRightMember
Tptcl * getRightMember() const
get right member
Definition: binary_tree.h:973
COMM::BinaryTree::getMemberIndex
int getMemberIndex(const size_t i) const
get member index
Definition: binary_tree.h:952
COMM::Binary::print
void print(std::ostream &_os) const
Definition: binary_tree.h:359
COMM::PI
const Float PI
Definition: binary_tree.h:9
COMM::Binary::calcMeanAnomaly
static Float calcMeanAnomaly(const Float _ecca, const Float _ecc)
calculate mean anomaly from eccentric anomaly (0-pi)
Definition: binary_tree.h:242
COMM::Binary::orbitToParticle
static void orbitToParticle(Tptcl &_p1, Tptcl &_p2, const Binary &_bin, const Float &_ecca, const Float _G)
Orbit to position and velocity.
Definition: binary_tree.h:37
COMM::BinaryTree::setMembers
void setMembers(Tptcl *_p1, Tptcl *_p2, const int _i1, const int _i2)
set members
Definition: binary_tree.h:931
COMM::BinaryTree< Tparticle >::ProcessFunctionRoot
T(*)(T &, BinaryTreeLocal &) ProcessFunctionRoot
iteratively root processing function template
Definition: binary_tree.h:843
COMM::BinaryTree::calcCenterOfMass
void calcCenterOfMass()
calculate center-of-mass from members
Definition: binary_tree.h:749
COMM::Binary::calcSemiFromPeriod
void calcSemiFromPeriod(const Float &_G)
from period calculate semi-major axis
Definition: binary_tree.h:319
COMM::Binary::particleToOrbit
static void particleToOrbit(Binary &_bin, const Tptcl &_p1, const Tptcl &_p2, const Float _G)
position velocity to orbit
Definition: binary_tree.h:101
to_int
#define to_int(x)
Definition: Float.h:26
COMM::BinaryTree
Binary tree cell.
Definition: binary_tree.h:492
COMM::Binary::r
Float r
Definition: binary_tree.h:24
COMM::Binary::calcEccAnomaly
static Float calcEccAnomaly(const Float _mean_anomaly, const Float _ecc)
calculate eccentric anomaly from mean anomaly
Definition: binary_tree.h:253
COMM::BinaryTree::processLeafIter
void processLeafIter(T &_dat, ProcessFunctionLeaf< T > _f)
Process (bottom) leaf data with extra dat iteratively (from left to right)
Definition: binary_tree.h:830
COMM::BinaryTree::stableCheckIter
static Float & stableCheckIter(BinaryTreeLocal &_bin, const Float _t_crit)
Three-body stability function for hierarchical binary tree.
Definition: binary_tree.h:573
COMM::Binary::particleToSemiEccPeriod
void particleToSemiEccPeriod(const Tpi &_p1, const Tpj &_p2, const Float _G)
calcualte semi, ecc and period
Definition: binary_tree.h:213
COMM::Vector3::z
T z
Definition: matrix.h:10
COMM::Binary::writeAscii
void writeAscii(std::ostream &_fout) const
write class data to file with ASCII format
Definition: binary_tree.h:443
COMM::Binary::am
Vector3< Float > am
Definition: binary_tree.h:25
COMM::Binary::readAscii
void readAscii(std::istream &_fin)
read class data to file with ASCII format
Definition: binary_tree.h:464
COMM::Binary::writeBinary
void writeBinary(FILE *_fp) const
write class data to file with binary format
Definition: binary_tree.h:425
COMM::Binary::ecc
Float ecc
Definition: binary_tree.h:15
COMM::BinaryTree< Tparticle >::ProcessFunctionLeaf
void(*)(T &, Tparticle *&) ProcessFunctionLeaf
iteratively leaf processing function template
Definition: binary_tree.h:821
COMM::Vector3::x
T x
Definition: matrix.h:10
COMM::BinaryTree::getMember
Tptcl * getMember(const size_t i) const
get member
Definition: binary_tree.h:940
COMM::BinaryTree::getLeftMemberAsTree
BinaryTreeLocal * getLeftMemberAsTree() const
get left member as tree
Definition: binary_tree.h:968
COMM::BinaryTree::generateLevelBranch
static void generateLevelBranch(BinaryTreeLocal &_bin, const bool _reset)
generate level and branch
Definition: binary_tree.h:616