PeTar
N-body code for collisional gravitational systems
ar_interaction.hpp
Go to the documentation of this file.
1 #pragma once
2 #include <cmath>
3 #include "Common/Float.h"
4 #include "Common/binary_tree.h"
5 #include "changeover.hpp"
6 #include "AR/force.h"
7 #include "hard_ptcl.hpp"
8 #include "Hermite/hermite_particle.h"
9 #include "ar_perturber.hpp"
10 #include "two_body_tide.hpp"
11 #ifdef BSE_BASE
12 #include "bse_interface.h"
13 #endif
14 
16 class ARInteraction{
17 public:
18  typedef H4::ParticleH4<PtclHard> H4Ptcl;
19  Float eps_sq;
20  Float gravitational_constant;
21 #ifdef STELLAR_EVOLUTION
22  int stellar_evolution_option;
23  bool stellar_evolution_write_flag;
24  Float time_interrupt_max;
25 #ifdef BSE_BASE
26  BSEManager bse_manager;
27  TwoBodyTide tide;
28  std::ofstream fout_sse;
29  std::ofstream fout_bse;
30 
31  ARInteraction(): eps_sq(Float(-1.0)), gravitational_constant(Float(-1.0)),
32  stellar_evolution_option(1), stellar_evolution_write_flag(true), time_interrupt_max(NUMERIC_FLOAT_MAX),
33  bse_manager(), fout_sse(), fout_bse() {}
34 #else
35  ARInteraction(): eps_sq(Float(-1.0)), gravitational_constant(Float(-1.0)),
36  stellar_evolution_option(0), stellar_evolution_write_flag(true), time_interrupt_max(NUMERIC_FLOAT_MAX){}
37 #endif
38 #else
39  ARInteraction(): eps_sq(Float(-1.0)), gravitational_constant(Float(-1.0)) {}
40 #endif
41 
43 
45  bool checkParams() {
46  ASSERT(eps_sq>=0.0);
47  ASSERT(gravitational_constant>0.0);
48 #ifdef STELLAR_EVOLUTION
49  ASSERT(time_interrupt_max>=0.0);
50 #ifdef BSE_BASE
51  ASSERT(stellar_evolution_option==0 || (stellar_evolution_option==1 && bse_manager.checkParams()) || (stellar_evolution_option==2 && bse_manager.checkParams() && tide.checkParams()));
52  ASSERT(!stellar_evolution_write_flag||(stellar_evolution_write_flag&&fout_sse.is_open()));
53  ASSERT(!stellar_evolution_write_flag||(stellar_evolution_write_flag&&fout_bse.is_open()));
54 #endif
55 #endif
56  return true;
57  }
58 
60  void print(std::ostream & _fout) const{
61  _fout<<"eps_sq : "<<eps_sq<<std::endl
62  <<"G : "<<gravitational_constant<<std::endl;
63 #ifdef STELLAR_EVOLUTION
64  _fout<<"SE_opt : "<<stellar_evolution_option<<std::endl;
65 #endif
66  }
67 
69 
77  inline Float calcInnerAccPotAndGTKickInvTwo(AR::Force& _f1, AR::Force& _f2, Float& _epot, const PtclHard& _p1, const PtclHard& _p2) {
78  // acceleration
79  const Float mass1 = _p1.mass;
80  const Float* pos1 = &_p1.pos.x;
81 
82  const Float mass2 = _p2.mass;
83  const Float* pos2 = &_p2.pos.x;
84 
85  Float gm1 = gravitational_constant*mass1;
86  Float gm2 = gravitational_constant*mass2;
87  Float gm1m2 = gm1*mass2;
88 
89  Float dr[3] = {pos2[0] -pos1[0],
90  pos2[1] -pos1[1],
91  pos2[2] -pos1[2]};
92  Float r2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
93  Float inv_r = 1.0/sqrt(r2);
94  Float inv_r3 = inv_r*inv_r*inv_r;
95 
96  Float* acc1 = _f1.acc_in;
97  Float* acc2 = _f2.acc_in;
98 
99 #ifdef AR_CHANGEOVER
100  auto& ch1 = _p1.changeover;
101  auto& ch2 = _p2.changeover;
102 
103  Float r = r2*invr;
104  Float k = ChangeOver::calcAcc0WTwo(ch1,ch2,r);
105  Float kpot = ChangeOver::calcPotWTwo(ch1,ch2,r);
106 
107  Float gmor3_1 = gm2*inv_r3*k;
108  Float gmor3_2 = gm1*inv_r3*k;
109 
110  Float inv_rk = inv_r*kpot;
111  Float gm1or = gm1*inv_rk;
112  Float gm2or = gm2*inv_rk;
113  Float gm1m2or = gm1m2*inv_rk;
114 #else
115  Float gmor3_1 = gm2*inv_r3;
116  Float gmor3_2 = gm1*inv_r3;
117 
118  Float gm1or = gm1*inv_r;
119  Float gm2or = gm2*inv_r;
120  Float gm1m2or = gm1m2*inv_r;
121 #endif
122 
123  acc1[0] = gmor3_1 * dr[0];
124  acc1[1] = gmor3_1 * dr[1];
125  acc1[2] = gmor3_1 * dr[2];
126 
127  _f1.pot_in = -gm2or;
128 
129  acc2[0] = - gmor3_2 * dr[0];
130  acc2[1] = - gmor3_2 * dr[1];
131  acc2[2] = - gmor3_2 * dr[2];
132 
133  _f2.pot_in = -gm1or;
134 
135 
136 #ifdef AR_TTL
137  // trans formation function gradient
138 #ifdef AR_CHANGEOVER
139  Float gm1m2or3 = gm1m2*inv_r3*k;
140 #else
141  Float gm1m2or3 = gm1m2*inv_r3;
142 #endif
143  Float* gtgrad1 = _f1.gtgrad;
144  Float* gtgrad2 = _f2.gtgrad;
145  gtgrad1[0] = gm1m2or3 * dr[0];
146  gtgrad1[1] = gm1m2or3 * dr[1];
147  gtgrad1[2] = gm1m2or3 * dr[2];
148 
149  gtgrad2[0] = - gtgrad1[0];
150  gtgrad2[1] = - gtgrad1[1];
151  gtgrad2[2] = - gtgrad1[2];
152 #endif
153 
154  // potential energy
155  _epot = - gm1m2or;
156 
157  // transformation factor for kick
158  Float gt_kick_inv = gm1m2or;
159 
160  return gt_kick_inv;
161  }
162 
164 
171  inline Float calcInnerAccPotAndGTKickInv(AR::Force* _force, Float& _epot, const PtclHard* _particles, const int _n_particle) {
172  _epot = Float(0.0);
173  Float gt_kick_inv = Float(0.0);
174 
175  for (int i=0; i<_n_particle; i++) {
176  const Float massi = _particles[i].mass;
177  const Float* posi = &_particles[i].pos.x;
178  Float* acci = _force[i].acc_in;
179  acci[0] = acci[1] = acci[2] = Float(0.0);
180 
181 #ifdef AR_TTL
182  Float* gtgradi = _force[i].gtgrad;
183  gtgradi[0] = gtgradi[1] = gtgradi[2] = Float(0.0);
184 #endif
185 
186  Float poti = Float(0.0);
187  Float gtki = Float(0.0);
188 
189  for (int j=0; j<_n_particle; j++) {
190  if (i==j) continue;
191  const Float massj = _particles[j].mass;
192  const Float* posj = &_particles[j].pos.x;
193  Float dr[3] = {posj[0] -posi[0],
194  posj[1] -posi[1],
195  posj[2] -posi[2]};
196  Float r2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
197  Float inv_r = 1.0/sqrt(r2);
198  Float inv_r3 = inv_r*inv_r*inv_r;
199 #ifdef AR_CHANGEOVER
200  Float r = r2*inv_r;
201  const Float kpot = ChangeOver::calcPotWTwo(_particles[i].changeover, _particles[j].changeover, r);
202  const Float k = ChangeOver::calcAcc0WTwo(_particles[i].changeover, _particles[j].changeover, r);
203 
204  Float gmor3 = gravitational_constant*massj*inv_r3*k;
205  Float gmor = gravitational_constant*massj*inv_r*kpot;
206 #else
207  Float gmor3 = gravitational_constant*massj*inv_r3;
208  Float gmor = gravitational_constant*massj*inv_r;
209 #endif
210  acci[0] += gmor3 * dr[0];
211  acci[1] += gmor3 * dr[1];
212  acci[2] += gmor3 * dr[2];
213 
214 #ifdef AR_TTL
215  Float gmimjor3 = massi*gmor3;
216  gtgradi[0] += gmimjor3 * dr[0];
217  gtgradi[1] += gmimjor3 * dr[1];
218  gtgradi[2] += gmimjor3 * dr[2];
219 #endif
220 
221  poti -= gmor;
222  gtki += gmor;
223 
224  }
225  _epot += poti * massi;
226  gt_kick_inv += gtki * massi;
227  }
228  _epot *= 0.5;
229  gt_kick_inv *= 0.5;
230 
231  return gt_kick_inv;
232  }
233 
235 
242  void calcAccPert(AR::Force* _force, const PtclHard* _particles, const int _n_particle, const H4Ptcl& _particle_cm, const ARPerturber& _perturber, const Float _time) {
243  static const Float inv3 = 1.0 / 3.0;
244 
245  // perturber force
246  const int n_pert = _perturber.neighbor_address.getSize();
247  const int n_pert_single = _perturber.n_neighbor_single;
248  const int n_pert_group = _perturber.n_neighbor_group;
249 
250  if (n_pert>0) {
251 
252  Float time = _time;
253 
254  auto* pert_adr = _perturber.neighbor_address.getDataAddress();
255 
256  Float xp[n_pert][3], xcm[3], m[n_pert];
257  ChangeOver* changeover[n_pert_single];
258  H4::NBAdr<PtclHard>::Group* ptclgroup[n_pert_group];
259 
260  int n_single_count=0;
261  int n_group_count=0;
262  for (int j=0; j<n_pert; j++) {
263  H4::NBAdr<PtclHard>::Single* pertj;
264  int k; // index of predicted data
265  if (pert_adr[j].type==H4::NBType::group) {
266  pertj = &(((H4::NBAdr<PtclHard>::Group*)pert_adr[j].adr)->cm);
267  k = n_group_count + n_pert_single;
268  ptclgroup[n_group_count] = (H4::NBAdr<PtclHard>::Group*)pert_adr[j].adr;
269  n_group_count++;
270  }
271  else {
272  pertj = (H4::NBAdr<PtclHard>::Single*)pert_adr[j].adr;
273  k = n_single_count;
274  changeover[n_single_count] = &pertj->changeover;
275  n_single_count++;
276  }
277 
278  Float dt = time - pertj->time;
279  //ASSERT(dt>=-1e-7);
280  xp[k][0] = pertj->pos[0] + dt*(pertj->vel[0] + 0.5*dt*(pertj->acc0[0] + inv3*dt*pertj->acc1[0]));
281  xp[k][1] = pertj->pos[1] + dt*(pertj->vel[1] + 0.5*dt*(pertj->acc0[1] + inv3*dt*pertj->acc1[1]));
282  xp[k][2] = pertj->pos[2] + dt*(pertj->vel[2] + 0.5*dt*(pertj->acc0[2] + inv3*dt*pertj->acc1[2]));
283 
284 
285  m[k] = pertj->mass;
286  }
287  ASSERT(n_single_count == n_pert_single);
288  ASSERT(n_group_count == n_pert_group);
289 
290  Float dt = time - _particle_cm.time;
291  //ASSERT(dt>=0.0);
292  xcm[0] = _particle_cm.pos[0] + dt*(_particle_cm.vel[0] + 0.5*dt*(_particle_cm.acc0[0] + inv3*dt*_particle_cm.acc1[0]));
293  xcm[1] = _particle_cm.pos[1] + dt*(_particle_cm.vel[1] + 0.5*dt*(_particle_cm.acc0[1] + inv3*dt*_particle_cm.acc1[1]));
294  xcm[2] = _particle_cm.pos[2] + dt*(_particle_cm.vel[2] + 0.5*dt*(_particle_cm.acc0[2] + inv3*dt*_particle_cm.acc1[2]));
295 
296 
297  Float acc_pert_cm[3]={0.0, 0.0, 0.0};
298  Float mcm = 0.0;
299  // if (_perturber.need_resolve_flag) {
300  // calculate component perturbation
301  for (int i=0; i<_n_particle; i++) {
302  Float* acc_pert = _force[i].acc_pert;
303  Float& pot_pert = _force[i].pot_pert;
304  const auto& pi = _particles[i];
305  auto& chi = pi.changeover;
306  acc_pert[0] = acc_pert[1] = acc_pert[2] = Float(0.0);
307  pot_pert = 0.0;
308 
309  Float xi[3];
310  xi[0] = pi.pos[0] + xcm[0];
311  xi[1] = pi.pos[1] + xcm[1];
312  xi[2] = pi.pos[2] + xcm[2];
313 
314  // single perturber
315  for (int j=0; j<n_pert_single; j++) {
316  Float dr[3] = {xp[j][0] - xi[0],
317  xp[j][1] - xi[1],
318  xp[j][2] - xi[2]};
319  Float r2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2] + eps_sq;
320  Float r = sqrt(r2);
321  Float k = ChangeOver::calcAcc0WTwo(chi, *changeover[j], r);
322  Float r3 = r*r2;
323  Float gm = gravitational_constant*m[j];
324  Float gmor3 = gm/r3 * k;
325 
326  acc_pert[0] += gmor3 * dr[0];
327  acc_pert[1] += gmor3 * dr[1];
328  acc_pert[2] += gmor3 * dr[2];
329 
330  }
331  // group perturber
332  for (int j=n_pert_single; j<n_pert; j++) {
333  Float dr[3] = {xp[j][0] - xi[0],
334  xp[j][1] - xi[1],
335  xp[j][2] - xi[2]};
336  Float r2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2] + eps_sq;
337  Float r = sqrt(r2);
338  const int jk = j-n_pert_single;
339  auto* ptcl_mem = ptclgroup[jk]->getDataAddress();
340  Float mk = 0.0;
341  for (int k=0; k<ptclgroup[jk]->getSize(); k++) {
342  mk += ptcl_mem[k].mass * ChangeOver::calcAcc0WTwo(chi, ptcl_mem[k].changeover, r);
343  }
344  Float r3 = r*r2;
345  Float gmor3 = gravitational_constant*mk/r3;
346 
347  acc_pert[0] += gmor3 * dr[0];
348  acc_pert[1] += gmor3 * dr[1];
349  acc_pert[2] += gmor3 * dr[2];
350 
351  }
352 
353  acc_pert_cm[0] += pi.mass *acc_pert[0];
354  acc_pert_cm[1] += pi.mass *acc_pert[1];
355  acc_pert_cm[2] += pi.mass *acc_pert[2];
356 
357  mcm += pi.mass;
358 
359  }
360 //#ifdef AR_DEBUG
361 // ASSERT(abs(mcm-_particle_cm.mass)<1e-10);
362 //#endif
363 
364  // get cm perturbation (exclude soft pert)
365  acc_pert_cm[0] /= mcm;
366  acc_pert_cm[1] /= mcm;
367  acc_pert_cm[2] /= mcm;
368 
369  // remove cm. perturbation
370  for (int i=0; i<_n_particle; i++) {
371  Float* acc_pert = _force[i].acc_pert;
372  Float& pot_pert = _force[i].pot_pert;
373  const auto& pi = _particles[i];
374  acc_pert[0] -= acc_pert_cm[0];
375  acc_pert[1] -= acc_pert_cm[1];
376  acc_pert[2] -= acc_pert_cm[2];
377 
378  pot_pert -= acc_pert[0]*pi.pos[0] + acc_pert[1]*pi.pos[1] + acc_pert[2]*pi.pos[2];
379 
380 #ifdef SOFT_PERT
381  if(_perturber.soft_pert!=NULL) {
382  // avoid too large perturbation force if system is disruptted
383  if (pi.pos*pi.pos<pi.changeover.getRout()*pi.changeover.getRout()) {
384  _perturber.soft_pert->eval(acc_pert, pi.pos);
385  pot_pert += _perturber.soft_pert->evalPot(pi.pos);
386  }
387  }
388 #endif
389  }
390 
391  }
392  else {
393 #ifdef SOFT_PERT
394  if(_perturber.soft_pert!=NULL) {
395  for(int i=0; i<_n_particle; i++) {
396  Float* acc_pert = _force[i].acc_pert;
397  Float& pot_pert = _force[i].pot_pert;
398  const auto& pi = _particles[i];
399  acc_pert[0] = acc_pert[1] = acc_pert[2] = pot_pert = Float(0.0);
400  // avoid too large perturbation force if system is disruptted
401  if (pi.pos*pi.pos<pi.changeover.getRout()*pi.changeover.getRout()) {
402  _perturber.soft_pert->eval(acc_pert, pi.pos);
403  pot_pert += _perturber.soft_pert->evalPot(pi.pos);
404  }
405  }
406  }
407 #endif
408  }
409  }
410 
412 
422  Float calcAccPotAndGTKickInv(AR::Force* _force, Float& _epot, const PtclHard* _particles, const int _n_particle, const H4Ptcl& _particle_cm, const ARPerturber& _perturber, const Float _time) {
423  // inner force
424  Float gt_kick_inv;
425  if (_n_particle==2) gt_kick_inv = calcInnerAccPotAndGTKickInvTwo(_force[0], _force[1], _epot, _particles[0], _particles[1]);
426  else gt_kick_inv = calcInnerAccPotAndGTKickInv(_force, _epot, _particles, _n_particle);
427 
428  calcAccPert(_force, _particles, _n_particle, _particle_cm, _perturber, _time);
429 
430  return gt_kick_inv;
431  }
432 
433 
435  Float calcPertFromForce(const Float* _force, const Float _mp, const Float _mpert) {
436  Float force2 = _force[0]*_force[0]+_force[1]*_force[1]+_force[2]*_force[2];
437 #ifdef AR_SLOWDOWN_PERT_R4
438  return force2/(gravitational_constant*_mp*_mpert);
439 #else
440  Float force = sqrt(force2)/gravitational_constant;
441  return sqrt(force/(_mp*_mpert))*force;
442 #endif
443  }
444 
446  static Float calcPertFromBinary(const COMM::Binary& _bin) {
447  Float apo = _bin.semi*(1.0+_bin.ecc);
448  Float apo2 = apo*apo;
449 #ifdef AR_SLOWDOWN_PERT_R4
450  return (_bin.m1*_bin.m2)/(apo2*apo2);
451 #else
452  return (_bin.m1*_bin.m2)/(apo2*apo);
453 #endif
454  }
455 
457  static Float calcPertFromMR(const Float _r, const Float _mp, const Float _mpert) {
458  Float r2 = _r*_r;
459 #ifdef AR_SLOWDOWN_PERT_R4
460  return _mp*_mpert/(r2*r2);
461 #else
462  return (_mp*_mpert)/(r2*_r);
463 #endif
464  }
465 
466 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
467 
469  void calcSlowDownTimeScale(Float& _t_min_sq, const Float dv[3], const Float dr[3], const Float& r, const Float& gm) {
470 
471  Float r2 = r*r;
472  Float v2 = dv[0]*dv[0] + dv[1]*dv[1] + dv[2]*dv[2];
473  Float drdv = dr[0]*dv[0] + dr[1]*dv[1] + dr[2]*dv[2];
474 
475  Float semi = 1.0/(2.0/r - v2/gm);
476  //hyperbolic, directly use velocity v
477  if (semi<0)
478  _t_min_sq = std::min(_t_min_sq, r2/v2);
479  else {
480  Float ra_fact = (1 - r/semi);
481  Float e2 = drdv*drdv/(gm*semi) + ra_fact*ra_fact; // ecc^2
482  Float r_vrmax = semi*(1-e2);
483  if (r<r_vrmax) {
484  // avoid decrese of vr once the orbit pass, calculate vr max at cos(E)=e (r==semi*(1-e^2))
485  // vr_max = sqrt(er*(drdv^2*er + r*vcr2^2))/(G(m1+m2)r)
486  // = e*sqrt[G(m1+m2)/(a*(1-e^2)]
487  Float vrmax_sq = e2*gm/r_vrmax;
488  //Float rv2 = r*v2;
489  //Float er = 2*gm - rv2;
490  //Float vcr2 = gm - rv2;
491  //Float vrmax_sq = er*(drdv*drdv*er + r*vcr2*vcr2)/(gm*gm*r2);
492  _t_min_sq = std::min(_t_min_sq, semi*semi/vrmax_sq);
493  }
494  else {
495  // r/vr
496  Float rovr = r2/abs(drdv);
497  _t_min_sq = std::min(_t_min_sq, rovr*rovr);
498  }
499  }
500  }
501 
503 
509  void calcSlowDownPertOne(Float& _pert_out, Float& _t_min_sq, const PtclHard& pi, const PtclHard& pj) {
510  Float dr[3] = {pj.pos[0] - pi.pos[0],
511  pj.pos[1] - pi.pos[1],
512  pj.pos[2] - pi.pos[2]};
513  Float r2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
514  Float r = sqrt(r2);
515  _pert_out += calcPertFromMR(r, pi.mass, pj.mass);
516 
517 #ifdef AR_SLOWDOWN_TIMESCALE
518  Float dv[3] = {pj.vel[0] - pi.vel[0],
519  pj.vel[1] - pi.vel[1],
520  pj.vel[2] - pi.vel[2]};
521 
522  // identify whether hyperbolic or closed orbit
523  Float gm = gravitational_constant*(pi.mass+pj.mass);
524 
525  calcSlowDownTimeScale(_t_min_sq, dv, dr, r, gm);
526  // force dependent method
527  // min sqrt(r^3/(G m))
528  //Float gmor3 = (mp+mcm)*r*r2/(sdt->G*mp*mcm);
529  //sdt->trf2_min = std::min(sdt->trf2_min, gmor3);
530 #endif
531  }
532 
534 
541  void calcSlowDownPert(Float& _pert_out, Float& _t_min_sq, const Float& _time, const H4Ptcl& _particle_cm, const ARPerturber& _perturber) {
542  static const Float inv3 = 1.0 / 3.0;
543 
544  const int n_pert = _perturber.neighbor_address.getSize();
545 
546  if (n_pert>0) {
547 
548  auto* pert_adr = _perturber.neighbor_address.getDataAddress();
549 
550  Float xp[3], xcm[3];
551  Float dt = _time - _particle_cm.time;
552  //ASSERT(dt>=0.0);
553  xcm[0] = _particle_cm.pos[0] + dt*(_particle_cm.vel[0] + 0.5*dt*(_particle_cm.acc0[0] + inv3*dt*_particle_cm.acc1[0]));
554  xcm[1] = _particle_cm.pos[1] + dt*(_particle_cm.vel[1] + 0.5*dt*(_particle_cm.acc0[1] + inv3*dt*_particle_cm.acc1[1]));
555  xcm[2] = _particle_cm.pos[2] + dt*(_particle_cm.vel[2] + 0.5*dt*(_particle_cm.acc0[2] + inv3*dt*_particle_cm.acc1[2]));
556 
557  Float mcm = _particle_cm.mass;
558  auto& chi = _particle_cm.changeover;
559 
560 #ifdef AR_SLOWDOWN_TIMESCALE
561  // velocity dependent method
562  Float vp[3], vcm[3];
563 
564  vcm[0] = _particle_cm.vel[0] + dt*(_particle_cm.acc0[0] + 0.5*dt*_particle_cm.acc1[0]);
565  vcm[1] = _particle_cm.vel[1] + dt*(_particle_cm.acc0[1] + 0.5*dt*_particle_cm.acc1[1]);
566  vcm[2] = _particle_cm.vel[2] + dt*(_particle_cm.acc0[2] + 0.5*dt*_particle_cm.acc1[2]);
567 #endif
568 
569  for (int j=0; j<n_pert; j++) {
570  H4::NBAdr<PtclHard>::Single* pertj;
571  if (pert_adr[j].type==H4::NBType::group) pertj = &(((H4::NBAdr<PtclHard>::Group*)pert_adr[j].adr)->cm);
572  else pertj = (H4::NBAdr<PtclHard>::Single*)pert_adr[j].adr;
573 
574  Float dt = _time - pertj->time;
575  //ASSERT(dt>=0.0);
576  xp[0] = pertj->pos[0] + dt*(pertj->vel[0] + 0.5*dt*(pertj->acc0[0] + inv3*dt*pertj->acc1[0]));
577  xp[1] = pertj->pos[1] + dt*(pertj->vel[1] + 0.5*dt*(pertj->acc0[1] + inv3*dt*pertj->acc1[1]));
578  xp[2] = pertj->pos[2] + dt*(pertj->vel[2] + 0.5*dt*(pertj->acc0[2] + inv3*dt*pertj->acc1[2]));
579 
580  Float mj = pertj->mass;
581 
582  auto& chj = pertj->changeover;
583 
584  Float dr[3] = {xp[0] - xcm[0],
585  xp[1] - xcm[1],
586  xp[2] - xcm[2]};
587 
588  Float r2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2] + eps_sq;
589  Float r = sqrt(r2);
590  Float k = ChangeOver::calcAcc0WTwo(chi, chj, r);
591  _pert_out += calcPertFromMR(r, mcm, k*mj);
592 
593 #ifdef AR_SLOWDOWN_TIMESCALE
594  // velocity dependent method
595  vp[0] = pertj->vel[0] + dt*(pertj->acc0[0] + 0.5*dt*pertj->acc1[0]);
596  vp[1] = pertj->vel[1] + dt*(pertj->acc0[1] + 0.5*dt*pertj->acc1[1]);
597  vp[2] = pertj->vel[2] + dt*(pertj->acc0[2] + 0.5*dt*pertj->acc1[2]);
598 
599  Float dv[3] = {vp[0] - vcm[0],
600  vp[1] - vcm[1],
601  vp[2] - vcm[2]};
602 
603  // identify whether hyperbolic or closed orbit
604  Float gm = gravitational_constant*(mcm+mj);
605 
606  calcSlowDownTimeScale(_t_min_sq, dv, dr, r, gm);
607 #endif
608  }
609  }
610 
611  // add soft perturbation
612  _pert_out += _perturber.soft_pert_min;
613 
614  }
615 #endif
616 
618 
624  template <class Tparticle>
625  int modifyOneParticle(Tparticle& _p, const Float& _time_now, const Float& _time_end) {
626 #ifdef STELLAR_EVOLUTION
627  // sample of mass loss
628  //if (_p.time_interrupt<_time_end) {
629  // _p.dm = -_p.mass*1e-4;
630  // _p.mass +=_p.dm;
631  // _p.time_interrupt = _time_end+1e-5;
632  // return true;
633  //}
634 #ifdef BSE_BASE
635  // SSE/BSE stellar evolution
636  if (_p.time_interrupt<=_time_end&&stellar_evolution_option>0) {
637  ASSERT(bse_manager.checkParams());
638 
639  int modify_flag = 1;
640 
641  // time_record and time_interrupt have offsets, thus use difference to obtain true dt
642  Float dt = _time_end - _p.time_record;
643 
644  // evolve star
645  StarParameterOut output;
646  StarParameter star_bk = _p.star;
647  int event_flag = bse_manager.evolveStar(_p.star, output, dt);
648 
649  // error
650  if (event_flag<0) {
651  std::cerr<<"SSE Error: ID= "<<_p.id;
652  _p.star.print(std::cerr);
653  output.print(std::cerr);
654  std::cerr<<std::endl;
655  DATADUMP("dump_sse_error");
656  std::cout<<std::flush;
657  std::cerr<<std::flush;
658  abort();
659  }
660 
661  // if expected time not reach, record actually evolved time
662  double dt_miss = bse_manager.getDTMiss(output);
663  _p.time_record += dt-dt_miss;
664 
665  // estimate next time to check
666  _p.time_interrupt = std::min(_p.time_record + bse_manager.getTimeStepStar(_p.star), time_interrupt_max);
667 
668  // record mass change (if loss, negative)
669  double dm = bse_manager.getMassLoss(output);
670  _p.dm += dm;
671  if (dm==0.0) modify_flag = 0;
672 
673  // change mass in main data
674  _p.mass = bse_manager.getMass(_p.star);
675 
676  // set merger check radius
677  _p.radius = bse_manager.getMergerRadius(_p.star);
678 
679  // type change
680  if (stellar_evolution_write_flag&&event_flag>=1) {
681 #pragma omp critical
682  {
683  fout_sse<<"Type_change ";
684  //bse_manager.printTypeChange(fout_sse, _p.star, output);
685  fout_sse<<std::setw(WRITE_WIDTH)<<_p.id;
686  star_bk.printColumn(fout_sse, WRITE_WIDTH);
687  _p.star.printColumn(fout_sse, WRITE_WIDTH);
688  //output.printColumn(fout_sse, WRITE_WIDTH);
689  fout_sse<<std::endl;
690  }
691  }
692 
693  // add velocity change if exist
694  if (event_flag==2) {
695  double dv[3];
696  double dvabs=bse_manager.getVelocityChange(dv, output);
697  assert(dvabs>0);
698  for (int k=0; k<3; k++) _p.vel[k] += dv[k];
699  modify_flag = 2;
700  if (stellar_evolution_write_flag) {
701 #pragma omp critical
702  {
703  fout_sse<<"SN_kick "
704  <<std::setw(WRITE_WIDTH)<<_p.id
705  <<std::setw(WRITE_WIDTH)<<dvabs*bse_manager.vscale;
706  _p.star.printColumn(fout_sse, WRITE_WIDTH);
707  fout_sse<<std::endl;
708  }
709  }
710  }
711  // if mass become zero, set to unused for removing
712  if (_p.mass==0.0) {
713  _p.group_data.artificial.setParticleTypeToUnused(); // necessary to identify particle to remove
714  modify_flag = 3;
715  }
716 
717  return modify_flag;
718  }
719 
720 #endif // BSE_BASE
721 #endif // STELLAR_EVOLUTION
722  return 0;
723  }
724 
726 
731  int modifyAndInterruptIter(AR::InterruptBinary<PtclHard>& _bin_interrupt, AR::BinaryTree<PtclHard>& _bin) {
732  int modify_return = 0;
733 #ifdef STELLAR_EVOLUTION
734  int modify_branch[2];
735  if (_bin.getMemberN()>2) {
736  for (int k=0; k<2; k++) {
737  if (_bin.isMemberTree(k)) {
738  modify_branch[k] = modifyAndInterruptIter(_bin_interrupt, *_bin.getMemberAsTree(k));
739  modify_return = std::max(modify_return, modify_branch[k]);
740  }
741 #ifdef BSE_BASE
742  // if member is star, evolve single star using SSE
743  else if (stellar_evolution_option>0) {
744  ASSERT(bse_manager.checkParams());
745  modify_branch[k] = modifyOneParticle(*_bin.getMember(k), _bin.getMember(k)->time_record, _bin_interrupt.time_now);
746  modify_return = std::max(modify_return, modify_branch[k]);
747  // if status not set, set to change
748  if (modify_branch[k]>0&&_bin_interrupt.status == AR::InterruptStatus::none) {
749  _bin_interrupt.status = AR::InterruptStatus::change;
750  _bin_interrupt.adr = &_bin;
751  }
752  }
753 #endif
754  }
755  // ensure to record the root binary tree to include all changed members
756  if (modify_branch[0]>0&&modify_branch[1]>0) {
757  _bin_interrupt.adr = &_bin;
758  }
759  if (_bin_interrupt.status == AR::InterruptStatus::destroy) {
760  // if both branch has destroyed, set destroy status, otherwise set merge status
761  if (!(modify_branch[0]==3&&modify_branch[1]==3))
762  _bin_interrupt.status = AR::InterruptStatus::merge;
763  }
764  }
765  else {
766  auto* p1 = _bin.getLeftMember();
767  auto* p2 = _bin.getRightMember();
768 
769  COMM::Vector3<Float> pos_red(p2->pos[0] - p1->pos[0], p2->pos[1] - p1->pos[1], p2->pos[2] - p1->pos[2]);
770  COMM::Vector3<Float> vel_red(p2->vel[0] - p1->vel[0], p2->vel[1] - p1->vel[1], p2->vel[2] - p1->vel[2]);
771  Float drdv = pos_red * vel_red;
772 
773 #ifdef BSE_BASE
774  auto postProcess =[&](StarParameterOut* out, Float* pos_cm, Float*vel_cm, Float& semi, Float& ecc, int binary_type_final) {
775  // if status not set, set to change
776  if (_bin_interrupt.status == AR::InterruptStatus::none)
777  _bin_interrupt.status = AR::InterruptStatus::change;
778 
779  // set return flag >0
780  modify_return = 2;
781 
782  p1->time_record = _bin_interrupt.time_now - bse_manager.getDTMiss(out[0]);
783  p2->time_record = _bin_interrupt.time_now - bse_manager.getDTMiss(out[1]);
784 
785  // estimate next time to check
786  p1->time_interrupt = std::min(p1->time_record + bse_manager.getTimeStepBinary(p1->star, p2->star, semi, ecc, binary_type_final), time_interrupt_max);
787  p2->time_interrupt = p1->time_interrupt;
788 
789  // reset collision state since binary orbit changes
790  if (p1->getBinaryInterruptState()== BinaryInterruptState::collision)
791  p1->setBinaryInterruptState(BinaryInterruptState::none);
792  if (p2->getBinaryInterruptState()== BinaryInterruptState::collision)
793  p2->setBinaryInterruptState(BinaryInterruptState::none);
794 
795  // set binary status
796  //p1->setBinaryPairID(p2->id);
797  //p2->setBinaryPairID(p1->id);
798  p1->setBinaryInterruptState(static_cast<BinaryInterruptState>(binary_type_final));
799  p2->setBinaryInterruptState(static_cast<BinaryInterruptState>(binary_type_final));
800 
801  // record mass change (if loss, negative)
802  // dm is used to correct energy, thus must be correctly set, use += since it may change mass before merge
803  p1->dm += bse_manager.getMassLoss(out[0]);
804  p2->dm += bse_manager.getMassLoss(out[1]);
805 
806  // update masses
807  p1->mass = bse_manager.getMass(p1->star);
808  p2->mass = bse_manager.getMass(p2->star);
809 
810  // set merger check radius
811  p1->radius = bse_manager.getMergerRadius(p1->star);
812  p2->radius = bse_manager.getMergerRadius(p2->star);
813 
814 
815  bool mass_zero_flag = false;
816 
817  // if both mass becomes zero, set destroy state
818  if (p1->mass==0.0&&p2->mass==0.0) {
819  p1->group_data.artificial.setParticleTypeToUnused(); // necessary to identify particle to remove
820  p2->group_data.artificial.setParticleTypeToUnused(); // necessary to identify particle to remove
821  _bin_interrupt.status = AR::InterruptStatus::destroy;
822  modify_return = 3;
823  mass_zero_flag = true;
824  }
825  else {
826  // if mass become zero, set to unused for removing and merger status
827  if (p1->mass==0.0) {
828  p1->group_data.artificial.setParticleTypeToUnused(); // necessary to identify particle to remove
829  _bin_interrupt.status = AR::InterruptStatus::merge;
830  p1->setBinaryInterruptState(BinaryInterruptState::none);
831  p2->setBinaryInterruptState(BinaryInterruptState::none);
832  // set new particle position and velocity to be the original cm
833  for (int k=0; k<3; k++) {
834  p2->pos[k] = pos_cm[k];
835  p2->vel[k] = vel_cm[k];
836  }
837  modifyOneParticle(*p2, p2->time_record, _bin_interrupt.time_now);
838  mass_zero_flag = true;
839  }
840 
841  if (p2->mass==0.0) {
842  p2->group_data.artificial.setParticleTypeToUnused(); // necessary to identify particle to remove
843  _bin_interrupt.status = AR::InterruptStatus::merge;
844  p1->setBinaryInterruptState(BinaryInterruptState::none);
845  p2->setBinaryInterruptState(BinaryInterruptState::none);
846  // set new particle position and velocity to be the original cm
847  for (int k=0; k<3; k++) {
848  p1->pos[k] = pos_cm[k];
849  p1->vel[k] = vel_cm[k];
850  }
851  modifyOneParticle(*p1, p1->time_record, _bin_interrupt.time_now);
852  mass_zero_flag = true;
853  }
854 
855  // case when SN kick appears
856  bool kick_flag=false;
857  for (int k=0; k<2; k++) {
858  double dv[4];
859  auto* pk = _bin.getMember(k);
860  dv[3] = bse_manager.getVelocityChange(dv,out[k]);
861  if (dv[3]>0) {
862  kick_flag=true;
863 #pragma omp critical
864  {
865  fout_bse<<"SN_kick "
866  <<std::setw(WRITE_WIDTH)<<p1->id
867  <<std::setw(WRITE_WIDTH)<<p2->id
868  <<std::setw(WRITE_WIDTH)<<k+1
869  <<std::setw(WRITE_WIDTH)<<dv[3]*bse_manager.vscale;
870  pk->star.printColumn(fout_bse, WRITE_WIDTH);
871  fout_bse<<std::endl;
872  }
873  for (int k=0; k<3; k++) pk->vel[k] += dv[k];
874  }
875  }
876 
877  if (!kick_flag && !mass_zero_flag) {
878  // case for elliptic case
879  if (ecc>=0.0&&ecc<=1.0) {
880  // obtain full orbital parameters
881  //_bin.calcOrbit(gravitational_constant);
882  // update new period, ecc
883 //#pragma omp critical
884 // std::cerr<<"Event: "<<event_flag<<" "<<_bin.period<<" "<<period<<" "<<_bin.ecc<<" "<<ecc<<std::endl;
885  _bin.semi = semi;
886  ASSERT(_bin.semi>0);
887  _bin.ecc = ecc;
888  _bin.m1 = p1->mass;
889  _bin.m2 = p2->mass;
890  //if (((ecc-ecc_bk)/(1-ecc)>0.01||(period-period_bk)/period>1e-2)) {
891  // kepler orbit to particles using the same ecc anomaly
892  _bin.calcParticles(gravitational_constant);
893  p1->pos += _bin.pos;
894  p2->pos += _bin.pos;
895  p1->vel += _bin.vel;
896  p2->vel += _bin.vel;
897  //}
898  }
899  // in case of disruption but no kick
900  else {
901  // obtain full orbital parameters
902  // _bin.calcOrbit(gravitational_constant);
903  // assume energy no change
904  if(_bin.semi>0) _bin.semi = -_bin.semi;
905  _bin.ecc = ecc;
906  _bin.m1 = p1->mass;
907  _bin.m2 = p2->mass;
908  ASSERT(ecc>=1.0);
909  // kepler orbit to particles using the same ecc anomaly
910  //if ((ecc-ecc_bk)/ecc>1e-6) {
911  _bin.calcParticles(gravitational_constant);
912  p1->pos += _bin.pos;
913  p2->pos += _bin.pos;
914  p1->vel += _bin.vel;
915  p2->vel += _bin.vel;
916  //}
917  }
918  }
919  }
920  if (mass_zero_flag) {
921  DATADUMP("dump_merger");
922  }
923 
924  };
925 
926  bool check_flag = false;
927  if (stellar_evolution_option>0) {
928  int binary_type_p1 = static_cast<int>(p1->getBinaryInterruptState());
929  int binary_type_p2 = static_cast<int>(p2->getBinaryInterruptState());
930  int binary_type_init = 0;
931  if (binary_type_p1==binary_type_p2) binary_type_init = binary_type_p1;
932 
933  // check whether need to update
934  double time_check = std::min(p1->time_interrupt, p2->time_interrupt);
935  Float dt = _bin_interrupt.time_now - std::max(p1->time_record,p2->time_record);
936  if (time_check<=_bin_interrupt.time_now&&dt>0) check_flag = true;
937 
938  if (!check_flag) {
939  // pre simple check whether calling BSE is needed
940  Float t_record_min = std::min(p1->time_record,p2->time_record);
941 
942 
943  if (t_record_min<_bin_interrupt.time_now&&_bin.semi>0 && (!bse_manager.isMassTransfer(binary_type_init)) && (!bse_manager.isDisrupt(binary_type_init))) {
944  check_flag=bse_manager.isCallBSENeeded(p1->star, p2->star, _bin.semi, _bin.ecc, dt);
945  }
946  }
947 
948  if (check_flag) {
949  ASSERT(bse_manager.checkParams());
950  // record address of modified binary
951  _bin_interrupt.adr = &_bin;
952 
953  // first evolve two components to the same starting time
954  if (p1->time_record!=p2->time_record) {
955  if (p1->time_record<p2->time_record) {
956  p1->time_interrupt = p1->time_record;
957  modifyOneParticle(*p1, p1->time_record, p2->time_record);
958  }
959  else {
960  p2->time_interrupt = p2->time_record;
961  modifyOneParticle(*p2, p2->time_record, p1->time_record);
962  }
963  }
964  Float dt = _bin_interrupt.time_now - std::max(p1->time_record,p2->time_record);
965  ASSERT(dt>0);
966 
967  StarParameterOut out[2];
968  _bin.calcOrbit(gravitational_constant);
969  Float ecc = _bin.ecc;
970  //Float ecc_bk = ecc;
971  Float semi = _bin.semi;
972  //Float semi_bk =semi;
973  Float mtot = p1->mass+p2->mass;
974  Float period = _bin.period;
975  //Float period_bk = period;
976 
977  // backup c.m. information
978  Float pos_cm[3], vel_cm[3];
979  for (int k=0; k<3; k++) {
980  pos_cm[k] = (p1->mass*p1->pos[k] + p2->mass*p2->pos[k])/mtot;
981  vel_cm[k] = (p1->mass*p1->vel[k] + p2->mass*p2->vel[k])/mtot;
982  }
983  // backup star
984  StarParameter p1_star_bk = p1->star;
985  StarParameter p2_star_bk = p2->star;
986 
987  BinaryEvent bin_event;
988  // loop until the time_end reaches
989  int event_flag = bse_manager.evolveBinary(p1->star, p2->star, out[0], out[1], semi, period, ecc, bin_event, binary_type_init, dt);
990 
991  // error
992  if (event_flag<0) {
993  std::cerr<<"BSE Error! ";
994  std::cerr<<" ID="<<p1->id<<" "<<p2->id<<" ";
995  std::cerr<<" semi[R*]: "
996  <<_bin.semi*bse_manager.rscale
997  <<" ecc: "<<_bin.ecc
998  <<" period[days]: "<<_bin.period*bse_manager.tscale*bse_manager.year_to_day;
999  std::cerr<<" Init: Star1: ";
1000  p1_star_bk.print(std::cerr);
1001  std::cerr<<" Star2: ";
1002  p2_star_bk.print(std::cerr);
1003  std::cerr<<" final: Star1: ";
1004  p1->star.print(std::cerr);
1005  out[0].print(std::cerr);
1006  std::cerr<<" Star2: ";
1007  p2->star.print(std::cerr);
1008  out[1].print(std::cerr);
1009  std::cerr<<std::endl;
1010  DATADUMP("dump_bse_error");
1011  bse_manager.dumpRandConstant("bse.rand.par");
1012  std::cout<<std::flush;
1013  std::cerr<<std::flush;
1014  abort();
1015  }
1016 
1017  // check binary type and print event information
1018  int binary_type_final=0;
1019  int nmax = bin_event.getEventNMax();
1020  int binary_type_init = bin_event.getType(bin_event.getEventIndexInit());
1021  for (int i=0; i<nmax; i++) {
1022  int binary_type = bin_event.getType(i);
1023  if (binary_type>0) {
1024  bool first_event = (i==0);
1025  if (stellar_evolution_write_flag) {
1026  if ((first_event&&binary_type_init!=binary_type)||!first_event) {
1027  //if (!(binary_type_init==11&&(binary_type==3||binary_type==11))) {// avoid repeating printing Start Roche and BSS
1028 #pragma omp critical
1029  {
1030  bse_manager.printBinaryEventColumnOne(fout_bse, bin_event, i, WRITE_WIDTH);
1031  fout_bse<<std::setw(WRITE_WIDTH)<<p1->id
1032  <<std::setw(WRITE_WIDTH)<<p2->id
1033  <<std::setw(WRITE_WIDTH)<<drdv*bse_manager.rscale*bse_manager.vscale
1034  <<std::setw(WRITE_WIDTH)<<_bin.r*bse_manager.rscale;
1035  fout_bse<<std::endl;
1036  }
1037  }
1038  }
1039  //if (binary_type==10) {
1040  // DATADUMP("co_dump");
1041  // abort();
1042  //}
1043  //if (binary_type==3&&p1->time_record>=2499.0114383817) {
1044  // DATADUMP("re_dump");
1045  // abort();
1046  //}
1047 
1048  //if (vkick[3]>0||vkick[7]>0) event_flag = 3; // kick
1049  if (binary_type>0) event_flag = std::max(event_flag, 1); // type change
1050  else if (bse_manager.isMassTransfer(binary_type)) event_flag = std::max(event_flag, 2); // orbit change
1051  else if (bse_manager.isDisrupt(binary_type)) event_flag = std::max(event_flag, 3); // disrupt
1052  else if (bse_manager.isMerger(binary_type)) event_flag = std::max(event_flag, 4); // Merger
1053  binary_type_final = binary_type;
1054  }
1055  else if(binary_type<0) break;
1056  }
1057 
1058  // update semi
1059  mtot = bse_manager.getMass(p1->star) + bse_manager.getMass(p2->star);
1060  semi = COMM::Binary::periodToSemi(period, mtot, gravitational_constant);
1061 
1062  postProcess(out, pos_cm, vel_cm, semi, ecc, binary_type_final);
1063  }
1064  }
1065 #endif // BSE_BASE
1066 
1067  // dynamical merger and tide check
1068  if (_bin_interrupt.status!=AR::InterruptStatus::merge&&_bin_interrupt.status!=AR::InterruptStatus::destroy) {
1069 
1070  auto merge = [&](const Float& dr, const Float& t_peri, const Float& sd_factor) {
1071  _bin_interrupt.adr = &_bin;
1072 
1073 #ifdef BSE_BASE
1074  //Float m1_bk = p1->mass;
1075  //Float m2_bk = p2->mass;
1076  // backup original data for print
1077 
1078  // first evolve two components to the current time
1079  if (stellar_evolution_option>0) {
1080  if (p1->time_record<_bin_interrupt.time_now) {
1081  p1->time_interrupt = p1->time_record; // force SSE evolution
1082  modifyOneParticle(*p1, p1->time_record, _bin_interrupt.time_now);
1083  }
1084  if (p2->time_record<_bin_interrupt.time_now) {
1085  p2->time_interrupt = p2->time_record; // force SSE evolution
1086  modifyOneParticle(*p2, p2->time_record, _bin_interrupt.time_now);
1087  }
1088 
1089  ASSERT(p1->star.tphys==p2->star.tphys);
1090  //ASSERT(p1->star.mass>0&&p2->star.mass>0); // one may have SNe before merge
1091 
1092  StarParameter p1_star_bk, p2_star_bk;
1093  StarParameterOut out[2];
1094  Float pos_cm[3], vel_cm[3];
1095  Float mtot = p1->mass+p2->mass;
1096 
1097  for (int k=0; k<3; k++) {
1098  pos_cm[k] = (p1->mass*p1->pos[k] + p2->mass*p2->pos[k])/mtot;
1099  vel_cm[k] = (p1->mass*p1->vel[k] + p2->mass*p2->vel[k])/mtot;
1100  }
1101 
1102  p1_star_bk = p1->star;
1103  p2_star_bk = p2->star;
1104  // call BSE function to merge two stars
1105  Float semi = _bin.semi;
1106  Float ecc = _bin.ecc;
1107  bse_manager.merge(p1->star, p2->star, out[0], out[1], semi, ecc);
1108 
1109  postProcess(out, pos_cm, vel_cm, semi, ecc, 0);
1110  if (stellar_evolution_write_flag&&(p1->mass==0.0||p2->mass==0.0)) {
1111 #pragma omp critical
1112  {
1113  fout_bse<<"Dynamic_merge: "
1114  <<std::setw(WRITE_WIDTH)<<p1->id
1115  <<std::setw(WRITE_WIDTH)<<p2->id
1116  <<std::setw(WRITE_WIDTH)<<_bin.period*bse_manager.tscale*bse_manager.year_to_day
1117  <<std::setw(WRITE_WIDTH)<<_bin.semi*bse_manager.rscale
1118  <<std::setw(WRITE_WIDTH)<<_bin.ecc;
1119 #ifndef DYNAMIC_MERGER_LESS_OUTPUT
1120  fout_bse<<std::setw(WRITE_WIDTH)<<dr*bse_manager.rscale
1121  <<std::setw(WRITE_WIDTH)<<t_peri*bse_manager.tscale*bse_manager.year_to_day
1122  <<std::setw(WRITE_WIDTH)<<sd_factor;
1123 #endif
1124  // before
1125  p1_star_bk.printColumn(fout_bse, WRITE_WIDTH);
1126  p2_star_bk.printColumn(fout_bse, WRITE_WIDTH);
1127  // after
1128  p1->star.printColumn(fout_bse, WRITE_WIDTH);
1129  p2->star.printColumn(fout_bse, WRITE_WIDTH);
1130  fout_bse<<std::endl;
1131 
1132  DATADUMP("dump_merger");
1133 
1134  }
1135  }
1136  }
1137 #else //not BSE_BASE
1138  // print data
1139  std::cerr<<"Binary Merge: time: "<<_bin_interrupt.time_now<<std::endl;
1140  _bin.Binary::printColumnTitle(std::cerr);
1141  //PtclHard::printColumnTitle(std::cerr);
1142  //PtclHard::printColumnTitle(std::cerr);
1143  std::cerr<<std::endl;
1144  _bin.Binary::printColumn(std::cerr);
1145  //p1->printColumn(std::cerr);
1146  //p2->printColumn(std::cerr);
1147  std::cerr<<std::endl;
1148 
1149  // set return flag >0
1150  modify_return = 2;
1151 
1152  p1->time_record = _bin_interrupt.time_now;
1153  p2->time_record = _bin_interrupt.time_now;
1154 
1155  // new particle data
1156  Float mcm = p1->mass + p2->mass;
1157  for (int k=0; k<3; k++) {
1158  p1->pos[k] = (p1->mass*p1->pos[k] + p2->mass*p2->pos[k])/mcm;
1159  p1->vel[k] = (p1->mass*p1->vel[k] + p2->mass*p2->vel[k])/mcm;
1160  }
1161  p1->dm += p2->mass;
1162  p2->dm -= p2->mass;
1163 
1164  p1->mass = mcm;
1165  p2->mass = 0.0;
1166 
1167  p2->radius = 0.0;
1168 
1169  if (_bin_interrupt.status == AR::InterruptStatus::none)
1170  _bin_interrupt.status = AR::InterruptStatus::merge;
1171 
1172  // reset collision state since binary orbit changes
1173  p1->setBinaryInterruptState(BinaryInterruptState::none);
1174  p2->setBinaryInterruptState(BinaryInterruptState::none);
1175 
1176  p2->group_data.artificial.setParticleTypeToUnused(); // necessary to identify particle to remove
1177 #endif
1178  //p1->setBinaryPairID(0);
1179  //p2->setBinaryPairID(0);
1180  };
1181 
1182  // delayed merger
1183  if (p1->getBinaryInterruptState()== BinaryInterruptState::collision &&
1184  p2->getBinaryInterruptState()== BinaryInterruptState::collision &&
1185  (p1->time_interrupt<_bin_interrupt.time_now && p2->time_interrupt<_bin_interrupt.time_now) &&
1186  (p1->getBinaryPairID()==p2->id||p2->getBinaryPairID()==p1->id)) {
1187  Float dr[3] = {p1->pos[0] - p2->pos[0],
1188  p1->pos[1] - p2->pos[1],
1189  p1->pos[2] - p2->pos[2]};
1190  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
1191  merge(std::sqrt(dr2), 0.0, 1.0);
1192  }
1193  else {
1194  // check merger
1195  Float radius = p1->radius + p2->radius;
1196 #ifndef BSE_BASE
1197  // slowdown case
1198  if (_bin.slowdown.getSlowDownFactor()>1.0) {
1199  ASSERT(_bin.semi>0.0);
1200  Float drdv;
1201  _bin.particleToSemiEcc(_bin.semi, _bin.ecc, _bin.r, drdv, *_bin.getLeftMember(), *_bin.getRightMember(), gravitational_constant);
1202  Float peri = _bin.semi*(1 - _bin.ecc);
1203  if (peri<radius && p1->getBinaryPairID()!=p2->id&&p2->getBinaryPairID()!=p1->id) {
1204  Float ecc_anomaly = _bin.calcEccAnomaly(_bin.r);
1205  Float mean_anomaly = _bin.calcMeanAnomaly(ecc_anomaly, _bin.ecc);
1206  Float mean_motion = sqrt(gravitational_constant*_bin.mass/(fabs(_bin.semi*_bin.semi*_bin.semi)));
1207  Float t_peri = mean_anomaly/mean_motion;
1208  if (drdv<0 && t_peri<_bin_interrupt.time_end-_bin_interrupt.time_now) {
1209  Float dr[3] = {p1->pos[0] - p2->pos[0],
1210  p1->pos[1] - p2->pos[1],
1211  p1->pos[2] - p2->pos[2]};
1212  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
1213  merge(std::sqrt(dr2), t_peri, _bin.slowdown.getSlowDownFactor());
1214  }
1215  else if (_bin.semi>0||(_bin.semi<0&&drdv<0)) {
1216  p1->setBinaryPairID(p2->id);
1217  p2->setBinaryPairID(p1->id);
1218  p1->setBinaryInterruptState(BinaryInterruptState::collision);
1219  p2->setBinaryInterruptState(BinaryInterruptState::collision);
1220  p1->time_interrupt = std::min(_bin_interrupt.time_now + drdv<0 ? t_peri : (_bin.period - t_peri), time_interrupt_max);
1221  p2->time_interrupt = p1->time_interrupt;
1222 
1223  }
1224  }
1225  }
1226  else { // no slowdown case, check separation directly
1227  Float dr[3] = {p1->pos[0] - p2->pos[0],
1228  p1->pos[1] - p2->pos[1],
1229  p1->pos[2] - p2->pos[2]};
1230  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
1231  if (dr2<radius*radius) merge(std::sqrt(dr2), 0.0, 1.0);
1232  }
1233 #else
1234  // in bse case, handle binary merger in bse, only check hyperbolic merger
1235  if (_bin.semi<0.0) {
1236  Float dr[3] = {p1->pos[0] - p2->pos[0],
1237  p1->pos[1] - p2->pos[1],
1238  p1->pos[2] - p2->pos[2]};
1239  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
1240  if (dr2<radius*radius) merge(std::sqrt(dr2), 0.0, 1.0);
1241  }
1242 #endif
1243  }
1244 
1245 
1246 #ifdef BSE_BASE
1247  // tide energy loss
1248  if (stellar_evolution_option==2 && p1->mass>0 && p2->mass>0) {
1249  if (drdv<0) { // when two star approach each other; reset tide status
1250  if (p1->getBinaryInterruptState() == BinaryInterruptState::tide) {
1251  p1->setBinaryInterruptState(BinaryInterruptState::none);
1252  }
1253  if (p2->getBinaryInterruptState() == BinaryInterruptState::tide) {
1254  p2->setBinaryInterruptState(BinaryInterruptState::none);
1255  }
1256  }
1257  else { // modify orbit based on energy loss
1258  int binary_type_p1 = static_cast<int>(p1->getBinaryInterruptState());
1259  int binary_type_p2 = static_cast<int>(p2->getBinaryInterruptState());
1260  long long int pair_id1 = p1->getBinaryPairID();
1261  long long int pair_id2 = p2->getBinaryPairID();
1262  bool tide_flag = true;
1263  if ((binary_type_p1 != binary_type_p2) || (pair_id1 != p2->id) || (pair_id2 != p1->id)) tide_flag = false;
1264  else if (bse_manager.isMassTransfer(binary_type_p1)
1265  || bse_manager.isMerger(binary_type_p1)
1266  || bse_manager.isDisrupt(binary_type_p1)
1267  || binary_type_p1 == 14)
1268  tide_flag = false;
1269 
1270  bool change_flag=false;
1271  if (tide_flag) {
1272  Float poly_type1=0, poly_type2=0;
1273  Float Etid=0, Ltid=0;
1274  Float semi = _bin.semi;
1275  Float ecc = _bin.ecc;
1276  Float rad1 = bse_manager.getStellarRadius(p1->star);
1277  Float rad2 = bse_manager.getStellarRadius(p2->star);
1278  if (p1->star.kw>=10&&p1->star.kw<15&&p2->star.kw>=10&&p2->star.kw<15) {
1279  if (_bin.semi<0) {
1280  bool merge_flag = tide.evolveOrbitHyperbolicGW(_bin, Etid, Ltid);
1281  if (merge_flag) {
1282  Float dr[3] = {p1->pos[0] - p2->pos[0],
1283  p1->pos[1] - p2->pos[1],
1284  p1->pos[2] - p2->pos[2]};
1285  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
1286  merge(std::sqrt(dr2), 0.0, 1.0);
1287  }
1288  else change_flag = true;
1289  }
1290  }
1291  else if (std::min(p1->star.kw, p2->star.kw)<13) {
1292  poly_type1 = (p1->star.kw<=2) ? 3.0 : 1.5;
1293  poly_type2 = (p2->star.kw<=2) ? 3.0 : 1.5;
1294  Etid = tide.evolveOrbitDynamicalTide(_bin, rad1, rad2, poly_type1, poly_type2);
1295  change_flag = (Etid>0);
1296  // for slowdown case, repeating tide effect based on slowdown factor
1297  Float sd_factor_ext = _bin.slowdown.getSlowDownFactor() - 1.5;
1298  if (change_flag && sd_factor_ext>0) {
1299  for (Float k=0; k<sd_factor_ext; k=k+1.0) {
1300  Float etid_k = tide.evolveOrbitDynamicalTide(_bin, rad1, rad2, poly_type1, poly_type2);
1301  if (etid_k==0) break;
1302  Etid += etid_k;
1303  }
1304  }
1305  }
1306 
1307  if (change_flag) {
1308 
1309  _bin_interrupt.adr = &_bin;
1310 
1311  // if status not set, set to change
1312  if (_bin_interrupt.status == AR::InterruptStatus::none)
1313  _bin_interrupt.status = AR::InterruptStatus::change;
1314  _bin.calcParticles(gravitational_constant);
1315  p1->pos += _bin.pos;
1316  p2->pos += _bin.pos;
1317  p1->vel += _bin.vel;
1318  p2->vel += _bin.vel;
1319 
1320  p1->setBinaryPairID(p2->id);
1321  p2->setBinaryPairID(p1->id);
1322  p1->setBinaryInterruptState(BinaryInterruptState::tide);
1323  p2->setBinaryInterruptState(BinaryInterruptState::tide);
1324 
1325  modify_return = 2;
1326 
1327 #pragma omp critical
1328  {
1329  fout_bse<<"Tide "
1330  <<std::setw(WRITE_WIDTH)<<_bin_interrupt.time_now
1331  <<std::setw(WRITE_WIDTH)<<p1->id
1332  <<std::setw(WRITE_WIDTH)<<p2->id
1333  <<std::setw(WRITE_WIDTH)<<pair_id1
1334  <<std::setw(WRITE_WIDTH)<<pair_id2
1335  <<std::setw(WRITE_WIDTH)<<binary_type_p1
1336  <<std::setw(WRITE_WIDTH)<<binary_type_p2
1337  <<std::setw(WRITE_WIDTH)<<poly_type1
1338  <<std::setw(WRITE_WIDTH)<<poly_type2
1339  <<std::setw(WRITE_WIDTH)<<drdv
1340  <<std::setw(WRITE_WIDTH)<<semi //old
1341  <<std::setw(WRITE_WIDTH)<<ecc //old
1342  <<std::setw(WRITE_WIDTH)<<Etid
1343  <<std::setw(WRITE_WIDTH)<<Ltid;
1344  _bin.BinarySlowDown::printColumn(fout_bse, WRITE_WIDTH);
1345  p1->star.printColumn(fout_bse, WRITE_WIDTH);
1346  p2->star.printColumn(fout_bse, WRITE_WIDTH);
1347  fout_bse<<std::endl;
1348  }
1349 
1350  }
1351  }
1352  }
1353  }
1354 #endif // BSE_BASE
1355  }
1356  }
1357 #endif
1358  return modify_return;
1359  }
1360 
1361 #ifndef AR_TTL
1362 
1366  Float calcGTDriftInv(Float _ekin_minus_etot) {
1367  return _ekin_minus_etot;
1368  }
1369 
1371 
1374  Float calcH(Float _ekin_minus_etot, Float _epot) {
1375  return log(_ekin_minus_etot) - log(-_epot);
1376  }
1377 #endif
1378 
1380 
1382  void writeBinary(FILE *_fp) const {
1383  fwrite(&eps_sq, sizeof(Float),1,_fp);
1384  fwrite(&gravitational_constant, sizeof(Float),1,_fp);
1385 #ifdef STELLAR_EVOLUTION
1386  fwrite(&stellar_evolution_option, sizeof(int),1,_fp);
1387  fwrite(&stellar_evolution_write_flag, sizeof(bool),1,_fp);
1388 #endif
1389  }
1390 
1392 
1394  void readBinary(FILE *_fin) {
1395  size_t rcount = fread(&eps_sq, sizeof(Float),1,_fin);
1396  rcount += fread(&gravitational_constant, sizeof(Float),1,_fin);
1397  if (rcount<2) {
1398  std::cerr<<"Error: Data reading fails! requiring data number is 2, only obtain "<<rcount<<".\n";
1399  abort();
1400  }
1401 #ifdef STELLAR_EVOLUTION
1402  rcount += fread(&stellar_evolution_option, sizeof(int),1,_fin);
1403  rcount += fread(&stellar_evolution_write_flag, sizeof(bool),1,_fin);
1404  if (rcount<4) {
1405  std::cerr<<"Error: Data reading fails! requiring data number is 4, only obtain "<<rcount<<".\n";
1406  abort();
1407  }
1408 #endif
1409  }
1410 };
1411 
ARInteraction::calcAccPert
void calcAccPert(AR::Force *_force, const PtclHard *_particles, const int _n_particle, const H4Ptcl &_particle_cm, const ARPerturber &_perturber, const Float _time)
(Necessary) calculate acceleration from perturber and the perturbation factor for slowdown calculatio...
Definition: ar_interaction.hpp:242
TwoBodyTide::checkParams
bool checkParams()
check whether parameters values are correct
Definition: two_body_tide.hpp:19
BSEManager::merge
void merge(StarParameter &_star1, StarParameter &_star2, StarParameterOut &_out1, StarParameterOut &_out2, double &_semi, double &_ecc)
merge two star using mix function, star 2 will becomes zero mass
Definition: bse_interface.h:1612
BSEManager::year_to_day
const double year_to_day
velocity scaling factor from NB to km/s
Definition: bse_interface.h:1063
TwoBodyTide::evolveOrbitHyperbolicGW
bool evolveOrbitHyperbolicGW(TBinary &_bin, Float &_Etid, Float &_Ltid)
evolve hyperbolic orbit based on GW effect
Definition: two_body_tide.hpp:31
BSEManager::rscale
double rscale
time scaling factor from NB to Myr (t[Myr]=t[NB]*tscale)
Definition: bse_interface.h:1060
ARInteraction::modifyAndInterruptIter
int modifyAndInterruptIter(AR::InterruptBinary< PtclHard > &_bin_interrupt, AR::BinaryTree< PtclHard > &_bin)
(Necessary) modify the orbits and interrupt check
Definition: ar_interaction.hpp:731
BSEManager::evolveBinary
int evolveBinary(StarParameter &_star1, StarParameter &_star2, StarParameterOut &_out1, StarParameterOut &_out2, double &_semi, double &_period, double &_ecc, BinaryEvent &_bse_event, const int &_binary_init_type, const double _dt_nb)
call evolv2 for a binary
Definition: bse_interface.h:1450
BSEManager::printBinaryEventColumnOne
void printBinaryEventColumnOne(std::ostream &_fout, const BinaryEvent &_bin_event, const int k, const int _width=20, const bool print_type_name=true)
print binary event one in column
Definition: bse_interface.h:1384
ARPerturber::soft_pert
TidalTensor * soft_pert
Definition: ar_perturber.hpp:14
BSEManager::isDisrupt
bool isDisrupt(const int _binary_type)
Definition: bse_interface.h:1221
bse_interface.h
BSEManager::isMassTransfer
bool isMassTransfer(const int _binary_type)
Definition: bse_interface.h:1208
galpy_pot_movie.dt
dt
Definition: galpy_pot_movie.py:134
BinaryInterruptState
BinaryInterruptState
Definition: particle_base.hpp:14
two_body_tide.hpp
StarParameter::printColumn
void printColumn(std::ostream &_fout, const int _width=20) const
print data of class members using column style
Definition: bse_interface.h:345
StarParameterOut::print
void print(std::ostream &fout) const
for print in one line
Definition: bse_interface.h:439
ARInteraction::eps_sq
Float eps_sq
Definition: ar_interaction.hpp:19
BinaryEvent
BSE based code event recorder class.
Definition: bse_interface.h:493
BSEManager::getStellarRadius
double getStellarRadius(StarParameter &_star)
get stellar radius in NB unit
Definition: bse_interface.h:1329
TwoBodyTide
Two Body Tide effct.
Definition: two_body_tide.hpp:9
BSEManager::getMass
double getMass(StarParameter &_star)
get current mass in NB unit
Definition: bse_interface.h:1313
BSEManager::getMassLoss
double getMassLoss(StarParameterOut &_out)
get mass loss in NB unit
Definition: bse_interface.h:1318
ARInteraction::print
void print(std::ostream &_fout) const
print parameters
Definition: ar_interaction.hpp:60
ParticleBase::pos
PS::F64vec pos
Definition: particle_base.hpp:24
ParticleBase::vel
PS::F64vec vel
Definition: particle_base.hpp:25
ARInteraction::calcPertFromMR
static Float calcPertFromMR(const Float _r, const Float _mp, const Float _mpert)
calculate perturbation from distance to perturber and masses of particle and perturber
Definition: ar_interaction.hpp:457
ARInteraction::checkParams
bool checkParams()
(Necessary) check whether publicly initialized parameters are correctly set
Definition: ar_interaction.hpp:45
PtclHard
Definition: hard_ptcl.hpp:5
ARInteraction::gravitational_constant
Float gravitational_constant
softening parameter
Definition: ar_interaction.hpp:20
BinaryInterruptState::none
@ none
PIKG::min
T min(const Vector3< T > &v)
Definition: pikg_vector.hpp:150
BinaryEvent::getEventNMax
int getEventNMax() const
Maximum event number that can be recored in one call of evolv2.
Definition: bse_interface.h:534
BSEManager::checkParams
bool checkParams()
Definition: bse_interface.h:1090
ARInteraction::calcH
Float calcH(Float _ekin_minus_etot, Float _epot)
(Necessary) calculate the time transformed Hamiltonian
Definition: ar_interaction.hpp:1374
BSEManager::getDTMiss
double getDTMiss(StarParameterOut &_out)
get the difference of required finishing time and actually evolved time in NB unit
Definition: bse_interface.h:1347
StarParameterOut
SSE/BSE based code star parameter for output.
Definition: bse_interface.h:393
DATADUMP
#define DATADUMP(expr)
Definition: hard_assert.hpp:227
ARInteraction::calcPertFromForce
Float calcPertFromForce(const Float *_force, const Float _mp, const Float _mpert)
calculate perturbation from c.m. acceleration
Definition: ar_interaction.hpp:435
ChangeOver
Changeover function class.
Definition: changeover.hpp:7
Ptcl::changeover
ChangeOver changeover
Definition: ptcl.hpp:41
TwoBodyTide::evolveOrbitDynamicalTide
Float evolveOrbitDynamicalTide(TBinary &_bin, const Float &rad1, const Float &rad2, const Float &poly_type1, const Float &poly_type2)
evolve two-body orbit based on dynamical tide implementation from Alessandro Alberto Trani
Definition: two_body_tide.hpp:102
BSEManager::getMergerRadius
double getMergerRadius(StarParameter &_star)
get merger radius in NB unit
Definition: bse_interface.h:1323
ARInteraction::calcAccPotAndGTKickInv
Float calcAccPotAndGTKickInv(AR::Force *_force, Float &_epot, const PtclHard *_particles, const int _n_particle, const H4Ptcl &_particle_cm, const ARPerturber &_perturber, const Float _time)
(Necessary) calculate acceleration from internal members and perturbers
Definition: ar_interaction.hpp:422
BSEManager::isMerger
bool isMerger(const int _binary_type)
notice kick priority is higher than others
Definition: bse_interface.h:1217
changeover.hpp
BSEManager::isCallBSENeeded
bool isCallBSENeeded(StarParameter &_star1, StarParameter &_star2, double &_semi, double &_ecc, double &_dt)
a simple check to determine whether call BSE is necessary by given a reference of time step
Definition: bse_interface.h:1772
ChangeOver::calcPotWTwo
static Float calcPotWTwo(const ChangeOver &_ch1, const ChangeOver &_ch2, const Float &_dr)
calculate changeover function Pot by selecting maximum rout
Definition: changeover.hpp:366
ChangeOver::calcAcc0WTwo
static Float calcAcc0WTwo(const ChangeOver &_ch1, const ChangeOver &_ch2, const Float &_dr)
calculate changeover function Acc0 by selecting maximum rout
Definition: changeover.hpp:372
BSEManager::evolveStar
int evolveStar(StarParameter &_star, StarParameterOut &_out, const double _dt, bool _unit_in_myr=false)
call SSE evolv1 for single star
Definition: bse_interface.h:1412
ARInteraction::readBinary
void readBinary(FILE *_fin)
read class data to file with binary format
Definition: ar_interaction.hpp:1394
BSEManager::tscale
double tscale
metallicity parameters
Definition: bse_interface.h:1059
ARInteraction::ARInteraction
ARInteraction()
Definition: ar_interaction.hpp:39
TidalTensor::eval
void eval(PS::F64 *acc, const PS::F64vec &pos) const
Definition: tidal_tensor.hpp:289
WRITE_WIDTH
const int WRITE_WIDTH
Definition: bse_test.cxx:10
BSEManager
SSE/BSE interface manager.
Definition: bse_interface.h:1053
ParticleBase::mass
PS::F64 mass
Definition: particle_base.hpp:23
PIKG::max
T max(const Vector3< T > &v)
Definition: pikg_vector.hpp:143
BSEManager::getVelocityChange
double getVelocityChange(double *dv, StarParameterOut &_out)
get velocity change in NB unit
Definition: bse_interface.h:1399
ARInteraction::calcPertFromBinary
static Float calcPertFromBinary(const COMM::Binary &_bin)
calculate perturbation from binary tree
Definition: ar_interaction.hpp:446
BinaryEvent::getEventIndexInit
int getEventIndexInit() const
The index of initial event in record array (last one)
Definition: bse_interface.h:539
ar_perturber.hpp
TidalTensor::evalPot
PS::F64 evalPot(const PS::F64vec &pos) const
Definition: tidal_tensor.hpp:344
ARPerturber
Perturber class for AR integration.
Definition: ar_perturber.hpp:11
BinaryInterruptState::collision
@ collision
BSEManager::vscale
double vscale
mass scaling factor from NB to Msun
Definition: bse_interface.h:1062
ARInteraction::calcInnerAccPotAndGTKickInvTwo
Float calcInnerAccPotAndGTKickInvTwo(AR::Force &_f1, AR::Force &_f2, Float &_epot, const PtclHard &_p1, const PtclHard &_p2)
(Necessary) calculate inner member acceleration, potential and inverse time transformation function g...
Definition: ar_interaction.hpp:77
ASSERT
#define ASSERT(expr)
Definition: hard_assert.hpp:238
ARInteraction
AR interaction clas.
Definition: ar_interaction.hpp:16
ARInteraction::H4Ptcl
H4::ParticleH4< PtclHard > H4Ptcl
Definition: ar_interaction.hpp:18
hard_ptcl.hpp
ARInteraction::writeBinary
void writeBinary(FILE *_fp) const
write class data to file with binary format
Definition: ar_interaction.hpp:1382
ARInteraction::calcGTDriftInv
Float calcGTDriftInv(Float _ekin_minus_etot)
(Necessary) calcualte the inverse time transformation factor for drift
Definition: ar_interaction.hpp:1366
StarParameter::print
void print(std::ostream &fout) const
for print in one line
Definition: bse_interface.h:309
BSEManager::getTimeStepBinary
double getTimeStepBinary(StarParameter &_star1, StarParameter &_star2, double &_semi, double &_ecc, int &_binary_type)
call BSE evolv2 for a binary
Definition: bse_interface.h:1719
ARInteraction::modifyOneParticle
int modifyOneParticle(Tparticle &_p, const Float &_time_now, const Float &_time_end)
(Necessary) modify one particle function
Definition: ar_interaction.hpp:625
BSEManager::dumpRandConstant
void dumpRandConstant(const char *_fname)
dump rand constant to file
Definition: bse_interface.h:1103
BSEManager::getTimeStepStar
double getTimeStepStar(StarParameter &_star)
get next time step to check in Myr
Definition: bse_interface.h:1689
ARInteraction::calcInnerAccPotAndGTKickInv
Float calcInnerAccPotAndGTKickInv(AR::Force *_force, Float &_epot, const PtclHard *_particles, const int _n_particle)
calculate inner member acceleration, potential and inverse time transformation function gradient and ...
Definition: ar_interaction.hpp:171
StarParameter
SSE/BSE based code star parameter for saving.
Definition: bse_interface.h:257
BinaryEvent::getType
int getType(const int index) const
Get the binary type.
Definition: bse_interface.h:544