SlowDown Algorithmic Regularization (SDAR)
Algorithmic Regularization with slowdown method for integrating few-body motions
symplectic_integrator.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <functional>
4 #include "Common/list.h"
5 #include "Common/particle_group.h"
6 #include "AR/symplectic_step.h"
7 #include "AR/force.h"
8 #include "AR/slow_down.h"
9 #include "AR/profile.h"
10 #include "AR/information.h"
11 #include "AR/interrupt.h"
12 
14 
17 namespace AR {
18 
20  void printFeatures(std::ostream & fout) {
21 #ifdef AR_TTL
22  fout<<"Use AR TTL method\n";
23 #else
24  fout<<"Use AR LogH method\n";
25 #endif
26 #ifdef AR_SLOWDOWN_TREE
27  fout<<"Use slowdown Tree method\n";
28 #endif
29 #ifdef AR_SLOWDOWN_ARRAY
30  fout<<"Use slowdown array method\n";
31 #endif
32 #ifdef AR_SLOWDOWN_TIMESCALE
33  fout<<"Use slowdown timescale criterion\n";
34 #endif
35 #ifdef AR_SLOWDOWN_MASSRATIO
36  fout<<"Use slowdown mass ratio criterion\n";
37 #endif
38  }
39 
41  void printDebugFeatures(std::ostream & fout) {
42 #ifdef AR_DEBUG
43  fout<<"Debug mode: AR\n";
44 #endif
45  }
46 
48  void printReference(std::ostream & fout, const int offset=4) {
49  for (int i=0; i<offset; i++) fout<<" ";
50  fout<<"SDAR: Wang L., Nitadori K., Makino J., 2020, MNRAS, 493, 3398"
51  <<std::endl;
52  }
53 
55 
57  template <class Tmethod>
58  class TimeTransformedSymplecticManager {
59  public:
60  Float time_error_max;
61  Float energy_error_relative_max;
62  Float time_step_min;
63  Float ds_scale;
64  Float slowdown_pert_ratio_ref;
65 #ifdef AR_SLOWDOWN_MASSRATIO
66  Float slowdown_mass_ref;
67 #endif
68  Float slowdown_timescale_max;
69  long long unsigned int step_count_max;
70  int interrupt_detection_option;
71 
72  Tmethod interaction;
73  SymplecticStep step;
74 
76  TimeTransformedSymplecticManager(): time_error_max(Float(-1.0)), energy_error_relative_max(Float(-1.0)), time_step_min(Float(-1.0)), ds_scale(1.0), slowdown_pert_ratio_ref(Float(-1.0)),
77 #ifdef AR_SLOWDOWN_MASSRATIO
78  slowdown_mass_ref(Float(-1.0)),
79 #endif
80  slowdown_timescale_max(0.0),
81  step_count_max(0), interrupt_detection_option(0), interaction(), step() {}
82 
84 
86  bool checkParams() {
87  //ASSERT(time_error_max>ROUND_OFF_ERROR_LIMIT);
88  ASSERT(time_error_max>0.0);
89  ASSERT(energy_error_relative_max>ROUND_OFF_ERROR_LIMIT);
90  //ASSERT(time_step_min>ROUND_OFF_ERROR_LIMIT);
91  ASSERT(time_step_min>0.0);
92  ASSERT(ds_scale>0.0);
93  ASSERT(slowdown_pert_ratio_ref>0.0);
94 #ifdef AR_SLOWDOWN_MASSRATIO
95  ASSERT(slowdown_mass_ref>0.0);
96 #endif
97  ASSERT(slowdown_timescale_max>0);
98  ASSERT(step_count_max>0);
99  ASSERT(step.getOrder()>0);
100  ASSERT(interaction.checkParams());
101  return true;
102  }
103 
105 
107  void writeBinary(FILE *_fout) {
108  size_t size = sizeof(*this) - sizeof(interaction) - sizeof(step);
109  fwrite(this, size, 1,_fout);
110  interaction.writeBinary(_fout);
111  step.writeBinary(_fout);
112  }
113 
115 
118  void readBinary(FILE *_fin, int _version=0) {
119  if (_version==0) {
120  size_t size = sizeof(*this) - sizeof(interaction) - sizeof(step);
121  size_t rcount = fread(this, size, 1, _fin);
122  if (rcount<1) {
123  std::cerr<<"Error: TimeTransformedSymplecticManager parameter reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
124  abort();
125  }
126  }
127  else if (_version==1) {
128  size_t rcount = fread(this, sizeof(Float), 3, _fin);
129  if (rcount<3) {
130  std::cerr<<"Error: TimeTransformedSymplecticManager parameter data reading fails! requiring data number is 3, only obtain "<<rcount<<".\n";
131  abort();
132  }
133  ds_scale=1.0;
134  size_t size = sizeof(*this) - sizeof(interaction) - sizeof(step) - 4*sizeof(Float);
135  rcount = fread(&slowdown_pert_ratio_ref, size, 1, _fin);
136  if (rcount<1) {
137  std::cerr<<"Error: TimeTransformedSymplecticManager parameter data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
138  abort();
139  }
140  }
141  else {
142  std::cerr<<"Error: TimeTransformedSymplecticManager.readBinary unknown version "<<_version<<", should be 0 or 1."<<std::endl;
143  abort();
144  }
145  interaction.readBinary(_fin);
146  step.readBinary(_fin);
147  }
148 
150  void print(std::ostream & _fout) const{
151  _fout<<"time_error_max : "<<time_error_max<<std::endl
152  <<"energy_error_relative_max : "<<energy_error_relative_max<<std::endl
153  <<"time_step_min : "<<time_step_min<<std::endl
154  <<"slowdown_pert_ratio_ref : "<<slowdown_pert_ratio_ref<<std::endl
155 #ifdef AR_SLOWDOWN_MASSRATIO
156  <<"slowdown_mass_ref : "<<slowdown_mass_ref<<std::endl
157 #endif
158  <<"slowdown_timescale_max : "<<slowdown_timescale_max<<std::endl
159  <<"step_count_max : "<<step_count_max<<std::endl
160  <<"ds_scale : "<<ds_scale<<std::endl;
161  interaction.print(_fout);
162  step.print(_fout);
163  }
164  };
165 
167 
174  template <class Tparticle, class Tpcm, class Tpert, class Tmethod, class Tinfo>
175  class TimeTransformedSymplecticIntegrator {
176  private:
177  // intergrated variables
178  Float time_;
179  Float etot_ref_;
180 
181  // calculated varaiables
182  Float ekin_;
183  Float epot_;
184 
185  // cumulative slowdown (inner + outer) energy change
186  Float de_change_interrupt_; // energy change due to interruption
187  Float dH_change_interrupt_; // hamiltonian change due to interruption
188 
189 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
190  Float ekin_sd_;
191  Float epot_sd_;
192  Float etot_sd_ref_;
193 
194  Float de_sd_change_cum_; // slowdown energy change
195  Float dH_sd_change_cum_; // slowdown Hamiltonian change
196  Float de_sd_change_interrupt_; // slowdown energy change due to interruption
197  Float dH_sd_change_interrupt_; // slowdown energy change due to interruption
198 #endif
199 
200 #ifdef AR_TTL
201  // transformation factors
202  Float gt_drift_inv_;
203  Float gt_kick_inv_;
204 #endif
205 
206  // force array
207  COMM::List<Force> force_;
208 
209  public:
210  TimeTransformedSymplecticManager<Tmethod>* manager;
211  COMM::ParticleGroup<Tparticle,Tpcm> particles;
212 #ifdef AR_SLOWDOWN_ARRAY
213  COMM::List<AR::BinaryTree<Tparticle>*> binary_slowdown;
214 #endif
215  Tpert perturber;
216  Tinfo info;
217  Profile profile;
218 
220  TimeTransformedSymplecticIntegrator(): time_(0), etot_ref_(0), ekin_(0), epot_(0), de_change_interrupt_(0), dH_change_interrupt_(0),
221 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
222  ekin_sd_(0), epot_sd_(0), etot_sd_ref_(0),
223  de_sd_change_cum_(0), dH_sd_change_cum_(0), de_sd_change_interrupt_(0), dH_sd_change_interrupt_(0),
224 #endif
225 #ifdef AR_TTL
226  gt_drift_inv_(0), gt_kick_inv_(0),
227 #endif
228  force_(), manager(NULL), particles(),
229 #ifdef AR_SLOWDOWN_ARRAY
230  binary_slowdown(),
231 #endif
232  perturber(), info(), profile() {}
233 
235 
237  bool checkParams() {
238  ASSERT(manager!=NULL);
239  ASSERT(manager->checkParams());
240  ASSERT(perturber.checkParams());
241  ASSERT(info.checkParams());
242  return true;
243  }
244 
246 
248  void reserveIntegratorMem() {
249  // force array always allocated local memory
250  int nmax = particles.getSizeMax();
251  ASSERT(nmax>0);
252  force_.setMode(COMM::ListMode::local);
253  force_.reserveMem(nmax);
254 #ifdef AR_SLOWDOWN_ARRAY
255  binary_slowdown.setMode(COMM::ListMode::local);
256  binary_slowdown.reserveMem(nmax/2+1);
257 #endif
258  }
259 
261 
263  void clear() {
264  time_ = 0.0;
265  etot_ref_ =0.0;
266  ekin_ = 0.0;
267  epot_ = 0.0;
268  de_change_interrupt_ = 0.0;
269  dH_change_interrupt_ = 0.0;
270 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
271  ekin_sd_ = 0.0;
272  epot_sd_ = 0.0;
273  etot_sd_ref_ = 0.0;
274  de_sd_change_cum_ = 0.0;
275  dH_sd_change_cum_ = 0.0;
276  de_sd_change_interrupt_ = 0.0;
277  dH_sd_change_interrupt_ = 0.0;
278 #endif
279 #ifdef AR_TTL
280  gt_drift_inv_ = 0.0;
281  gt_kick_inv_ = 0.0;
282 #endif
283  force_.clear();
284  particles.clear();
285 #ifdef AR_SLOWDOWN_ARRAY
286  binary_slowdown.clear();
287 #endif
288  perturber.clear();
289  info.clear();
290  profile.clear();
291  }
292 
294  ~TimeTransformedSymplecticIntegrator() {
295  clear();
296  }
297 
299 
301  TimeTransformedSymplecticIntegrator& operator = (const TimeTransformedSymplecticIntegrator& _sym) {
302  clear();
303  time_ = _sym.time_;
304  etot_ref_ = _sym.etot_ref_;
305 
306  ekin_ = _sym.ekin_;
307  epot_ = _sym.epot_;
308  de_change_interrupt_= _sym.de_change_interrupt_;
309  dH_change_interrupt_= _sym.dH_change_interrupt_;
310 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
311  ekin_sd_= _sym.ekin_sd_;
312  epot_sd_= _sym.epot_sd_;
313  etot_sd_ref_= _sym.etot_sd_ref_;
314  de_sd_change_cum_= _sym.de_sd_change_cum_;
315  dH_sd_change_cum_= _sym.dH_sd_change_cum_;
316  de_sd_change_interrupt_= _sym.de_sd_change_interrupt_;
317  dH_sd_change_interrupt_= _sym.dH_sd_change_interrupt_;
318 #endif
319 #ifdef AR_TTL
320  gt_drift_inv_ = _sym.gt_drift_inv_;
321  gt_kick_inv_ = _sym.gt_kick_inv_;
322 #endif
323  force_ = _sym.force_;
324  manager = _sym.manager;
325 #ifdef AR_SLOWDOWN_ARRAY
326  binary_slowdown = _sym.binary_slowdown;
327 #endif
328  particles = _sym.particles;
329  info = _sym.binarytree;
330  profile = _sym.profile;
331 
332  return *this;
333  }
334 
335  private:
336 
337 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
338 
345  void calcSlowDownPertInnerBinaryIter(Float& _pert_out, Float& _t_min_sq, AR::BinaryTree<Tparticle>& _bini, AR::BinaryTree<Tparticle>& _binj) {
346  ASSERT(&_bini != &_binj);
347 // ASSERT(_bini.getMemberIndex(0)!=_binj.getMemberIndex(0));
348 
349  for (int k=0; k<2; k++) {
350  if (_binj.isMemberTree(k)) {
351  auto* bink = _binj.getMemberAsTree(k);
352  int check_flag = _bini.isSameBranch(*bink);
353  if (check_flag==0) {// no relation
354  if (bink->semi>0.0) manager->interaction.calcSlowDownPertOne(_pert_out, _t_min_sq, _bini, *bink);
355  else calcSlowDownPertInnerBinaryIter(_pert_out, _t_min_sq, _bini, *bink);
356  }
357  else if (check_flag==-2) { // _binj is the upper root tree
358  calcSlowDownPertInnerBinaryIter(_pert_out, _t_min_sq, _bini, *bink);
359  }
360  // other cases (same or sub branch), stop iteration.
361  }
362  else {
363  auto* pk = _binj.getMember(k);
364  if (pk->mass>0.0) manager->interaction.calcSlowDownPertOne(_pert_out, _t_min_sq, _bini, *pk);
365  }
366  }
367  }
368 
370 
373  void calcSlowDownInnerBinary(BinaryTree<Tparticle>& _bin) {
374  _bin.slowdown.pert_in = manager->interaction.calcPertFromBinary(_bin);
375 
376  Float pert_out = 0.0;
377  Float t_min_sq = NUMERIC_FLOAT_MAX;
378  auto& bin_root = info.getBinaryTreeRoot();
379 
380  calcSlowDownPertInnerBinaryIter(pert_out, t_min_sq, _bin, bin_root);
381 
382  _bin.slowdown.pert_out = pert_out + bin_root.slowdown.pert_out;
383 
384 #ifdef AR_SLOWDOWN_TIMESCALE
385  // velocity dependent method
386  //Float trv_ave = sdtdat.mtot/sqrt(sdtdat.mvor[0]*sdtdat.mvor[0] + sdtdat.mvor[1]*sdtdat.mvor[1] + sdtdat.mvor[2]*sdtdat.mvor[2]);
387  // get min of velocity and force dependent values
388  //Float t_min = std::min(trv_ave, sqrt(sdtdat.trf2_min));
389  _bin.slowdown.timescale = std::min(_bin.slowdown.getTimescaleMax(), sqrt(t_min_sq));
390 
391  // stablility criterion
392  // The slowdown factor should not make the system unstable, thus the Qst/Q set the limitation of the increasing of inner semi-major axis.
393  if (_bin.stab>0 && _bin.stab != NUMERIC_FLOAT_MAX) {
394  Float semi_amplify_max = std::max(Float(1.0),1.0/_bin.stab);
395  Float period_amplify_max = pow(semi_amplify_max,3.0/2.0);
396  Float timescale_stab = period_amplify_max*_bin.period;
397  _bin.slowdown.timescale = std::min(_bin.slowdown.timescale, timescale_stab);
398  }
399 
400 #else
401  _bin.slowdown.timescale = _bin.slowdown.getTimescaleMax();
402 #endif
403 
404  // only set slowdown if semi > 0 and stable
405  bool set_sd_flag = true;
406  if (_bin.semi>0) {
407  for (int k=0; k<2; k++) {
408  if (_bin.isMemberTree(k)) {
409  auto* bink = _bin.getMemberAsTree(k);
410  if (bink->stab>1.0) set_sd_flag = false;
411  }
412  }
413  }
414  else set_sd_flag = false;
415 
416  if (set_sd_flag) {
417  _bin.slowdown.period = _bin.period;
418  _bin.slowdown.calcSlowDownFactor();
419  }
420  else {
421  _bin.slowdown.setSlowDownFactor(1.0);
422  }
423  }
424 #endif // AR_SLOWDOWN_TREE || AR_SLOWDOWN_ARRAY
425 
426 #ifdef AR_SLOWDOWN_TREE
427 
429 
433  void calcTwoEKinIter(const Float& _inv_nest_sd_up, AR::BinaryTree<Tparticle>& _bin){
434  Float inv_nest_sd = _inv_nest_sd_up/_bin.slowdown.getSlowDownFactor();
435  Float* vel_cm = _bin.getVel();
436  for (int k=0; k<2; k++) {
437  Float* vk;
438  Float mk;
439  if (_bin.isMemberTree(k)) {
440  auto* bink = _bin.getMemberAsTree(k);
441  vk = bink->getVel();
442  mk = bink->mass;
443  calcTwoEKinIter(inv_nest_sd, *bink);
444  }
445  else {
446  auto* pk = _bin.getMember(k);
447  vk = pk->getVel();
448  mk = pk->mass;
449  ekin_ += mk * (vk[0]*vk[0]+vk[1]*vk[1]+vk[2]*vk[2]);
450  }
451  Float vrel[3] = {vk[0] - vel_cm[0],
452  vk[1] - vel_cm[1],
453  vk[2] - vel_cm[2]};
454  ekin_sd_ += mk * inv_nest_sd * (vrel[0]*vrel[0] + vrel[1]*vrel[1] + vrel[2]*vrel[2]);
455  }
456  }
457 
459  void calcEKin() {
460  ekin_ = ekin_sd_ = 0.0;
461  auto& bin_root=info.getBinaryTreeRoot();
462  Float sd_factor=1.0;
463 
464  calcTwoEKinIter(sd_factor, bin_root);
465  Float* vcm = bin_root.getVel();
466  // notice the cm velocity may not be zero after interruption, thus need to be added
467  ekin_sd_ += bin_root.mass*(vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2]);
468 
469  ekin_ *= 0.5;
470  ekin_sd_ *= 0.5;
471  }
472 
473 
475 
478  void kickVel(const Float& _dt) {
479  const int num = particles.getSize();
480  Tparticle* pdat = particles.getDataAddress();
481  Force* force = force_.getDataAddress();
482  for (int i=0; i<num; i++) {
483  // kick velocity
484  Float* vel = pdat[i].getVel();
485  Float* acc = force[i].acc_in;
486  Float* pert= force[i].acc_pert;
487  // half dv
488  vel[0] += _dt * (acc[0] + pert[0]);
489  vel[1] += _dt * (acc[1] + pert[1]);
490  vel[2] += _dt * (acc[2] + pert[2]);
491  }
492  // update binary c.m. velocity interation
493  updateBinaryVelIter(info.getBinaryTreeRoot());
494  }
495 
497 
503  void driftPosTreeIter(const Float& _dt, const Float* _vel_sd_up, const Float& _inv_nest_sd_up, AR::BinaryTree<Tparticle>& _bin) {
504  // current nested sd factor
505  Float inv_nest_sd = _inv_nest_sd_up/_bin.slowdown.getSlowDownFactor();
506 
507  Float* vel_cm = _bin.getVel();
508 
509  auto driftPos=[&](Float* pos, Float* vel, Float* vel_sd) {
510  //scale velocity referring to binary c.m.
511  vel_sd[0] = (vel[0] - vel_cm[0]) * inv_nest_sd + _vel_sd_up[0];
512  vel_sd[1] = (vel[1] - vel_cm[1]) * inv_nest_sd + _vel_sd_up[1];
513  vel_sd[2] = (vel[2] - vel_cm[2]) * inv_nest_sd + _vel_sd_up[2];
514 
515  pos[0] += _dt * vel_sd[0];
516  pos[1] += _dt * vel_sd[1];
517  pos[2] += _dt * vel_sd[2];
518  };
519 
520  for (int k=0; k<2; k++) {
521  if (_bin.isMemberTree(k)) {
522  auto* pj = _bin.getMemberAsTree(k);
523  Float* pos = pj->getPos();
524  Float* vel = pj->getVel();
525  Float vel_sd[3];
526  driftPos(pos, vel, vel_sd);
527  driftPosTreeIter(_dt, vel_sd, inv_nest_sd, *pj);
528 //#ifdef AR_DEBUG
529 // auto& pos1= pj->getLeftMember()->pos;
530 // auto& pos2= pj->getRightMember()->pos;
531 // Float m1 = pj->getLeftMember()->mass;
532 // Float m2 = pj->getRightMember()->mass;
533 // Float pos_cm[3] = {(m1*pos1[0]+m2*pos2[0])/(m1+m2),
534 // (m1*pos1[1]+m2*pos2[1])/(m1+m2),
535 // (m1*pos1[2]+m2*pos2[2])/(m1+m2)};
536 // Float dpos[3] = {pos[0]-pos_cm[0],
537 // pos[1]-pos_cm[1],
538 // pos[2]-pos_cm[2]};
539 // Float dr2 = dpos[0]*dpos[0]+dpos[1]*dpos[1]+dpos[2]*dpos[2];
540 // ASSERT(dr2<1e-10);
541 //#endif
542  }
543  else {
544  auto* pj = _bin.getMember(k);
545  Float* pos = pj->getPos();
546  Float* vel = pj->getVel();
547  Float vel_sd[3];
548  driftPos(pos, vel, vel_sd);
549  }
550  }
551  }
552 
554 
557  void driftTimeAndPos(const Float& _dt) {
558  // drift time
559  time_ += _dt;
560 
561  // the particle cm velocity is zero (assume in rest-frame)
562  ASSERT(!particles.isOriginFrame());
563  auto& bin_root=info.getBinaryTreeRoot();
564  Float vel_cm[3] = {0.0,0.0,0.0};
565  Float sd_factor=1.0;
566 
567  driftPosTreeIter(_dt, vel_cm, sd_factor, bin_root);
568  }
569 
570 #ifdef AR_TIME_FUNCTION_MULTI_R
571  Float calcInvR(Tparticle& _p1, Tparticle& _p2) {
573  Float dr[3] = {_p1.pos[0] - _p2.pos[0],
574  _p1.pos[1] - _p2.pos[1],
575  _p1.pos[2] - _p2.pos[2]};
576  Float r2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
577  Float invr = 1/sqrt(r2);
578  return invr;
579  }
580 
582  Float calcMultiInvRIter(AR::BinaryTree<Tparticle>& _bin){
583  Float gt_kick_inv= calcInvR(*_bin.getLeftMember(), *_bin.getRightMember());
584  for (int k=0; k<2; k++) {
585  if (_bini.isMemberTree(k)) // tree - tree
586  gt_kick_inv *= calcMultiInvRIter(*(_bini.getMemberAsTree(k)));
587  }
588  return gt_kick_inv;
589  }
590 
592  Float calcGtDriftInvIter(AR::BinaryTree<Tparticle>& _bin){
593  Float gt_drift_inv= calInvR(*_bin.getLeftMember(), *_bin.getRightMember());
594  for (int k=0; k<2; k++) {
595  if (_bini.isMemberTree(k)) // tree - tree
596  gt_kick_inv *= calcMultiInvRIter(*(_bini.getMemberAsTree(k)));
597  }
598  return gt_kick_inv;
599  }
600 #endif
601 
603 
609  Float calcAccPotAndGTKickInvTwo(const Float& _inv_nest_sd, const int _i, const int _j) {
610  ASSERT(_i>=0&&_i<particles.getSize());
611  ASSERT(_j>=0&&_j<particles.getSize());
612 
613  // calculate pair interaction
614  Force fij[2];
615  Float epotij;
616  Float gt_kick_inv = manager->interaction.calcInnerAccPotAndGTKickInvTwo(fij[0], fij[1], epotij, particles[_i], particles[_j]);
617 
618  // scale binary pair force with slowdown
619  force_[_i].acc_in[0] += fij[0].acc_in[0]*_inv_nest_sd;
620  force_[_i].acc_in[1] += fij[0].acc_in[1]*_inv_nest_sd;
621  force_[_i].acc_in[2] += fij[0].acc_in[2]*_inv_nest_sd;
622  force_[_j].acc_in[0] += fij[1].acc_in[0]*_inv_nest_sd;
623  force_[_j].acc_in[1] += fij[1].acc_in[1]*_inv_nest_sd;
624  force_[_j].acc_in[2] += fij[1].acc_in[2]*_inv_nest_sd;
625 
626  epot_ += epotij;
627  epot_sd_ += epotij*_inv_nest_sd;
628 
629 #ifdef AR_TTL
630  // scale gtgrad with slowdown
631  force_[_i].gtgrad[0] += fij[0].gtgrad[0]*_inv_nest_sd;
632  force_[_i].gtgrad[1] += fij[0].gtgrad[1]*_inv_nest_sd;
633  force_[_i].gtgrad[2] += fij[0].gtgrad[2]*_inv_nest_sd;
634  force_[_j].gtgrad[0] += fij[1].gtgrad[0]*_inv_nest_sd;
635  force_[_j].gtgrad[1] += fij[1].gtgrad[1]*_inv_nest_sd;
636  force_[_j].gtgrad[2] += fij[1].gtgrad[2]*_inv_nest_sd;
637 #endif
638 
639  return gt_kick_inv * _inv_nest_sd;
640  }
641 
643 
649  Float calcAccPotAndGTKickInvOneTreeIter(const Float& _inv_nest_sd, const int _i, AR::BinaryTree<Tparticle>& _bin) {
650  Float gt_kick_inv=0.0;
651  for (int k=0; k<2; k++) {
652  if (_bin.isMemberTree(k)) // particle - tree
653  gt_kick_inv += calcAccPotAndGTKickInvOneTreeIter(_inv_nest_sd, _i, *(_bin.getMemberAsTree(k)));
654  else // particle - particle
655  gt_kick_inv += calcAccPotAndGTKickInvTwo(_inv_nest_sd, _i, _bin.getMemberIndex(k));
656  }
657  return gt_kick_inv;
658  }
659 
661 
667  Float calcAccPotAndGTKickInvCrossTreeIter(const Float& _inv_nest_sd, AR::BinaryTree<Tparticle>& _bini, AR::BinaryTree<Tparticle>& _binj) {
668  ASSERT(&_bini!=&_binj);
669  Float gt_kick_inv=0.0;
670  for (int k=0; k<2; k++) {
671  if (_bini.isMemberTree(k)) // tree - tree
672  gt_kick_inv += calcAccPotAndGTKickInvCrossTreeIter(_inv_nest_sd, *(_bini.getMemberAsTree(k)), _binj);
673  else // particle - tree
674  gt_kick_inv += calcAccPotAndGTKickInvOneTreeIter(_inv_nest_sd, _bini.getMemberIndex(k), _binj);
675  }
676  return gt_kick_inv;
677  }
678 
680 
685  Float calcAccPotAndGTKickInvTreeIter(const Float& _inv_nest_sd_up, AR::BinaryTree<Tparticle>& _bin) {
686  // current nested sd factor
687  Float inv_nest_sd = _inv_nest_sd_up/_bin.slowdown.getSlowDownFactor();
688  Float gt_kick_inv = 0.0;
689 
690  // check left
691  if (_bin.isMemberTree(0)) { // left is tree
692  auto* bin_left = _bin.getMemberAsTree(0);
693  // inner interaction of left tree
694  gt_kick_inv += calcAccPotAndGTKickInvTreeIter(inv_nest_sd, *bin_left);
695 
696  if (_bin.isMemberTree(1)) { // right is tree
697  auto* bin_right = _bin.getMemberAsTree(1);
698 
699  // inner interaction of right tree
700  gt_kick_inv += calcAccPotAndGTKickInvTreeIter(inv_nest_sd, *bin_right);
701 
702  // cross interaction
703  gt_kick_inv += calcAccPotAndGTKickInvCrossTreeIter(inv_nest_sd, *bin_left, *bin_right);
704  }
705  else { // right is particle
706  // cross interaction from particle j to tree left
707  gt_kick_inv += calcAccPotAndGTKickInvOneTreeIter(inv_nest_sd, _bin.getMemberIndex(1), *bin_left);
708  }
709  }
710  else { // left is particle
711  if (_bin.isMemberTree(1)) { // right is tree
712  auto* bin_right = _bin.getMemberAsTree(1);
713  // inner interaction of right tree
714  gt_kick_inv += calcAccPotAndGTKickInvTreeIter(inv_nest_sd, *bin_right);
715 
716  // cross interaction from particle i to tree right
717  gt_kick_inv += calcAccPotAndGTKickInvOneTreeIter(inv_nest_sd, _bin.getMemberIndex(0), *bin_right);
718  }
719  else { // right is particle
720  // particle - particle interaction
721  gt_kick_inv += calcAccPotAndGTKickInvTwo(inv_nest_sd, _bin.getMemberIndex(0), _bin.getMemberIndex(1));
722  }
723  }
724 
725 
726 #ifdef AR_TIME_FUNCTION_MULTI_R
727  gt_kick_inv = calcMultiInvRIter(_bin);
728 #endif
729 
730  return gt_kick_inv;
731  }
732 
733 
735 
738  inline Float calcAccPotAndGTKickInv() {
739  epot_ = 0.0;
740  epot_sd_ = 0.0;
741  for (int i=0; i<force_.getSize(); i++) force_[i].clear();
742 
743  Float gt_kick_inv = calcAccPotAndGTKickInvTreeIter(1.0, info.getBinaryTreeRoot());
744  // pertuber force
745  manager->interaction.calcAccPert(force_.getDataAddress(), particles.getDataAddress(), particles.getSize(), particles.cm, perturber, getTime());
746  return gt_kick_inv;
747  }
748 
749 #ifdef AR_TTL
750 
758  Float kickEtotAndGTDriftTreeIter(const Float& _dt, const Float* _vel_sd_up, const Float& _inv_nest_sd_up, AR::BinaryTree<Tparticle>& _bin) {
759  // current nested sd factor
760  Float inv_nest_sd = _inv_nest_sd_up/_bin.slowdown.getSlowDownFactor();
761  Float dgt_drift_inv = 0.0;
762  Float de = 0.0;
763 
764  Float* vel_cm = _bin.getVel();
765 
766  for (int k=0; k<2; k++) {
767  if (_bin.isMemberTree(k)) {
768  auto* bink = _bin.getMemberAsTree(k);
769  Float* vel = bink->getVel();
770  Float vel_sd[3] = {(vel[0] - vel_cm[0]) * inv_nest_sd + _vel_sd_up[0],
771  (vel[1] - vel_cm[1]) * inv_nest_sd + _vel_sd_up[1],
772  (vel[2] - vel_cm[2]) * inv_nest_sd + _vel_sd_up[2]};
773  dgt_drift_inv += kickEtotAndGTDriftTreeIter(_dt, vel_sd, inv_nest_sd, *bink);
774  }
775  else {
776  int i = _bin.getMemberIndex(k);
777  ASSERT(i>=0&&i<particles.getSize());
778  ASSERT(&particles[i]==_bin.getMember(k));
779 
780  Float* gtgrad = force_[i].gtgrad;
781  Float* pert = force_[i].acc_pert;
782  Float* vel = particles[i].getVel();
783  Float vel_sd[3] = {(vel[0] - vel_cm[0]) * inv_nest_sd + _vel_sd_up[0],
784  (vel[1] - vel_cm[1]) * inv_nest_sd + _vel_sd_up[1],
785  (vel[2] - vel_cm[2]) * inv_nest_sd + _vel_sd_up[2]};
786 
787  de += particles[i].mass * (vel[0] * pert[0] +
788  vel[1] * pert[1] +
789  vel[2] * pert[2]);
790 #ifndef AR_TIME_FUNCTION_MULTI_R
791  dgt_drift_inv += (vel_sd[0] * gtgrad[0] +
792  vel_sd[1] * gtgrad[1] +
793  vel_sd[2] * gtgrad[2]);
794 #endif
795  }
796  }
797  etot_ref_ += _dt * de;
798  etot_sd_ref_ += _dt * de;
799 
800  return dgt_drift_inv;
801  }
802 
804 
807  void kickEtotAndGTDrift(const Float _dt) {
808  // the particle cm velocity is zero (assume in rest-frame)
809  ASSERT(!particles.isOriginFrame());
810  auto& bin_root=info.getBinaryTreeRoot();
811  Float vel_cm[3] = {0.0,0.0,0.0};
812  Float sd_factor=1.0;
813 
814  Float dgt_drift_inv = kickEtotAndGTDriftTreeIter(_dt, vel_cm, sd_factor, bin_root);
815 #ifdef AR_TIME_FUNCTION_MULTI_R
816  dgt_drift_inv = calcGTDriftInvIter(bin_root);
817 #endif
818  gt_drift_inv_ += dgt_drift_inv*_dt;
819  }
820 
821 #else
822 
826  inline void kickEtot(const Float _dt) {
827  Float de = 0.0;
828  const int num = particles.getSize();
829  Tparticle* pdat = particles.getDataAddress();
830  Force* force = force_.getDataAddress();
831  for (int i=0;i<num;i++) {
832  Float mass= pdat[i].mass;
833  Float* vel = pdat[i].getVel();
834  Float* pert= force[i].acc_pert;
835  de += mass * (vel[0] * pert[0] +
836  vel[1] * pert[1] +
837  vel[2] * pert[2]);
838  }
839  etot_sd_ref_ += _dt * de;
840  etot_ref_ += _dt * de;
841  }
842 #endif // AR_TTL
843 
844 #else
845 
846 #ifdef AR_SLOWDOWN_ARRAY
847 
851  inline void correctAccPotGTKickInvSlowDownInner(Float& _gt_kick_inv) {
852  int n = binary_slowdown.getSize();
853  Float gt_kick_inv_cor = 0.0;
854  Float de = 0.0;
855  for (int i=1; i<n; i++) {
856  auto& sdi = binary_slowdown[i];
857  ASSERT(sdi!=NULL);
858  int i1 = sdi->getMemberIndex(0);
859  int i2 = sdi->getMemberIndex(1);
860  Float kappa = sdi->slowdown.getSlowDownFactor();
861  Float kappa_inv = 1.0/kappa;
862  if (i1>=0) {
863  ASSERT(i2>=0);
864  ASSERT(i1!=i2);
865  ASSERT(i1<particles.getSize());
866  ASSERT(i2<particles.getSize());
867 
868  // calculate pair interaction
869  Force fi[2];
870  Float epoti;
871  Float gt_kick_inv_i = manager->interaction.calcInnerAccPotAndGTKickInvTwo(fi[0], fi[1], epoti, particles[i1], particles[i2]);
872  // scale binary pair force with slowdown
873  force_[i1].acc_in[0] += fi[0].acc_in[0]*kappa_inv - fi[0].acc_in[0];
874  force_[i1].acc_in[1] += fi[0].acc_in[1]*kappa_inv - fi[0].acc_in[1];
875  force_[i1].acc_in[2] += fi[0].acc_in[2]*kappa_inv - fi[0].acc_in[2];
876  force_[i2].acc_in[0] += fi[1].acc_in[0]*kappa_inv - fi[1].acc_in[0];
877  force_[i2].acc_in[1] += fi[1].acc_in[1]*kappa_inv - fi[1].acc_in[1];
878  force_[i2].acc_in[2] += fi[1].acc_in[2]*kappa_inv - fi[1].acc_in[2];
879 
880  de += epoti*kappa_inv - epoti;
881  // scale gtgrad with slowdown
882 #ifdef AR_TTL
883  force_[i1].gtgrad[0] += fi[0].gtgrad[0]*kappa_inv - fi[0].gtgrad[0];
884  force_[i1].gtgrad[1] += fi[0].gtgrad[1]*kappa_inv - fi[0].gtgrad[1];
885  force_[i1].gtgrad[2] += fi[0].gtgrad[2]*kappa_inv - fi[0].gtgrad[2];
886  force_[i2].gtgrad[0] += fi[1].gtgrad[0]*kappa_inv - fi[1].gtgrad[0];
887  force_[i2].gtgrad[1] += fi[1].gtgrad[1]*kappa_inv - fi[1].gtgrad[1];
888  force_[i2].gtgrad[2] += fi[1].gtgrad[2]*kappa_inv - fi[1].gtgrad[2];
889 #endif
890 
891  // gt kick
892  gt_kick_inv_cor += gt_kick_inv_i*(kappa_inv - 1.0);
893  }
894  }
895 
896  // global slowdown
897  const Float kappa_inv_global = 1.0/binary_slowdown[0]->slowdown.getSlowDownFactor();
898  for (int i=0; i<force_.getSize(); i++) {
899  force_[i].acc_in[0] *= kappa_inv_global;
900  force_[i].acc_in[1] *= kappa_inv_global;
901  force_[i].acc_in[2] *= kappa_inv_global;
902 #ifdef AR_TTL
903  force_[i].gtgrad[0] *= kappa_inv_global;
904  force_[i].gtgrad[1] *= kappa_inv_global;
905  force_[i].gtgrad[2] *= kappa_inv_global;
906 #endif
907  }
908 
909  epot_sd_ = (epot_ + de)*kappa_inv_global;
910  _gt_kick_inv = (_gt_kick_inv + gt_kick_inv_cor)*kappa_inv_global;
911  }
912 
914 
918  inline void correctPosSlowDownInner(const Float _dt, const Float _sd_global_inv) {
919  int n = binary_slowdown.getSize();
920  for (int i=1; i<n; i++) {
921  auto& sdi = binary_slowdown[i];
922  ASSERT(sdi!=NULL);
923  Float kappa = sdi->slowdown.getSlowDownFactor();
924  Float kappa_inv_m_one = (1.0/kappa - 1.0)*_sd_global_inv;
925  Float* velcm = sdi->getVel();
926  for (int k=0; k<2; k++) {
927  int j = sdi->getMemberIndex(k);
928  ASSERT(j>=0&&j<particles.getSize());
929  Float* pos = particles[j].getPos();
930  Float* vel = particles[j].getVel();
931 
932  // only scale velocity referring to binary c.m.
933  Float vrel[3] = { vel[0] - velcm[0],
934  vel[1] - velcm[1],
935  vel[2] - velcm[2]};
936  pos[0] += _dt * vrel[0] * kappa_inv_m_one;
937  pos[1] += _dt * vrel[1] * kappa_inv_m_one;
938  pos[2] += _dt * vrel[2] * kappa_inv_m_one;
939  }
940  }
941  }
942 
944  inline void updateCenterOfMassForBinaryWithSlowDownInner() {
945  int n = binary_slowdown.getSize();
946  for (int i=1; i<n; i++) {
947  auto& sdi = binary_slowdown[i];
948  int i1 = sdi->getMemberIndex(0);
949  int i2 = sdi->getMemberIndex(1);
950  ASSERT(i1>=0&&i1<particles.getSize());
951  ASSERT(i2>=0&&i2<particles.getSize());
952 
953  Float m1 = particles[i1].mass;
954  Float* pos1 = particles[i1].getPos();
955  Float* vel1 = particles[i1].getVel();
956  Float m2 = particles[i2].mass;
957  Float* pos2 = particles[i2].getPos();
958  Float* vel2 = particles[i2].getVel();
959  Float mcm = m1+m2;
960 
961  // first obtain the binary c.m. velocity
962  Float mcminv = 1.0/mcm;
963 
964  sdi->mass = mcm;
965  Float* pos = sdi->getPos();
966  pos[0] = (m1*pos1[0] + m2*pos2[0])*mcminv;
967  pos[1] = (m1*pos1[1] + m2*pos2[1])*mcminv;
968  pos[2] = (m1*pos1[2] + m2*pos2[2])*mcminv;
969 
970  Float* vel = sdi->getVel();
971  vel[0] = (m1*vel1[0] + m2*vel2[0])*mcminv;
972  vel[1] = (m1*vel1[1] + m2*vel2[1])*mcminv;
973  vel[2] = (m1*vel1[2] + m2*vel2[2])*mcminv;
974  }
975  }
976 
978 
980  inline void calcEkinSlowDownInner(const Float& _ekin) {
981  int n = binary_slowdown.getSize();
982  if (n==0) return;
983  const Float kappa_inv_global = 1.0/binary_slowdown[0]->slowdown.getSlowDownFactor();
984  Float de = Float(0.0);
985  for (int i=1; i<n; i++) {
986  auto& sdi = binary_slowdown[i];
987  ASSERT(sdi!=NULL);
988  Float kappa = sdi->slowdown.getSlowDownFactor();
989  Float kappa_inv_m_one = 1.0/kappa - 1.0;
990  Float* velcm = sdi->getVel();
991  for (int k=0; k<2; k++) {
992  int j = sdi->getMemberIndex(k);
993  ASSERT(j>=0&&j<particles.getSize());
994  Float* vel = particles[j].getVel();
995 
996  // only scale velocity referring to binary c.m.
997  Float vrel[3] = { vel[0] - velcm[0],
998  vel[1] - velcm[1],
999  vel[2] - velcm[2]};
1000 
1001  de += kappa_inv_m_one * particles[j].mass * (vrel[0]*vrel[0] + vrel[1]*vrel[1] + vrel[2]*vrel[2]);
1002  }
1003  }
1004  ekin_sd_ = (_ekin + 0.5*de)*kappa_inv_global;
1005  }
1006 
1008  static int findSlowDownInnerBinaryIter(COMM::List<AR::BinaryTree<Tparticle>*>& _binary_slowdown, const int& _c1, const int& _c2, AR::BinaryTree<Tparticle>& _bin) {
1009  // find leaf binary
1010  if (_bin.getMemberN()==2 && _bin.semi>0.0) {
1011  _binary_slowdown.increaseSizeNoInitialize(1);
1012  auto& sdi_new = _binary_slowdown.getLastMember();
1013  sdi_new = &_bin;
1014  return _c1+_c2+1;
1015  }
1016  else return _c1+_c2;
1017  }
1018 
1020 
1023  void findSlowDownInner(const Float _time) {
1024  auto& bin_root = info.getBinaryTreeRoot();
1025  binary_slowdown.resizeNoInitialize(1);
1026  int ncount[2]={0,0};
1027  int nsd = bin_root.processTreeIter(binary_slowdown, ncount[0], ncount[1], findSlowDownInnerBinaryIter);
1028  ASSERT(nsd==binary_slowdown.getSize()-1);
1029 
1030  for (int i=1; i<=nsd; i++) {
1031 #ifdef AR_SLOWDOWN_MASSRATIO
1032  const Float mass_ratio = manager->slowdown_mass_ref/binary_slowdown[i].bin->mass;
1033 #else
1034  const Float mass_ratio = 1.0;
1035 #endif
1036  binary_slowdown[i]->slowdown.initialSlowDownReference(mass_ratio*manager->slowdown_pert_ratio_ref, manager->slowdown_timescale_max);
1037  binary_slowdown[i]->slowdown.setUpdateTime(time_);
1038  }
1039  }
1040 
1041 #endif // AR_SLOWDOWN_ARRAY
1042 
1043 
1045  inline void calcEKin(){
1046  ekin_ = Float(0.0);
1047  const int num = particles.getSize();
1048  Tparticle* pdat = particles.getDataAddress();
1049  for (int i=0; i<num; i++) {
1050  const Float *vi=pdat[i].getVel();
1051  ekin_ += 0.5 * pdat[i].mass * (vi[0]*vi[0]+vi[1]*vi[1]+vi[2]*vi[2]);
1052  }
1053 #ifdef AR_SLOWDOWN_ARRAY
1054  calcEkinSlowDownInner(ekin_);
1055 #endif
1056  }
1057 
1059 
1062  inline void kickVel(const Float _dt) {
1063  const int num = particles.getSize();
1064  Tparticle* pdat = particles.getDataAddress();
1065  Force* force = force_.getDataAddress();
1066  for (int i=0; i<num; i++) {
1067  // kick velocity
1068  Float* vel = pdat[i].getVel();
1069  Float* acc = force[i].acc_in;
1070  Float* pert= force[i].acc_pert;
1071  // half dv
1072  vel[0] += _dt * (acc[0] + pert[0]);
1073  vel[1] += _dt * (acc[1] + pert[1]);
1074  vel[2] += _dt * (acc[2] + pert[2]);
1075  }
1076 #ifdef AR_SLOWDOWN_ARRAY
1077  // update c.m. of binaries
1078  updateCenterOfMassForBinaryWithSlowDownInner();
1079 #endif
1080  }
1081 
1083 
1086  inline void driftTimeAndPos(const Float _dt) {
1087  // drift time
1088  time_ += _dt;
1089 
1090  // drift position
1091  const int num = particles.getSize();
1092  Tparticle* pdat = particles.getDataAddress();
1093 #ifdef AR_SLOWDOWN_ARRAY
1094  const Float kappa_inv = 1.0/binary_slowdown[0]->slowdown.getSlowDownFactor();
1095  const Float dt_sd = _dt * kappa_inv;
1096 
1097  for (int i=0; i<num; i++) {
1098  Float* pos = pdat[i].getPos();
1099  Float* vel = pdat[i].getVel();
1100  pos[0] += dt_sd * vel[0];
1101  pos[1] += dt_sd * vel[1];
1102  pos[2] += dt_sd * vel[2];
1103  }
1104 
1105  // correct postion drift due to inner binary slowdown
1106  correctPosSlowDownInner(_dt, kappa_inv);
1107 #else
1108  for (int i=0; i<num; i++) {
1109  Float* pos = pdat[i].getPos();
1110  Float* vel = pdat[i].getVel();
1111  pos[0] += _dt * vel[0];
1112  pos[1] += _dt * vel[1];
1113  pos[2] += _dt * vel[2];
1114  }
1115 #endif
1116  }
1117 
1119 
1122  inline Float calcAccPotAndGTKickInv() {
1123  Float gt_kick_inv = manager->interaction.calcAccPotAndGTKickInv(force_.getDataAddress(), epot_, particles.getDataAddress(), particles.getSize(), particles.cm, perturber, getTime());
1124 
1125 //#ifdef AR_DEBUG
1126 // // check c.m. force
1127 // Force fcm;
1128 // Float mcm=0.0;
1129 // for (int i=0; i<particles.getSize(); i++) {
1130 // for (int k=0; k<3; k++) {
1131 // fcm.acc_in[k] += particles[i].mass * force_[i].acc_in[k];
1132 // }
1133 // mcm += particles[i].mass;
1134 // }
1135 // for (int k=0; k<3; k++) {
1136 // fcm.acc_in[k] /= mcm;
1137 // ASSERT(abs(fcm.acc_in[k])<1e-10);
1138 // }
1139 //#endif
1140 #ifdef AR_SLOWDOWN_ARRAY
1141  // slowdown binary acceleration
1142  correctAccPotGTKickInvSlowDownInner(gt_kick_inv);
1143 //#ifdef AR_DEBUG
1144 // // check c.m. force
1145 // fcm.acc_in[0] = fcm.acc_in[1] = fcm.acc_in[2] = 0.0;
1146 // for (int i=0; i<particles.getSize(); i++) {
1147 // for (int k=0; k<3; k++) {
1148 // fcm.acc_in[k] += particles[i].mass * force_[i].acc_in[k];
1149 // }
1150 // }
1151 // for (int k=0; k<3; k++) {
1152 // fcm.acc_in[k] /= mcm;
1153 // ASSERT(abs(fcm.acc_in[k])<1e-10);
1154 // }
1155 //#endif
1156 #endif
1157  return gt_kick_inv;
1158  }
1159 
1160 
1161 #ifdef AR_TTL
1162 
1167  inline void kickEtotAndGTDrift(const Float _dt) {
1168  Float de = Float(0.0);
1169  Float dg = Float(0.0);
1170  const int num = particles.getSize();
1171  Tparticle* pdat = particles.getDataAddress();
1172  Force* force = force_.getDataAddress();
1173  for (int i=0;i<num;i++) {
1174  Float mass= pdat[i].mass;
1175  Float* vel = pdat[i].getVel();
1176  Float* pert= force[i].acc_pert;
1177  Float* gtgrad=force[i].gtgrad;
1178  de += mass * (vel[0] * pert[0] +
1179  vel[1] * pert[1] +
1180  vel[2] * pert[2]);
1181  dg += (vel[0] * gtgrad[0] +
1182  vel[1] * gtgrad[1] +
1183  vel[2] * gtgrad[2]);
1184  }
1185  etot_ref_ += _dt * de;
1186 
1187 #ifdef AR_SLOWDOWN_ARRAY
1188  etot_sd_ref_ += _dt * de;
1189 
1190  // correct gt_drift_inv
1191  const Float kappa_inv_global = 1.0/binary_slowdown[0]->slowdown.getSlowDownFactor();
1192  int n = binary_slowdown.getSize();
1193  for (int i=1; i<n; i++) {
1194  auto& sdi = binary_slowdown[i];
1195  ASSERT(sdi!=NULL);
1196  Float kappa = sdi->slowdown.getSlowDownFactor();
1197  Float kappa_inv = 1.0/kappa;
1198  Float* velcm = sdi->getVel();
1199  for (int k=0; k<2; k++) {
1200  int j = sdi->getMemberIndex(k);
1201  if (j>=0) {
1202  ASSERT(j<particles.getSize());
1203  Float* gtgrad=force_[j].gtgrad;
1204  Float* vel = particles[j].getVel();
1205  Float vrel[3] = { vel[0] - velcm[0],
1206  vel[1] - velcm[1],
1207  vel[2] - velcm[2]};
1208  dg += (vrel[0] * (kappa_inv-1)* gtgrad[0] +
1209  vrel[1] * (kappa_inv-1)* gtgrad[1] +
1210  vrel[2] * (kappa_inv-1)* gtgrad[2]);
1211  }
1212  }
1213  }
1214  gt_drift_inv_ += _dt * dg *kappa_inv_global;
1215 #else
1216  gt_drift_inv_ += _dt * dg;
1217 #endif
1218  }
1219 
1220 #else
1221 
1225  inline void kickEtot(const Float _dt) {
1226  Float de = 0.0;
1227  const int num = particles.getSize();
1228  Tparticle* pdat = particles.getDataAddress();
1229  Force* force = force_.getDataAddress();
1230  for (int i=0;i<num;i++) {
1231  Float mass= pdat[i].mass;
1232  Float* vel = pdat[i].getVel();
1233  Float* pert= force[i].acc_pert;
1234  de += mass * (vel[0] * pert[0] +
1235  vel[1] * pert[1] +
1236  vel[2] * pert[2]);
1237  }
1238 #ifdef AR_SLOWDOWN_ARRAY
1239  etot_sd_ref_ += _dt * de;
1240 #endif
1241  etot_ref_ += _dt * de;
1242  }
1243 #endif //AR_TTL
1244 
1245 #endif // AR_SLOWDOWN_TREE
1246 
1248  void updateBinaryVelIter(AR::BinaryTree<Tparticle>& _bin) {
1249  _bin.vel[0]= _bin.vel[1] = _bin.vel[2] = 0.0;
1250  Float mcm_inv = 1.0/_bin.mass;
1251  for (int k=0; k<2; k++) {
1252  if (_bin.isMemberTree(k)) {
1253  auto* bink = _bin.getMemberAsTree(k);
1254  updateBinaryVelIter(*bink);
1255  _bin.vel[0] += bink->mass*bink->vel[0];
1256  _bin.vel[1] += bink->mass*bink->vel[1];
1257  _bin.vel[2] += bink->mass*bink->vel[2];
1258  }
1259  else {
1260  auto* pk = _bin.getMember(k);
1261  _bin.vel[0] += pk->mass*pk->vel[0];
1262  _bin.vel[1] += pk->mass*pk->vel[1];
1263  _bin.vel[2] += pk->mass*pk->vel[2];
1264  }
1265  }
1266  _bin.vel[0] *= mcm_inv;
1267  _bin.vel[1] *= mcm_inv;
1268  _bin.vel[2] *= mcm_inv;
1269  }
1270 
1272  void updateBinaryCMIter(AR::BinaryTree<Tparticle>& _bin) {
1273  _bin.pos[0]= _bin.pos[1] = _bin.pos[2] = 0.0;
1274  _bin.vel[0]= _bin.vel[1] = _bin.vel[2] = 0.0;
1275  Float mcm_member = 0.0;
1276  for (int k=0; k<2; k++) {
1277  if (_bin.isMemberTree(k)) {
1278  auto* bink = _bin.getMemberAsTree(k);
1279  updateBinaryCMIter(*bink);
1280  mcm_member += bink->mass;
1281  _bin.pos[0] += bink->mass*bink->pos[0];
1282  _bin.pos[1] += bink->mass*bink->pos[1];
1283  _bin.pos[2] += bink->mass*bink->pos[2];
1284  _bin.vel[0] += bink->mass*bink->vel[0];
1285  _bin.vel[1] += bink->mass*bink->vel[1];
1286  _bin.vel[2] += bink->mass*bink->vel[2];
1287  }
1288  else {
1289  auto* pk = _bin.getMember(k);
1290  mcm_member += pk->mass;
1291  _bin.pos[0] += pk->mass*pk->pos[0];
1292  _bin.pos[1] += pk->mass*pk->pos[1];
1293  _bin.pos[2] += pk->mass*pk->pos[2];
1294  _bin.vel[0] += pk->mass*pk->vel[0];
1295  _bin.vel[1] += pk->mass*pk->vel[1];
1296  _bin.vel[2] += pk->mass*pk->vel[2];
1297  }
1298  }
1299  Float mcm_inv = 1.0/mcm_member;
1300  _bin.mass = mcm_member;
1301  _bin.m1 = _bin.getLeftMember()->mass;
1302  _bin.m2 = _bin.getRightMember()->mass;
1303  _bin.vel[0] *= mcm_inv;
1304  _bin.vel[1] *= mcm_inv;
1305  _bin.vel[2] *= mcm_inv;
1306  _bin.pos[0] *= mcm_inv;
1307  _bin.pos[1] *= mcm_inv;
1308  _bin.pos[2] *= mcm_inv;
1309  }
1310 
1312  void setBinaryCMZeroIter(AR::BinaryTree<Tparticle>& _bin) {
1313  _bin.mass = 0.0;
1314  for (int k=0; k<2; k++) {
1315  if (_bin.isMemberTree(k)) {
1316  auto* bink = _bin.getMemberAsTree(k);
1317  setBinaryCMZeroIter(*bink);
1318  }
1319  }
1320  }
1321 
1323  bool updateBinarySemiEccPeriodIter(AR::BinaryTree<Tparticle>& _bin, const Float& _G, const Float _time, const bool _check=false) {
1324  bool check = _check;
1325  for (int k=0; k<2; k++) {
1326  if (_bin.isMemberTree(k)) {
1327  auto* bink = _bin.getMemberAsTree(k);
1328  check=updateBinarySemiEccPeriodIter(*bink, _G, _time, _check);
1329  }
1330  }
1331  if ((_time>_bin.stab_check_time&&_bin.stab>1.0&&_bin.m1>0.0&&_bin.m2>0.0)||check) {
1332  _bin.calcSemiEccPeriod(_G);
1333  _bin.stab_check_time = _time + _bin.period;
1334  return true;
1335  }
1336  return false;
1337  }
1338 
1339  public:
1340 #ifdef AR_SLOWDOWN_TREE
1341 
1346  void updateSlowDownAndCorrectEnergy(const bool _update_energy_flag, const bool _stable_check_flag) {
1347  auto& bin_root = info.getBinaryTreeRoot();
1348  auto& sd_root = bin_root.slowdown;
1349 
1350 #ifdef AR_TTL
1351  Float sd_backup = sd_root.getSlowDownFactor();
1352 #endif
1353 
1354  // when the maximum inner slowdown is large, the outer should not be slowed down since the system may not be stable.
1355  //if (inner_sd_change_flag&&sd_org_inner_max<1000.0*manager->slowdown_pert_ratio_ref) sd_root.setSlowDownFactor(1.0);
1356  //if (time_>=sd_root.getUpdateTime()) {
1357  sd_root.pert_in = manager->interaction.calcPertFromBinary(bin_root);
1358  sd_root.pert_out = 0.0;
1359  Float t_min_sq= NUMERIC_FLOAT_MAX;
1360  manager->interaction.calcSlowDownPert(sd_root.pert_out, t_min_sq, getTime(), particles.cm, perturber);
1361  sd_root.timescale = std::min(sd_root.getTimescaleMax(), sqrt(t_min_sq));
1362 
1363  //Float period_amplify_max = NUMERIC_FLOAT_MAX;
1364  if (_stable_check_flag) {
1365  // check whether the system is stable for 10000 out period and the apo-center is below break criterion
1366  Float stab = bin_root.stableCheckIter(bin_root,10000*bin_root.period);
1367  Float apo = bin_root.semi*(1+bin_root.ecc);
1368  if (stab<1.0 && apo<info.r_break_crit) {
1369  sd_root.period = bin_root.period;
1370  sd_root.calcSlowDownFactor();
1371  }
1372  else sd_root.setSlowDownFactor(1.0);
1373 
1374  // stablility criterion
1375  // The slowdown factor should not make the system unstable, thus the Qst/Q set the limitation of the increasing of inner semi-major axis.
1376  //if (stab>0 && stab != NUMERIC_FLOAT_MAX) {
1377  // Float semi_amplify_max = std::max(Float(1.0),1.0/stab);
1378  // period_amplify_max = pow(semi_amplify_max,3.0/2.0);
1379  //}
1380  }
1381  else if (bin_root.semi>0) {
1382  sd_root.period = bin_root.period;
1383  sd_root.calcSlowDownFactor();
1384  }
1385  else sd_root.setSlowDownFactor(1.0);
1386 
1387  //sd_root.increaseUpdateTimeOnePeriod();
1388  //}
1389 
1390  // inner binary slowdown
1391  Float sd_org_inner_max = 0.0;
1392  bool inner_sd_change_flag=false;
1393  int n_bin = info.binarytree.getSize();
1394  for (int i=0; i<n_bin-1; i++) {
1395  auto& bini = info.binarytree[i];
1396  //if (time_>=bini.slowdown.getUpdateTime()) {
1397  bini.calcCenterOfMass();
1398  calcSlowDownInnerBinary(bini);
1399 
1400  //sdi->slowdown.increaseUpdateTimeOnePeriod();
1401  sd_org_inner_max = std::max(bini.slowdown.getSlowDownFactorOrigin(),sd_org_inner_max);
1402  inner_sd_change_flag=true;
1403  //}
1404  }
1405 
1406 
1407  if (_update_energy_flag) {
1408  Float ekin_sd_bk = ekin_sd_;
1409  Float epot_sd_bk = epot_sd_;
1410  Float H_sd_bk = getHSlowDown();
1411  if(inner_sd_change_flag) {
1412 #ifdef AR_TTL
1413  Float gt_kick_inv_new = calcAccPotAndGTKickInv();
1414  gt_drift_inv_ += gt_kick_inv_new - gt_kick_inv_;
1415  gt_kick_inv_ = gt_kick_inv_new;
1416 #else
1417  calcAccPotAndGTKickInv();
1418 #endif
1419  calcEKin();
1420  }
1421  else {
1422  Float kappa_inv = 1.0/sd_root.getSlowDownFactor();
1423 #ifdef AR_TTL
1424  Float gt_kick_inv_new = gt_kick_inv_*sd_backup*kappa_inv;
1425  gt_drift_inv_ += gt_kick_inv_new - gt_kick_inv_;
1426  gt_kick_inv_ = gt_kick_inv_new;
1427 #endif
1428  ekin_sd_ = ekin_*kappa_inv;
1429  epot_sd_ = epot_*kappa_inv;
1430  }
1431  Float de_sd = (ekin_sd_ - ekin_sd_bk) + (epot_sd_ - epot_sd_bk);
1432  etot_sd_ref_ += de_sd;
1433 
1434  Float dH_sd = getHSlowDown() - H_sd_bk;
1435 
1436  // add slowdown change to the global slowdown energy
1437  de_sd_change_cum_ += de_sd;
1438  dH_sd_change_cum_ += dH_sd;
1439  }
1440  }
1441 #elif AR_SLOWDOWN_ARRAY
1442 
1447  void updateSlowDownAndCorrectEnergy(const bool _update_energy_flag, const bool _stable_check_flag) {
1448  auto& bin_root = *binary_slowdown[0];
1449  auto& sd_root = bin_root.slowdown;
1450 #ifdef AR_TTL
1451  Float sd_backup = sd_root.getSlowDownFactor();
1452 #endif
1453 
1454  // when the maximum inner slowdown is large, the outer should not be slowed down since the system may not be stable.
1455  //if (inner_sd_change_flag&&sd_org_inner_max<1000.0*manager->slowdown_pert_ratio_ref) sd_root.setSlowDownFactor(1.0);
1456  //if (time_>=sd_root.getUpdateTime()) {
1457  sd_root.pert_in = manager->interaction.calcPertFromBinary(bin_root);
1458  sd_root.pert_out = 0.0;
1459  Float t_min_sq= NUMERIC_FLOAT_MAX;
1460  manager->interaction.calcSlowDownPert(sd_root.pert_out, t_min_sq, getTime(), particles.cm, perturber);
1461  sd_root.timescale = std::min(sd_root.getTimescaleMax(), sqrt(t_min_sq));
1462 
1463  //Float period_amplify_max = NUMERIC_FLOAT_MAX;
1464  if (_stable_check_flag) {
1465  // check whether the system is stable for 10000 out period
1466  Float stab = bin_root.stableCheckIter(bin_root,10000*bin_root.period);
1467  if (stab<1.0) {
1468  sd_root.period = bin_root.period;
1469  sd_root.calcSlowDownFactor();
1470  }
1471  else sd_root.setSlowDownFactor(1.0);
1472  // stablility criterion
1473  // The slowdown factor should not make the system unstable, thus the Qst/Q set the limitation of the increasing of inner semi-major axis.
1474  //if (stab>0) {
1475  // Float semi_amplify_max = std::max(Float(1.0),1.0/stab);
1476  // period_amplify_max = pow(semi_amplify_max,3.0/2.0);
1477  //}
1478  }
1479  else if (bin_root.semi>0) {
1480  sd_root.period = bin_root.period;
1481  sd_root.calcSlowDownFactor();
1482  }
1483  else sd_root.setSlowDownFactor(1.0);
1484 
1485  //sd_root.increaseUpdateTimeOnePeriod();
1486  //}
1487 
1488 
1489  // inner binary slowdown
1490  int n = binary_slowdown.getSize();
1491  bool modified_flag=false;
1492  for (int i=1; i<n; i++) {
1493  auto* sdi = binary_slowdown[i];
1494  //if (time_>=sdi->slowdown.getUpdateTime()) {
1495  sdi->calcCenterOfMass();
1496  calcSlowDownInnerBinary(*sdi);
1497  //sdi->slowdown.increaseUpdateTimeOnePeriod();
1498  modified_flag=true;
1499  //}
1500  }
1501 
1502  if (_update_energy_flag) {
1503  Float ekin_sd_bk = ekin_sd_;
1504  Float epot_sd_bk = epot_sd_;
1505  Float H_sd_bk = getHSlowDown();
1506 
1507  if (modified_flag) {
1508  // initialize the gt_drift_inv_ with new slowdown factor
1509  Float gt_kick_inv_sdi = manager->interaction.calcAccPotAndGTKickInv(force_.getDataAddress(), epot_, particles.getDataAddress(), particles.getSize(), particles.cm, perturber, getTime());
1510  correctAccPotGTKickInvSlowDownInner(gt_kick_inv_sdi);
1511 #ifdef AR_TTL
1512  gt_drift_inv_ += gt_kick_inv_sdi - gt_kick_inv_;
1513  gt_kick_inv_ = gt_kick_inv_sdi;
1514 #endif
1515  // correct etot_sd_ref_ with new slowdown
1516  calcEkinSlowDownInner(ekin_);
1517  }
1518  else {
1519  Float kappa_inv = 1.0/sd_root.getSlowDownFactor();
1520 #ifdef AR_TTL
1521  Float gt_kick_inv_sdi = gt_kick_inv_*sd_backup*kappa_inv;
1522  gt_drift_inv_ += gt_kick_inv_sdi - gt_kick_inv_;
1523  gt_kick_inv_ = gt_kick_inv_sdi;
1524 #endif
1525  // only need to correct the total value
1526  ekin_sd_ = ekin_*kappa_inv;
1527  epot_sd_ = epot_*kappa_inv;
1528  }
1529 
1530  Float de_sd = (ekin_sd_ - ekin_sd_bk) + (epot_sd_ - epot_sd_bk);
1531  etot_sd_ref_ += de_sd;
1532 
1533  Float dH_sd = getHSlowDown() - H_sd_bk;
1534 
1535  // add slowdown change to the global slowdown energy
1536  de_sd_change_cum_ += de_sd;
1537  dH_sd_change_cum_ += dH_sd;
1538  }
1539  }
1540 #endif
1541 
1543 
1546  void initialIntegration(const Float _time) {
1547  ASSERT(checkParams());
1548 
1549  // particle number and data address
1550  const int n_particle = particles.getSize();
1551 
1552  // resize force array
1553  force_.resizeNoInitialize(n_particle);
1554 
1555  // Initial intgrt value t (avoid confusion of real time when slowdown is used)
1556  time_ = _time;
1557 
1558  // check particle number
1559  ASSERT(particles.getSize()>=2);
1560 
1561  // reset particle modification flag
1562  particles.setModifiedFalse();
1563 
1564  // check the center-of-mass initialization
1565  if(particles.isOriginFrame()) {
1566  particles.calcCenterOfMass();
1567  particles.shiftToCenterOfMassFrame();
1568  for (int i=0; i<info.binarytree.getSize(); i++) {
1569  auto& bin = info.binarytree[i];
1570  bin.pos[0] -= particles.cm.pos[0];
1571  bin.pos[1] -= particles.cm.pos[1];
1572  bin.pos[2] -= particles.cm.pos[2];
1573  bin.vel[0] -= particles.cm.vel[0];
1574  bin.vel[1] -= particles.cm.vel[1];
1575  bin.vel[2] -= particles.cm.vel[2];
1576  }
1577  }
1578  ASSERT(info.getBinaryTreeRoot().pos[0]*info.getBinaryTreeRoot().pos[0]<1e-10);
1579  ASSERT(info.getBinaryTreeRoot().vel[0]*info.getBinaryTreeRoot().vel[0]<1e-10);
1580 
1581 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
1582 #ifdef AR_SLOWDOWN_TREE
1583  for (int i=0; i<info.binarytree.getSize(); i++)
1584  info.binarytree[i].slowdown.initialSlowDownReference(manager->slowdown_pert_ratio_ref, manager->slowdown_timescale_max);
1585 #else
1586  binary_slowdown.increaseSizeNoInitialize(1);
1587  binary_slowdown[0] = &info.getBinaryTreeRoot();
1588 
1589  // set slowdown reference
1590  SlowDown& slowdown_root = info.getBinaryTreeRoot().slowdown;
1591 
1592  // slowdown for the system
1593  slowdown_root.initialSlowDownReference(manager->slowdown_pert_ratio_ref, manager->slowdown_timescale_max);
1594 
1595  if (particles.getSize()>2) {
1596  findSlowDownInner(time_);
1597  // update c.m. of binaries
1598  //updateCenterOfMassForBinaryWithSlowDownInner();
1599  }
1600 #endif
1601 
1602  updateSlowDownAndCorrectEnergy(false,true);
1603 
1604 #ifdef AR_TTL
1605  gt_kick_inv_ = calcAccPotAndGTKickInv();
1606 
1607  // initially gt_drift
1608  gt_drift_inv_ = gt_kick_inv_;
1609 
1610 #else
1611  calcAccPotAndGTKickInv();
1612 #endif
1613 
1614  calcEKin();
1615 
1616  etot_ref_ = ekin_ + epot_;
1617  etot_sd_ref_ = ekin_sd_ + epot_sd_;
1618 
1619  Float de_sd = etot_sd_ref_ - etot_ref_;
1620 
1621  // add slowdown change to the global slowdown energy
1622  de_sd_change_cum_ += de_sd;
1623  dH_sd_change_cum_ = 0.0;
1624 
1625 #else // No SLOWDOWN
1626  Tparticle* particle_data = particles.getDataAddress();
1627  Force* force_data = force_.getDataAddress();
1628 
1629 #ifdef AR_TTL
1630  gt_kick_inv_ = manager->interaction.calcAccPotAndGTKickInv(force_data, epot_, particle_data, n_particle, particles.cm, perturber, _time);
1631 
1632  // initially gt_drift
1633  gt_drift_inv_ = gt_kick_inv_;
1634 
1635 #else
1636  manager->interaction.calcAccPotAndGTKickInv(force_data, epot_, particle_data, n_particle, particles.cm, perturber, _time);
1637 #endif
1638 
1639  // calculate kinetic energy
1640  calcEKin();
1641 
1642  // initial total energy
1643  etot_ref_ = ekin_ + epot_;
1644 
1645 #endif
1646  }
1647 
1649 
1653  void integrateOneStep(const Float _ds, Float _time_table[]) {
1654  ASSERT(checkParams());
1655 
1656  ASSERT(!particles.isModified());
1657  ASSERT(_ds>0);
1658 
1659  // symplectic step coefficent group n_particleber
1660  const int nloop = manager->step.getCDPairSize();
1661 
1662  for (int i=0; i<nloop; i++) {
1663  // step for drift
1664  Float ds_drift = manager->step.getCK(i)*_ds;
1665 
1666  // inverse time transformation factor for drift
1667 #ifdef AR_TTL
1668  Float gt_drift_inv = gt_drift_inv_;
1669 #else
1670 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
1671  Float gt_drift_inv = manager->interaction.calcGTDriftInv(ekin_sd_-etot_sd_ref_); // pt = -etot
1672 #else
1673  Float gt_drift_inv = manager->interaction.calcGTDriftInv(ekin_-etot_ref_); // pt = -etot
1674 #endif
1675 #endif
1676 
1677  // drift
1678  Float dt_drift = ds_drift/gt_drift_inv;
1679 
1680  // drift time and postion
1681  driftTimeAndPos(dt_drift);
1682  _time_table[i] = time_;
1683 
1684  // step for kick
1685  Float ds_kick = manager->step.getDK(i)*_ds;
1686 
1688  Float gt_kick_inv = calcAccPotAndGTKickInv();
1689 
1690  // time step for kick
1691  Float dt_kick = ds_kick/gt_kick_inv;
1692 
1693  // kick half step for velocity
1694  kickVel(0.5*dt_kick);
1695 
1696 #ifdef AR_TTL
1697  // back up gt_kick
1698  gt_kick_inv_ = gt_kick_inv;
1699  // kick total energy and inverse time transformation factor for drift
1700  kickEtotAndGTDrift(dt_kick);
1701 #else
1702  // kick total energy
1703  kickEtot(dt_kick);
1704 #endif
1705  // kick half step for velocity
1706  kickVel(0.5*dt_kick);
1707 
1708  // calculate kinetic energy
1709  calcEKin();
1710  }
1711  }
1712 
1713 
1715 
1720  void integrateTwoOneStep(const Float _ds, Float _time_table[]) {
1721  ASSERT(checkParams());
1722 
1723  ASSERT(!particles.isModified());
1724  ASSERT(_ds>0);
1725 
1726  // symplectic step coefficent group number
1727  const int nloop = manager->step.getCDPairSize();
1728 
1729  const int n_particle = particles.getSize();
1730  ASSERT(n_particle==2);
1731 
1732 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
1733  const Float kappa_inv = 1.0/info.getBinaryTreeRoot().slowdown.getSlowDownFactor();
1734 #endif
1735 
1736  Tparticle* particle_data = particles.getDataAddress();
1737  Float mass1 = particle_data[0].mass;
1738  Float* pos1 = particle_data[0].getPos();
1739  Float* vel1 = particle_data[0].getVel();
1740 
1741  Float mass2 = particle_data[1].mass;
1742  Float* pos2 = particle_data[1].getPos();
1743  Float* vel2 = particle_data[1].getVel();
1744 
1745  Force* force_data = force_.getDataAddress();
1746  Float* acc1 = force_data[0].acc_in;
1747  Float* pert1= force_data[0].acc_pert;
1748 
1749  Float* acc2 = force_data[1].acc_in;
1750  Float* pert2= force_data[1].acc_pert;
1751 #ifdef AR_TTL
1752  Float* gtgrad1 = force_data[0].gtgrad;
1753  Float* gtgrad2 = force_data[1].gtgrad;
1754 #endif
1755 
1756 #ifdef AR_DEBUG_PRINT_DKD
1757  std::cout<<"K "<<time_<<" "
1758  <<pos2[0]-pos1[0]<<" "<<pos2[1]-pos1[1]<<" "<<pos2[2]-pos1[2]<<" "
1759  <<vel2[0]-vel1[0]<<" "<<vel2[1]-vel1[1]<<" "<<vel2[2]-vel1[2]<<" "
1760  <<ekin_<<" "<<epot_<<" "<<etot_ref_<<std::endl;
1761 #endif
1762 
1763  for (int i=0; i<nloop; i++) {
1764  // step for drift
1765  Float ds = manager->step.getCK(i)*_ds;
1766  // inverse time transformation factor for drift
1767 #ifdef AR_TTL
1768  Float gt_inv = gt_drift_inv_;
1769 #else
1770 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
1771  Float gt_inv = manager->interaction.calcGTDriftInv(ekin_sd_-etot_sd_ref_); // pt = -etot_sd
1772 #else
1773  Float gt_inv = manager->interaction.calcGTDriftInv(ekin_-etot_ref_); // pt = -etot
1774 #endif
1775 #endif
1776  // drift
1777  Float dt = ds/gt_inv;
1778  ASSERT(!ISNAN(dt));
1779 
1780  // drift time
1781  time_ += dt;
1782 
1783  // update real time
1784  _time_table[i] = time_;
1785 
1786 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
1787  Float dt_sd = dt*kappa_inv;
1788 
1789  // drift position
1790  pos1[0] += dt_sd * vel1[0];
1791  pos1[1] += dt_sd * vel1[1];
1792  pos1[2] += dt_sd * vel1[2];
1793 
1794  pos2[0] += dt_sd * vel2[0];
1795  pos2[1] += dt_sd * vel2[1];
1796  pos2[2] += dt_sd * vel2[2];
1797 #else
1798  // drift position
1799  pos1[0] += dt * vel1[0];
1800  pos1[1] += dt * vel1[1];
1801  pos1[2] += dt * vel1[2];
1802 
1803  pos2[0] += dt * vel2[0];
1804  pos2[1] += dt * vel2[1];
1805  pos2[2] += dt * vel2[2];
1806 #endif
1807 
1808  // step for kick
1809  ds = manager->step.getDK(i)*_ds;
1810 
1811  gt_inv = manager->interaction.calcAccPotAndGTKickInv(force_data, epot_, particle_data, n_particle, particles.cm, perturber, _time_table[i]);
1812 
1813  ASSERT(!ISNAN(epot_));
1814 
1815 #ifdef AR_DEBUG_PRINT_DKD
1816  if (i>0)
1817  std::cout<<"K "<<time_<<" "
1818  <<pos2[0]-pos1[0]<<" "<<pos2[1]-pos1[1]<<" "<<pos2[2]-pos1[2]<<" "
1819  <<vel2[0]-vel1[0]<<" "<<vel2[1]-vel1[1]<<" "<<vel2[2]-vel1[2]<<" "
1820  <<ekin_<<" "<<epot_<<" "<<etot_ref_<<std::endl;
1821 #endif
1822 
1823  // kick half step for velocity
1824  Float dvel1[3], dvel2[3];
1825 
1826 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
1827  // time step for kick
1828  gt_inv *= kappa_inv;
1829 
1830  dt = 0.5*ds/gt_inv;
1831 
1832  dvel1[0] = dt * (acc1[0]*kappa_inv + pert1[0]);
1833  dvel1[1] = dt * (acc1[1]*kappa_inv + pert1[1]);
1834  dvel1[2] = dt * (acc1[2]*kappa_inv + pert1[2]);
1835 
1836  dvel2[0] = dt * (acc2[0]*kappa_inv + pert2[0]);
1837  dvel2[1] = dt * (acc2[1]*kappa_inv + pert2[1]);
1838  dvel2[2] = dt * (acc2[2]*kappa_inv + pert2[2]);
1839 #else
1840  dt = 0.5*ds/gt_inv;
1841 
1842  dvel1[0] = dt * (acc1[0] + pert1[0]);
1843  dvel1[1] = dt * (acc1[1] + pert1[1]);
1844  dvel1[2] = dt * (acc1[2] + pert1[2]);
1845 
1846  dvel2[0] = dt * (acc2[0] + pert2[0]);
1847  dvel2[1] = dt * (acc2[1] + pert2[1]);
1848  dvel2[2] = dt * (acc2[2] + pert2[2]);
1849 #endif
1850 
1851  vel1[0] += dvel1[0];
1852  vel1[1] += dvel1[1];
1853  vel1[2] += dvel1[2];
1854 
1855  vel2[0] += dvel2[0];
1856  vel2[1] += dvel2[1];
1857  vel2[2] += dvel2[2];
1858 
1859 #ifdef AR_DEBUG_PRINT_DKD
1860  std::cout<<"D "<<time_<<" "
1861  <<pos2[0]-pos1[0]<<" "<<pos2[1]-pos1[1]<<" "<<pos2[2]-pos1[2]<<" "
1862  <<vel2[0]-vel1[0]<<" "<<vel2[1]-vel1[1]<<" "<<vel2[2]-vel1[2]<<" "
1863  <<ekin_<<" "<<epot_<<" "<<etot_ref_<<std::endl;
1864 #endif
1865 
1866  // kick total energy and time transformation factor for drift
1867  etot_ref_ += 2.0*dt * (mass1* (vel1[0] * pert1[0] +
1868  vel1[1] * pert1[1] +
1869  vel1[2] * pert1[2]) +
1870  mass2* (vel2[0] * pert2[0] +
1871  vel2[1] * pert2[1] +
1872  vel2[2] * pert2[2]));
1873 
1874 #ifdef AR_TTL
1875  // back up gt_kick_inv
1876  gt_kick_inv_ = gt_inv;
1877 
1878 #if (defined AR_SLOWDOWN_ARRAY ) || (defined AR_SLOWDOWN_TREE)
1879  // integrate gt_drift_inv
1880  gt_drift_inv_ += 2.0*dt*kappa_inv*kappa_inv* (vel1[0] * gtgrad1[0] +
1881  vel1[1] * gtgrad1[1] +
1882  vel1[2] * gtgrad1[2] +
1883  vel2[0] * gtgrad2[0] +
1884  vel2[1] * gtgrad2[1] +
1885  vel2[2] * gtgrad2[2]);
1886 #else
1887  // integrate gt_drift_inv
1888  gt_drift_inv_ += 2.0*dt* (vel1[0] * gtgrad1[0] +
1889  vel1[1] * gtgrad1[1] +
1890  vel1[2] * gtgrad1[2] +
1891  vel2[0] * gtgrad2[0] +
1892  vel2[1] * gtgrad2[1] +
1893  vel2[2] * gtgrad2[2]);
1894 #endif
1895 
1896 #endif // AR_TTL
1897 
1898  // kick half step for velocity
1899  vel1[0] += dvel1[0];
1900  vel1[1] += dvel1[1];
1901  vel1[2] += dvel1[2];
1902 
1903  vel2[0] += dvel2[0];
1904  vel2[1] += dvel2[1];
1905  vel2[2] += dvel2[2];
1906 
1907  // calculate kinetic energy
1908  ekin_ = 0.5 * (mass1 * (vel1[0]*vel1[0]+vel1[1]*vel1[1]+vel1[2]*vel1[2]) +
1909  mass2 * (vel2[0]*vel2[0]+vel2[1]*vel2[1]+vel2[2]*vel2[2]));
1910 
1911  }
1912 
1913 #if (defined AR_SLOWDOWN_ARRAY ) || (defined AR_SLOWDOWN_TREE)
1914  // make consistent slowdown inner energy
1915  etot_sd_ref_ = etot_ref_*kappa_inv;
1916  ekin_sd_ = ekin_*kappa_inv;
1917  epot_sd_ = epot_*kappa_inv;
1918 #endif
1919  }
1920 
1921  // Integrate the system to a given time
1926  InterruptBinary<Tparticle> integrateToTime(const Float _time_end) {
1927  ASSERT(checkParams());
1928 
1929  // real full time step
1930  const Float dt_full = _time_end - time_;
1931 
1932  // time error
1933  const Float time_error = manager->time_error_max;
1934 
1935  // energy error limit
1936  const Float energy_error_rel_max = manager->energy_error_relative_max;
1937  // expect energy_error using half step if energy_error_rel_max reached
1938  //const Float energy_error_rel_max_half_step = energy_error_rel_max * manager->step.calcErrorRatioFromStepModifyFactor(0.5);
1939  //const Float dt_min = manager->time_step_min;
1940 
1941  // backup data size
1942  const int bk_data_size = getBackupDataSize();
1943 
1944  Float backup_data[bk_data_size]; // for backup chain data
1945 #ifdef AR_DEBUG_DUMP
1946  Float backup_data_init[bk_data_size]; // for backup initial data
1947 #endif
1948  bool backup_flag=true; // flag for backup or restore
1949 
1950  // time table
1951  const int cd_pair_size = manager->step.getCDPairSize();
1952  Float time_table[cd_pair_size]; // for storing sub-integrated time
1953 
1954  // two switch step control
1955  Float ds[2] = {info.ds,info.ds}; // step with a buffer
1956  Float ds_init = info.ds; //backup initial step
1957  int ds_switch=0; // 0 or 1
1958 
1959  // reduce ds control, three level
1960  const int n_reduce_level_max=10;
1961 
1962  struct DsBackupManager{
1963  Float ds_backup[n_reduce_level_max+1];
1964  int n_step_wait_recover_ds[n_reduce_level_max+1];
1965  int n_reduce_level;
1966 
1967  DsBackupManager(const Float _ds) { initial(_ds), n_reduce_level = -1; }
1968 
1969  void initial(const Float _ds) {
1970  for (int i=0; i<n_reduce_level_max; i++) {
1971  ds_backup[i] = _ds;
1972  n_step_wait_recover_ds[i] = 0;
1973  }
1974  }
1975 
1976  // shift backup by one level, if successful, return true
1977  bool shiftReduceLevel() {
1978  if (n_reduce_level>0) {
1979  for (int i=0; i<n_reduce_level-1; i++) {
1980  ds_backup[i] =ds_backup[i+1];
1981  n_step_wait_recover_ds[i] = n_step_wait_recover_ds[i+1];
1982  }
1983  n_reduce_level--;
1984  return true;
1985  }
1986  return false;
1987  }
1988 
1989  // record ds to backup
1990  void backup(const Float _ds, const Float _modify_factor) {
1991  //if (n_reduce_level==n_reduce_level_max) shiftReduceLevel(); // not converse sometimes, becomes infinite small steps
1992  if (n_reduce_level==n_reduce_level_max)
1993  n_step_wait_recover_ds[n_reduce_level] += 2*to_int(1.0/_modify_factor);
1994  else {
1995  n_reduce_level++;
1996  ds_backup[n_reduce_level] = _ds;
1997  n_step_wait_recover_ds[n_reduce_level] = 2*to_int(1.0/_modify_factor);
1998  }
1999  }
2000 
2001  // count step (return false) and recover ds if necessary (return true)
2002  bool countAndRecover(Float &_ds, Float &_modify_factor, const bool _recover_flag) {
2003  if (n_reduce_level>=0) {
2004  if (n_step_wait_recover_ds[n_reduce_level] ==0) {
2005  if (_recover_flag) {
2006  _modify_factor = ds_backup[n_reduce_level]/_ds;
2007  _ds = ds_backup[n_reduce_level];
2008  n_step_wait_recover_ds[n_reduce_level] = -1;
2009  n_reduce_level--;
2010  return true;
2011  }
2012  }
2013  else {
2014  n_step_wait_recover_ds[n_reduce_level]--;
2015  return false;
2016  }
2017  }
2018  return false;
2019  }
2020 
2021  } ds_backup(info.ds);
2022 
2023  int reduce_ds_count=0; // number of reduce ds (ignore first few steps)
2024  Float step_modify_factor=1.0; // step modify factor
2025  Float previous_step_modify_factor=1.0; // step modify factor
2026  Float previous_error_ratio=-1; // previous error ratio when ds is reduced
2027  bool previous_is_restore=false; // previous step is reduced or not
2028 
2029  // time end control
2030  int n_step_end=0; // number of steps integrated to reach the time end for one during the time sychronization sub steps
2031  bool time_end_flag=false; // indicate whether time reach the end
2032 
2033  // step count
2034  long long unsigned int step_count=0; // integration step
2035  long long unsigned int step_count_tsyn=0; // time synchronization step
2036  InterruptBinary<Tparticle> bin_interrupt;
2037  bin_interrupt.time_now=time_ + info.time_offset;
2038  bin_interrupt.time_end=_time_end + info.time_offset;
2039  InterruptBinary<Tparticle> bin_interrupt_return = bin_interrupt;
2040 
2041  // particle data
2042  const int n_particle = particles.getSize();
2043 
2044 /* This must suppress since after findslowdowninner, slowdown inner is reset to 1.0, recalculate ekin_sdi give completely wrong value for energy correction for slowdown change later
2045 #ifdef AR_DEBUG
2046  Float ekin_check = ekin_;
2047  calcEKin();
2048  ASSERT(abs(ekin_check-ekin_)<1e-10);
2049  ekin_ = ekin_check;
2050 #endif
2051 */
2052  // warning print flag
2053  bool warning_print_once=true;
2054 
2055 #ifdef AR_DEBUG_DUMP
2056  // back up initial data
2057  backupIntData(backup_data_init);
2058 #endif
2059 
2060 #ifdef AR_SLOWDOWN_ARRAY
2061  // find new inner slowdown binaries, the binary tree data may be modified, thus it is safer to recheck slowdown inner binary at beginning to avoid memory issue (bin is pointer).
2062 #ifdef AR_TTL
2063  if (n_particle >2) {
2064  int nold = binary_slowdown.getSize();
2065  findSlowDownInner(time_);
2066  int nnew = binary_slowdown.getSize();
2067  // in case slowdown is disabled in the next step, gt_drift_inv_ should be re-initialized
2068  if (nold>0&&nnew==0) {
2069  gt_kick_inv_ = manager->interaction.calcAccPotAndGTKickInv(force_.getDataAddress(), epot_, particles.getDataAddress(), particles.getSize(), particles.cm, perturber, time_);
2070  gt_drift_inv_ = gt_kick_inv_;
2071  }
2072  }
2073 #else
2074  if (n_particle >2) findSlowDownInner(time_);
2075 #endif
2076 #endif
2077 
2078 #ifdef AR_SLOWDOWN_TREE
2079  // update slowdown and correct slowdown energy and gt_inv
2080  updateSlowDownAndCorrectEnergy(true, true);
2081 #endif
2082 
2083 
2084  // reset binary stab_check_time
2085  for (int i=0; i<info.binarytree.getSize(); i++)
2086  info.binarytree[i].stab_check_time = time_;
2087 
2088  // integration loop
2089  while(true) {
2090  // backup data
2091  bool binary_update_flag=false;
2092  auto& bin_root = info.getBinaryTreeRoot();
2093  auto& G = manager->interaction.gravitational_constant;
2094 
2095  if(backup_flag) {
2096  // check interrupt condiction, ensure that time end not reach
2097  if (manager->interrupt_detection_option>0 && !time_end_flag) {
2098  bin_interrupt.time_now = time_ + info.time_offset;
2099  bin_interrupt.time_end = _time_end + info.time_offset;
2100  // calc perturbation energy
2101  //Float epert=0.0;
2102  //for (int i=0; i<n_particle; i++) {
2103  // epert += force_[i].pot_pert*particles[i].mass;
2104  //}
2105  manager->interaction.modifyAndInterruptIter(bin_interrupt, bin_root);
2106  //InterruptBinary<Tparticle>* bin_intr_ptr = &bin_interrupt;
2107  //bin_intr_ptr = bin_root.processRootIter(bin_intr_ptr, Tmethod::modifyAndInterruptIter);
2108  ASSERT(bin_interrupt.checkParams());
2109  if (bin_interrupt.status!=InterruptStatus::none) {
2110  // the mode return back to the root scope
2111  if (manager->interrupt_detection_option==2) {
2112  // cumulative step count
2113  profile.step_count = step_count;
2114  profile.step_count_tsyn = step_count_tsyn;
2115  profile.step_count_sum += step_count;
2116  profile.step_count_tsyn_sum += step_count_tsyn;
2117 
2118  return bin_interrupt;
2119  }
2120  else {
2121 
2122  // check whether destroy appears (all masses becomes zero)
2123  if (bin_interrupt.status==InterruptStatus::destroy) {
2124  // all particles become zero masses
2125 #ifdef AR_DEBUG
2126  for (int j=0; j<n_particle; j++) {
2127  ASSERT(particles[j].mass==0.0);
2128  }
2129 #endif
2130  de_change_interrupt_ -= etot_ref_;
2131  dH_change_interrupt_ -= getH();
2132  ekin_ = epot_ = etot_ref_ = 0.0;
2133 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
2134  de_sd_change_cum_ -= etot_sd_ref_;
2135  dH_sd_change_interrupt_ -= getHSlowDown();
2136  ekin_sd_ = epot_sd_ = etot_sd_ref_ = 0.0;
2137 #endif
2138 #ifdef AR_DEBUG_PRINT
2139  std::cerr<<"Interrupt condition triggered! Destroy";
2140  std::cerr<<" Time: "<<time_;
2141  bin_interrupt.adr->printColumnTitle(std::cerr);
2142  std::cerr<<std::endl;
2143  bin_interrupt.adr->printColumn(std::cerr);
2144  std::cerr<<std::endl;
2145  Tparticle::printColumnTitle(std::cerr);
2146  std::cerr<<std::endl;
2147  for (int j=0; j<2; j++) {
2148  bin_interrupt.adr->getMember(j)->printColumn(std::cerr);
2149  std::cerr<<std::endl;
2150  }
2151 #endif
2152 
2153  // set binary tree mass to zero
2154  setBinaryCMZeroIter(bin_root);
2155 
2156  // cumulative step count
2157  profile.step_count = step_count;
2158  profile.step_count_tsyn = step_count_tsyn;
2159  profile.step_count_sum += step_count;
2160  profile.step_count_tsyn_sum += step_count_tsyn;
2161 
2162  Float dt = _time_end - time_;
2163  time_ += dt;
2164 
2165  return bin_interrupt;
2166  }
2167 
2168  Float ekin_bk = ekin_;
2169  Float epot_bk = epot_;
2170  Float H_bk = getH();
2171 
2172 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
2173  Float ekin_sd_bk = ekin_sd_;
2174  Float epot_sd_bk = epot_sd_;
2175  Float H_sd_bk = getHSlowDown();
2176 #endif
2177 
2178  // update binary tree mass
2179  info.generateBinaryTree(particles, G);
2180  //updateBinaryCMIter(bin_root);
2181  //updateBinarySemiEccPeriodIter(bin_root, G, time_, true);
2182  binary_update_flag = true;
2183  //bool stable_check=
2184  //if (stable_check) bin_root.stableCheckIter(bin_root, 10000*bin_root.period);
2185 
2186  // should do later, original mass still needed
2187  //particles.cm.mass += bin_interrupt.dm;
2188 
2189 #ifdef AR_TTL
2190  Float gt_kick_inv_new = calcAccPotAndGTKickInv();
2191  Float d_gt_kick_inv = gt_kick_inv_new - gt_kick_inv_;
2192  // when the change is large, initialize gt_drift_inv_ to avoid large error
2193  if (fabs(d_gt_kick_inv)/std::max(fabs(gt_kick_inv_),fabs(gt_kick_inv_new)) >1e-3)
2194  gt_drift_inv_ = gt_kick_inv_new;
2195  else
2196  gt_drift_inv_ += d_gt_kick_inv;
2197  gt_kick_inv_ = gt_kick_inv_new;
2198 #else
2199  calcAccPotAndGTKickInv();
2200 #endif
2201  // calculate kinetic energy
2202  calcEKin();
2203 
2204  // Notice initially etot_ref_ does not include epert. The perturbation effect is accumulated in the integration. Here instance change of mass does not create any work. So no need to add de_pert
2205  // get perturbation energy change due to mass change
2206  //Float epert_new = 0.0;
2207  //for (int i=0; i<n_particle; i++) {
2208  // epert_new += force_[i].pot_pert*particles[i].mass;
2209  //}
2210  //Float de_pert = epert_new - epert; // notice this is double perturbation potential
2211 
2212  // get energy change
2213  Float de = (ekin_ - ekin_bk) + (epot_ - epot_bk); //+ de_pert;
2214  etot_ref_ += de;
2215  de_change_interrupt_ += de;
2216  dH_change_interrupt_ += getH() - H_bk;
2217 
2218 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
2219  Float de_sd = (ekin_sd_ - ekin_sd_bk) + (epot_sd_ - epot_sd_bk);// + de_pert;
2220  etot_sd_ref_ += de_sd;
2221 
2222  Float dH_sd = getHSlowDown() - H_sd_bk;
2223 
2224  // add slowdown change to the global slowdown energy
2225  de_sd_change_interrupt_ += de_sd;
2226  dH_sd_change_interrupt_ += dH_sd;
2227  de_sd_change_cum_ += de_sd;
2228  dH_sd_change_cum_ += dH_sd;
2229 #endif //SLOWDOWN
2230 
2231 #ifdef AR_DEBUG_PRINT
2232  std::cerr<<"Interrupt condition triggered!";
2233  std::cerr<<" Time: "<<time_;
2234 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
2235  std::cerr<<" Energy change: dE_SD: "<<de_sd<<" dH_SD: "<<dH_sd;
2236  std::cerr<<" Slowdown: "<<bin_root.slowdown.getSlowDownFactor()<<std::endl;
2237 #endif
2238  bin_interrupt.adr->printColumnTitle(std::cerr);
2239  std::cerr<<std::endl;
2240  bin_interrupt.adr->printColumn(std::cerr);
2241  std::cerr<<std::endl;
2242  Tparticle::printColumnTitle(std::cerr);
2243  std::cerr<<std::endl;
2244  for (int j=0; j<2; j++) {
2245  bin_interrupt.adr->getMember(j)->printColumn(std::cerr);
2246  std::cerr<<std::endl;
2247  }
2248 #endif
2249 
2250  // change fix step option to make safety if energy change is large
2251  //info.fix_step_option=FixStepOption::none;
2252 
2253  // if time_end flag set, reset it to be safety
2254  //time_end_flag = false;
2255 
2256  // check merger case
2257  if (bin_interrupt.status==InterruptStatus::merge) {
2258  // count particle having mass
2259  int count_mass=0;
2260  int index_mass_last=-1;
2261  for (int j=0; j<n_particle; j++) {
2262  if (particles[j].mass>0.0) {
2263  count_mass++;
2264  index_mass_last=j;
2265  }
2266  }
2267  // only one particle has mass, drift directly
2268  if (count_mass==1) {
2269  ASSERT(index_mass_last<n_particle&&index_mass_last>=0);
2270  auto& p = particles[index_mass_last];
2271  Float dt = _time_end - time_;
2272  p.pos[0] += dt * p.vel[0];
2273  p.pos[1] += dt * p.vel[1];
2274  p.pos[2] += dt * p.vel[2];
2275 
2276  // cumulative step count
2277  profile.step_count = step_count;
2278  profile.step_count_tsyn = step_count_tsyn;
2279  profile.step_count_sum += step_count;
2280  profile.step_count_tsyn_sum += step_count_tsyn;
2281 
2282  time_ += dt;
2283 
2284  return bin_interrupt;
2285  }
2286  // if only two particles have mass, switch off auto ds adjustment
2287  if (count_mass==2) {
2288  info.fix_step_option=FixStepOption::later;
2289  }
2290  //else {
2291  // info.generateBinaryTree(particles, G);
2292  //}
2293  }
2294 
2295 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
2296  updateSlowDownAndCorrectEnergy(true, true);
2297 #endif
2298 
2299  info.ds = info.calcDsKeplerBinaryTree(*bin_interrupt.adr, manager->step.getOrder(), G, manager->ds_scale);
2300  Float ds_max = manager->step.calcStepModifyFactorFromErrorRatio(2.0)*ds_init;
2301  Float ds_min = manager->step.calcStepModifyFactorFromErrorRatio(0.5)*ds_init;
2302  if (info.ds>ds_max || info.ds<ds_min) {
2303 #ifdef AR_DEBUG_PRINT
2304  std::cerr<<"Change ds after interruption: ds(init): "<<ds_init<<" ds(new): "<<info.ds<<" ds(now): "<<ds[0]<<std::endl;
2305 #endif
2306  ASSERT(info.ds>0);
2307  ds[0] = std::min(ds[0], info.ds);
2308  ds[1] = std::min(ds[1], info.ds);
2309  ds_backup.initial(info.ds);
2310  ds_init = info.ds;
2311  }
2312  else info.ds = ds_init;
2313 
2314  // return one should be the top root
2315  if (bin_interrupt_return.status!=InterruptStatus::none) {
2316  if (bin_interrupt_return.adr!= bin_interrupt.adr) {
2317  // give root address if interrupted binaries are different from previous one
2318  bin_interrupt_return.adr = &(info.getBinaryTreeRoot());
2319  }
2320  if (bin_interrupt.status==InterruptStatus::merge)
2321  bin_interrupt_return.status = InterruptStatus::merge;
2322  }
2323  else bin_interrupt_return = bin_interrupt;
2324  }
2325  bin_interrupt.clear();
2326  }
2327  }
2328 
2329 
2330  // update binary orbit and ds if unstable
2331  if (!time_end_flag&&!binary_update_flag) {
2332  bool update_flag=updateBinarySemiEccPeriodIter(bin_root, G, time_);
2333 
2334 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
2335  updateSlowDownAndCorrectEnergy(true, true);
2336 #endif
2337 
2338  if (update_flag) {
2339  // update slowdown and correct slowdown energy and gt_inv
2340 
2341 #ifdef AR_DEBUG_PRINT
2342  std::cerr<<"Update binary tree orbits, time= "<<time_<<"\n";
2343 #endif
2344  info.ds = info.calcDsKeplerBinaryTree(bin_root, manager->step.getOrder(), G, manager->ds_scale);
2345  if (abs(ds_init-info.ds)/ds_init>0.1) {
2346 #ifdef AR_DEBUG_PRINT
2347  std::cerr<<"Change ds after update binary orbit: ds(init): "<<ds_init<<" ds(new): "<<info.ds<<" ds(now): "<<ds[0]<<std::endl;
2348 #endif
2349  ASSERT(info.ds>0);
2350  ds[0] = std::min(ds[0], info.ds);
2351  ds[1] = std::min(ds[1], info.ds);
2352  ds_backup.initial(info.ds);
2353  ds_init = info.ds;
2354  }
2355  }
2356  }
2357 
2358  int bk_return_size = backupIntData(backup_data);
2359  ASSERT(bk_return_size == bk_data_size);
2360  (void)bk_return_size;
2361 
2362  }
2363  else { //restore data
2364  int bk_return_size = restoreIntData(backup_data);
2365  ASSERT(bk_return_size == bk_data_size);
2366  (void)bk_return_size;
2367 //#ifdef AR_SLOWDOWN_ARRAY
2368  // update c.m. of binaries
2369  // binary c.m. is not backup, thus recalculate to get correct c.m. velocity for position drift correction due to slowdown inner (the first drift in integrateonestep assume c.m. vel is up to date)
2370 // updateCenterOfMassForBinaryWithSlowDownInner();
2371 //#elif AR_SLOWDOWN_TREE
2372  updateBinaryCMIter(info.getBinaryTreeRoot());
2373 //#endif
2374  }
2375 
2376  // get real time
2377  Float dt = time_;
2378 
2379  // integrate one step
2380  ASSERT(!ISINF(ds[ds_switch]));
2381  if(n_particle==2) integrateTwoOneStep(ds[ds_switch], time_table);
2382  else integrateOneStep(ds[ds_switch], time_table);
2383  //info.generateBinaryTree(particles, G);
2384 
2385  // real step size
2386  dt = time_ - dt;
2387 // ASSERT(dt>0.0);
2388 
2389  step_count++;
2390 
2391  // energy check
2392 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
2393  Float energy_error_bk = getEnergyErrorSlowDownFromBackup(backup_data);
2394  Float etot_ref_bk = getEtotSlowDownRefFromBackup(backup_data);
2395  Float energy_error = getEnergyErrorSlowDown();
2396  Float H_bk = getHSlowDownFromBackup(backup_data);
2397  Float H = getHSlowDown();
2398 #else
2399  Float energy_error_bk = getEnergyErrorFromBackup(backup_data);
2400  Float etot_ref_bk = getEtotRefFromBackup(backup_data);
2401  Float energy_error = getEnergyError();
2402  Float H_bk = getHFromBackup(backup_data);
2403  Float H = getH();
2404 #endif
2405  Float energy_error_diff = energy_error - energy_error_bk;
2406 
2407  Float energy_error_rel_abs = abs(energy_error_diff/etot_ref_bk);
2408 
2409  // get integration error for extended Hamiltonian
2410  Float integration_error_rel_abs = abs(H-H_bk);
2411  // H should be zero initially
2412  Float integration_error_rel_cum_abs = abs(H);
2413 
2414  Float integration_error_ratio = energy_error_rel_max/integration_error_rel_abs;
2415 
2416  // time error
2417  Float time_diff_rel = (_time_end - time_)/dt_full;
2418 
2420  auto regularStepFactor = [](const Float _fac) {
2421  Float fac = 1.0;
2422  if (_fac<1) while (fac>_fac) fac *= 0.5;
2423  else {
2424  while (fac<=_fac) fac *= 2.0;
2425  fac *= 0.5;
2426  }
2427  return fac;
2428  };
2429 
2430  Float error_increase_ratio_regular = manager->step.calcErrorRatioFromStepModifyFactor(2.0);
2431 
2432  // error message print
2433  auto printMessage = [&](const char* message) {
2434  std::cerr<<message<<std::endl;
2435  std::cerr<<" T: "<<time_
2436  <<" dT_err/T: "<<time_diff_rel
2437  <<" ds: "<<ds[ds_switch]
2438  <<" ds_init: "<<ds_init
2439  <<" |Int_err/E|: "<<integration_error_rel_abs
2440  <<" |Int_err_cum/E|: "<<integration_error_rel_cum_abs
2441  <<" |dE/E|: "<<energy_error_rel_abs
2442  <<" dE_cum: "<<energy_error
2443  <<" Etot_sd: "<<etot_ref_bk
2444  <<" T_end_flag: "<<time_end_flag
2445  <<" Step_count: "<<step_count;
2446  switch (info.fix_step_option) {
2447  case FixStepOption::always:
2448  std::cerr<<" Fix: always"<<std::endl;
2449  break;
2450  case FixStepOption::later:
2451  std::cerr<<" Fix: later"<<std::endl;
2452  break;
2453  case FixStepOption::none:
2454  std::cerr<<" Fix: none"<<std::endl;
2455  break;
2456  default:
2457  break;
2458  }
2459  };
2460 
2461 #ifdef AR_COLLECT_DS_MODIFY_INFO
2462  auto collectDsModifyInfo = [&](const char* error_message) {
2463  std::cerr<<error_message<<": "
2464  <<"time "<<time_<<" "
2465  <<"ds_new "<<ds[1-ds_switch]<<" "
2466  <<"ds_init "<<ds_init<<" "
2467  <<"modify "<<step_modify_factor<<" "
2468  <<"steps "<<step_count<<" "
2469  <<"n_mods "<<reduce_ds_count<<" "
2470  <<"err "<<integration_error_rel_abs<<" "
2471  <<"err/max "<<1.0/integration_error_ratio<<" "
2472  <<"errcum/E "<<integration_error_rel_cum_abs<<" "
2473  <<"dt "<<dt<<" "
2474  <<"n_ptcl "<<n_particle<<" ";
2475  for (int i=0; i<info.binarytree.getSize(); i++) {
2476  auto& bini = info.binarytree[i];
2477  std::cerr<<"semi "<<bini.semi<<" "
2478  <<"ecc "<<bini.ecc<<" "
2479  <<"period "<<bini.period<<" "
2480  <<"m1 "<<bini.m1<<" "
2481  <<"m2 "<<bini.m2<<" "
2482  <<"stab "<<bini.stab<<" "
2483  <<"sd "<<bini.slowdown.getSlowDownFactor()<<" "
2484  <<"sd_org "<<bini.slowdown.getSlowDownFactorOrigin()<<" "
2485  <<"pert_in "<<bini.slowdown.getPertIn()<<" "
2486  <<"pert_out "<<bini.slowdown.getPertOut()<<" ";
2487  }
2488  std::cerr<<std::endl;
2489  };
2490 #endif
2491 
2492 //#ifdef AR_WARN
2493  // warning for large number of steps
2494  if(warning_print_once&&step_count>=manager->step_count_max) {
2495  if(step_count%manager->step_count_max==0) {
2496  printMessage("Warning: step count is signficiant large");
2497  for (int i=0; i<info.binarytree.getSize(); i++){
2498  auto& bin = info.binarytree[i];
2499  std::cerr<<" Binary["<<i<<"]: "
2500  <<" i1="<<bin.getMemberIndex(0)
2501  <<" i2="<<bin.getMemberIndex(1)
2502  <<" m1="<<bin.m1
2503  <<" m2="<<bin.m2
2504  <<" semi= "<<bin.semi
2505  <<" ecc= "<<bin.ecc
2506  <<" period= "<<bin.period
2507  <<" stab= "<<bin.stab
2508  <<" SD= "<<bin.slowdown.getSlowDownFactor()
2509  <<" SD_org= "<<bin.slowdown.getSlowDownFactorOrigin()
2510  <<" Tscale= "<<bin.slowdown.timescale
2511  <<" pert_in= "<<bin.slowdown.pert_in
2512  <<" pert_out= "<<bin.slowdown.pert_out;
2513  std::cerr<<std::endl;
2514  warning_print_once = false;
2515  }
2516  //printColumnTitle(std::cerr,20,info.binarytree.getSize());
2517  //std::cerr<<std::endl;
2518  //printColumn(std::cerr,20,info.binarytree.getSize());
2519  //std::cerr<<std::endl;
2520 #ifdef AR_DEBUG_DUMP
2521  if (!info.dump_flag) {
2522  DATADUMP("dump_large_step");
2523  info.dump_flag=true;
2524  }
2525 #endif
2526 
2527 // // increase step size if energy error is small, not works correctly, suppress
2528 // if(integration_error_rel_abs<energy_error_rel_max) {
2529 // Float integration_error_ratio = energy_error_rel_max/integration_error_rel_abs;
2530 // Float step_modify_factor = manager->step.calcStepModifyFactorFromErrorRatio(integration_error_ratio);
2531 // ASSERT(step_modify_factor>0.0);
2532 // ds[ds_switch] *= step_modify_factor;
2533 // info.ds = ds[ds_switch];
2534 // ds[1-ds_switch] = ds[ds_switch];
2535 // ds_backup.initial(info.ds);
2536 // ds_init = info.ds;
2537 // ASSERT(!ISINF(ds[ds_switch]));
2538 //#ifdef AR_DEBUG_PRINT
2539 // std::cerr<<"Energy error is small enough for increase step, integration_error_rel_abs="<<integration_error_rel_abs
2540 // <<" energy_error_rel_max="<<energy_error_rel_max<<" step_modify_factor="<<step_modify_factor<<" new ds="<<ds[1-ds_switch]<<std::endl;
2541 //#endif
2542 // }
2543  }
2544  }
2545 //#endif
2546 
2547  // When time sychronization steps too large, abort
2548  if(step_count_tsyn>manager->step_count_max) {
2549  printMessage("Error! step count after time synchronization is too large");
2550  printColumnTitle(std::cerr,20,info.binarytree.getSize());
2551  std::cerr<<std::endl;
2552  printColumn(std::cerr,20,info.binarytree.getSize());
2553  std::cerr<<std::endl;
2554 // restoreIntData(backup_data_init);
2555 #ifdef AR_DEBUG_DUMP
2556  if (!info.dump_flag) {
2557  DATADUMP("dump_large_step");
2558  info.dump_flag=true;
2559  }
2560 #endif
2561  abort();
2562  }
2563 
2564 
2565 #ifdef AR_DEEP_DEBUG
2566  printMessage("");
2567  std::cerr<<"Timetable: ";
2568  for (int i=0; i<cd_pair_size; i++) std::cerr<<" "<<time_table[manager->step.getSortCumSumCKIndex(i)];
2569  std::cerr<<std::endl;
2570 #endif
2571 
2572  ASSERT(!ISNAN(integration_error_rel_abs));
2573 
2574  // modify step if energy error is large
2575  if(integration_error_rel_abs>energy_error_rel_max && info.fix_step_option!=FixStepOption::always) {
2576 
2577  bool check_flag = true;
2578 
2579  // check whether already modified
2580  if (previous_step_modify_factor!=1.0) {
2581  ASSERT(previous_error_ratio>0.0);
2582 
2583  // if error does not reduce much, do not modify step anymore
2584  if (integration_error_ratio>0.5*previous_error_ratio) check_flag=false;
2585  }
2586 
2587  if (check_flag) {
2588  // for initial steps, reduce step permanently
2589  if(step_count<5) {
2590 
2591  // estimate the modification factor based on the symplectic order
2592  // limit step_modify_factor to 0.125
2593  step_modify_factor = std::max(regularStepFactor(manager->step.calcStepModifyFactorFromErrorRatio(integration_error_ratio)), Float(0.125));
2594  ASSERT(step_modify_factor>0.0);
2595 
2596  previous_step_modify_factor = step_modify_factor;
2597  previous_error_ratio = integration_error_ratio;
2598 
2599  ds[ds_switch] *= step_modify_factor;
2600  ds[1-ds_switch] = ds[ds_switch];
2601  // permanently reduce ds
2602  // info.ds = ds[ds_switch];
2603  // ASSERT(!ISINF(info.ds));
2604  ds_backup.initial(info.ds);
2605 
2606  backup_flag = false;
2607 #ifdef AR_COLLECT_DS_MODIFY_INFO
2608  collectDsModifyInfo("Large_energy_error");
2609 #endif
2610  continue;
2611  }
2612  // for big energy error, reduce step temparely
2613  else if (info.fix_step_option==FixStepOption::none) {
2614 
2615  // estimate the modification factor based on the symplectic order
2616  // limit step_modify_factor to 0.125
2617  step_modify_factor = std::max(regularStepFactor(manager->step.calcStepModifyFactorFromErrorRatio(integration_error_ratio)), Float(0.125));
2618  ASSERT(step_modify_factor>0.0);
2619 
2620  previous_step_modify_factor = step_modify_factor;
2621  previous_error_ratio = integration_error_ratio;
2622 
2623  ds_backup.backup(ds[ds_switch], step_modify_factor);
2624  if(previous_is_restore) reduce_ds_count++;
2625 
2626  ds[ds_switch] *= step_modify_factor;
2627  ds[1-ds_switch] = ds[ds_switch];
2628  ASSERT(!ISINF(ds[ds_switch]));
2629 
2630  // if multiple times reduction happens, permanently reduce ds
2631  //if (reduce_ds_count>3) {
2632  // bool shift_flag = ds_backup.shiftReduceLevel();
2633  // if (!shift_flag) ds_backup.initial(ds[ds_switch]);
2634  // info.ds = ds_backup.ds_backup[0];
2635  // reduce_ds_count=0;
2636  //}
2637 
2638  backup_flag = false;
2639 #ifdef AR_COLLECT_DS_MODIFY_INFO
2640  collectDsModifyInfo("Large_energy_error");
2641 #endif
2642  continue;
2643  }
2644  }
2645  }
2646 // too much output
2647 //#ifdef AR_WARN
2648 // if(integration_error_rel_abs>100.0*energy_error_rel_max) {
2649 // std::cerr<<"Warning: symplectic integrator error > 100*criterion:"<<integration_error_rel_abs<<std::endl;
2650 // }
2651 //#endif
2652 
2653  // if negative step, reduce step size
2654  if(!time_end_flag&&dt<0) {
2655  // limit step_modify_factor to 0.125
2656  step_modify_factor = std::min(std::max(regularStepFactor(manager->step.calcStepModifyFactorFromErrorRatio(abs(_time_end/dt))), Float(0.0625)),Float(0.5));
2657  ASSERT(step_modify_factor>0.0);
2658  previous_step_modify_factor = step_modify_factor;
2659  previous_error_ratio = integration_error_ratio;
2660 
2661  ds[ds_switch] *= step_modify_factor;
2662  ds[1-ds_switch] = ds[ds_switch];
2663  ASSERT(!ISINF(ds[ds_switch]));
2664 
2665  // for initial steps, reduce step permanently
2666  if (step_count<5) {
2667  //info.ds = ds[ds_switch];
2668  ds_backup.initial(info.ds);
2669  }
2670  else { // reduce step temparely
2671  ds_backup.backup(ds[ds_switch], step_modify_factor);
2672  }
2673 
2674  backup_flag = false;
2675 
2676 #ifdef AR_COLLECT_DS_MODIFY_INFO
2677  collectDsModifyInfo("Negative_step");
2678 #endif
2679  continue;
2680 // std::cerr<<"Error! symplectic integrated time step ("<<dt<<") < minimum step ("<<dt_min<<")!\n";
2681 // printMessage();
2682 //#ifdef AR_DEBUG_DUMP
2683 // DATADUMP("dump_negative_time");
2684 //#endif
2685 // abort();
2686  }
2687 
2688  // if no modification, reset previous values
2689  previous_step_modify_factor = 1.0;
2690  previous_error_ratio = -1.0;
2691 
2692  // check integration time
2693  if(time_ < _time_end - time_error){
2694  // step increase depend on n_step_wait_recover_ds
2695  if(info.fix_step_option==FixStepOption::none && !time_end_flag) {
2696  // waiting step count reach
2697  previous_is_restore=ds_backup.countAndRecover(ds[1-ds_switch], step_modify_factor, integration_error_ratio>error_increase_ratio_regular);
2698  if (previous_is_restore) {
2699  //previous_error_ratio = -1;
2700  //previous_step_modify_factor = 1.0;
2701 #ifdef AR_COLLECT_DS_MODIFY_INFO
2702  collectDsModifyInfo("Reuse_backup_ds");
2703 #endif
2704  }
2705  // increase step size if energy error is small, not works correctly, integration error may not increase when ds becomes larger, then a very large ds may appear after several iterations. suppress
2706  else if(integration_error_rel_abs<0.5*energy_error_rel_max && dt>0.0 && dt_full/dt>std::max(100.0,0.02*manager->step_count_max)) {
2707  Float integration_error_ratio = energy_error_rel_max/integration_error_rel_abs;
2708  Float step_modify_factor = std::min(Float(100.0),manager->step.calcStepModifyFactorFromErrorRatio(integration_error_ratio));
2709  ASSERT(step_modify_factor>0.0);
2710  ds[1-ds_switch] *= step_modify_factor;
2711  info.ds = ds[1-ds_switch];
2712  ASSERT(!ISINF(ds[1-ds_switch]));
2713 #ifdef AR_DEBUG_PRINT
2714  std::cerr<<"Energy error is small enough for increase step, integration_error_rel_abs="<<integration_error_rel_abs
2715  <<" energy_error_rel_max="<<energy_error_rel_max<<" step_modify_factor="<<step_modify_factor<<" new ds="<<ds[1-ds_switch]<<std::endl;
2716 #endif
2717  }
2718  }
2719 
2720  // time sychronization on case, when step size too small to reach time end, increase step size
2721  if(time_end_flag && ds[ds_switch]==ds[1-ds_switch]) {
2722  step_count_tsyn++;
2723 
2724  Float dt_end = _time_end - time_;
2725  if (dt<0) {
2726  // limit step_modify_factor to 0.125
2727  step_modify_factor = std::min(std::max(regularStepFactor(manager->step.calcStepModifyFactorFromErrorRatio(abs(_time_end/dt))), Float(0.0625)),Float(0.5));
2728  ASSERT(step_modify_factor>0.0);
2729 
2730  ds[ds_switch] *= step_modify_factor;
2731  ds[1-ds_switch] = ds[ds_switch];
2732  ASSERT(!ISINF(ds[ds_switch]));
2733  }
2734  else if (n_step_end>1 && dt<0.3*dt_end) {
2735  // dt should be >0.0
2736  // ASSERT(dt>0.0);
2737  ds[1-ds_switch] = ds[ds_switch] * dt_end/dt;
2738  ASSERT(!ISINF(ds[1-ds_switch]));
2739 #ifdef AR_DEEP_DEBUG
2740  std::cerr<<"Time step dt(real) "<<dt<<" <0.3*(time_end-time)(real) "<<dt_end<<" enlarge step factor: "<<dt_end/dt<<" new ds: "<<ds[1-ds_switch]<<std::endl;
2741 #endif
2742  }
2743  else n_step_end++;
2744  }
2745 
2746  // when used once, update to the new step
2747  ds[ds_switch] = ds[1-ds_switch];
2748  ASSERT(!ISINF(ds[ds_switch]));
2749  ds_switch = 1-ds_switch;
2750 
2751  if (dt>0) backup_flag = true;
2752  else backup_flag = false;
2753  }
2754  else if(time_ > _time_end + time_error) {
2755  time_end_flag = true;
2756  backup_flag = false;
2757 
2758  step_count_tsyn++;
2759  n_step_end=0;
2760 
2761  // check timetable
2762  int i=-1,k=0; // i indicate the increasing time index, k is the corresponding index in time_table
2763  for(i=0; i<cd_pair_size; i++) {
2764  k = manager->step.getSortCumSumCKIndex(i);
2765  if(_time_end<time_table[k]) break;
2766  }
2767  if (i==0) { // first step case
2768  ASSERT(time_table[k]>0.0);
2769  ds[ds_switch] *= manager->step.getSortCumSumCK(i)*_time_end/time_table[k];
2770  ds[1-ds_switch] = ds[ds_switch];
2771  ASSERT(!ISINF(ds[ds_switch]));
2772 #ifdef AR_DEEP_DEBUG
2773  std::cerr<<"Time_end reach, time[k]= "<<time_table[k]<<" time= "<<time_<<" time_end/time[k]="<<_time_end/time_table[k]<<" CumSum_CK="<<manager->step.getSortCumSumCK(i)<<" ds(next) = "<<ds[ds_switch]<<" ds(next_next) = "<<ds[1-ds_switch]<<"\n";
2774 #endif
2775  }
2776  else { // not first step case, get the interval time
2777  // previous integrated sub time in time table
2778  Float time_prev = time_table[manager->step.getSortCumSumCKIndex(i-1)];
2779  Float dt_k = time_table[k] - time_prev;
2780  Float ds_tmp = ds[ds_switch];
2781  // get cumsum CK factor for two steps near the time_end
2782  Float cck_prev = manager->step.getSortCumSumCK(i-1);
2783  Float cck = manager->step.getSortCumSumCK(i);
2784  // in case the time is between two sub step, first scale the next step with the previous step CumSum CK cck(i-1)
2785  ASSERT(!ISINF(cck_prev));
2786  ds[ds_switch] *= cck_prev;
2787  ASSERT(!ISINF(ds[ds_switch]));
2788  // then next next step, scale with the CumSum CK between two step: cck(i) - cck(i-1)
2789  ASSERT(dt_k>0.0);
2790  ds[1-ds_switch] = ds_tmp*(cck-cck_prev)*std::min(Float(1.0),(_time_end-time_prev+time_error)/dt_k);
2791  ASSERT(!ISINF(ds[1-ds_switch]));
2792 
2793 #ifdef AR_DEEP_DEBUG
2794  std::cerr<<"Time_end reach, time_prev= "<<time_prev<<" time[k]= "<<time_table[k]<<" time= "<<time_<<" (time_end-time_prev)/dt="<<(_time_end-time_prev)/dt<<" CumSum_CK="<<cck<<" CumSum_CK(prev)="<<cck_prev<<" ds(next) = "<<ds[ds_switch]<<" ds(next_next) = "<<ds[1-ds_switch]<<" \n";
2795 #endif
2796  }
2797  }
2798  else {
2799 #ifdef AR_DEEP_DEBUG
2800  std::cerr<<"Finish, time_diff_rel = "<<time_diff_rel<<" integration_error_rel_abs = "<<integration_error_rel_abs<<std::endl;
2801 #endif
2802 //#ifdef AR_WARN
2803 // if (integration_error_rel_cum_abs>energy_error_rel_max) {
2804 // std::cerr<<"AR large energy error at the end! ";
2805 // printMessage();
2806 //#ifdef AR_DEBUG_DUMP
2808 // DATADUMP("dump_large_error");
2809 //#endif
2810 // }
2811 //#endif
2812  break;
2813  }
2814  }
2815 
2816  // cumulative step count
2817  profile.step_count = step_count;
2818  profile.step_count_tsyn = step_count_tsyn;
2819  profile.step_count_sum += step_count;
2820  profile.step_count_tsyn_sum += step_count_tsyn;
2821 
2822  return bin_interrupt_return;
2823  }
2824 
2826 
2829  void correctCenterOfMassDrift() {
2830  ASSERT(!particles.isOriginFrame());
2831  Float mcm=0.0, pos_cm[3]={0.0,0.0,0.0}, vel_cm[3]={0.0,0.0,0.0};
2832  auto* particle_data= particles.getDataAddress();
2833  for (int i=0; i<particles.getSize(); i++) {
2834  const Float *ri = particle_data[i].pos;
2835  const Float *vi = particle_data[i].getVel();
2836  const Float mi = particle_data[i].mass;
2837 
2838  pos_cm[0] += ri[0] * mi;
2839  pos_cm[1] += ri[1] * mi;
2840  pos_cm[2] += ri[2] * mi;
2841 
2842  vel_cm[0] += vi[0] * mi;
2843  vel_cm[1] += vi[1] * mi;
2844  vel_cm[2] += vi[2] * mi;
2845  mcm += mi;
2846  }
2847  pos_cm[0] /= mcm;
2848  pos_cm[1] /= mcm;
2849  pos_cm[2] /= mcm;
2850  vel_cm[0] /= mcm;
2851  vel_cm[1] /= mcm;
2852  vel_cm[2] /= mcm;
2853 
2854  for (int i=0; i<particles.getSize(); i++) {
2855  Float *ri = particle_data[i].pos;
2856  Float *vi = particle_data[i].getVel();
2857 
2858  ri[0] -= pos_cm[0];
2859  ri[1] -= pos_cm[1];
2860  ri[2] -= pos_cm[2];
2861  vi[0] -= vel_cm[0];
2862  vi[1] -= vel_cm[1];
2863  vi[2] -= vel_cm[2];
2864  }
2865  }
2866 
2867 #ifdef AR_SLOWDOWN_ARRAY
2868 
2872  template <class Tptcl>
2873  void writeBackSlowDownParticles(const Tptcl& _particle_cm) {
2874  ASSERT(particles.getMode()==COMM::ListMode::copy);
2875  ASSERT(!particles.isOriginFrame());
2876  auto* particle_adr = particles.getOriginAddressArray();
2877  auto* particle_data= particles.getDataAddress();
2878  const Float kappa_inv = 1.0/binary_slowdown[0]->slowdown.getSlowDownFactor();
2879  for (int i=0; i<particles.getSize(); i++) {
2880  //ASSERT(particle_adr[i]->mass == particle_data[i].mass);
2881  particle_adr[i]->mass = particle_data[i].mass;
2882 
2883  particle_adr[i]->pos[0] = particle_data[i].pos[0] + _particle_cm.pos[0];
2884  particle_adr[i]->pos[1] = particle_data[i].pos[1] + _particle_cm.pos[1];
2885  particle_adr[i]->pos[2] = particle_data[i].pos[2] + _particle_cm.pos[2];
2886 
2887  particle_adr[i]->vel[0] = particle_data[i].vel[0]*kappa_inv + _particle_cm.vel[0];
2888  particle_adr[i]->vel[1] = particle_data[i].vel[1]*kappa_inv + _particle_cm.vel[1];
2889  particle_adr[i]->vel[2] = particle_data[i].vel[2]*kappa_inv + _particle_cm.vel[2];
2890  }
2891 
2892  // correct inner slowdown velocity
2893  int nsd= binary_slowdown.getSize();
2894  for (int i=1; i<nsd; i++) {
2895  auto& sdi = binary_slowdown[i];
2896  ASSERT(sdi!=NULL);
2897  Float kappa = sdi->slowdown.getSlowDownFactor();
2898  Float kappa_inv_m_one = (1.0/kappa - 1.0)*kappa_inv;
2899  Float* velcm = sdi->getVel();
2900  for (int k=0; k<2; k++) {
2901  int j = sdi->getMemberIndex(k);
2902  ASSERT(j>=0&&j<particles.getSize());
2903  Float* vel = particle_data[j].getVel();
2904 
2905  // only scale velocity referring to binary c.m.
2906  Float vrel[3] = { vel[0] - velcm[0],
2907  vel[1] - velcm[1],
2908  vel[2] - velcm[2]};
2909  particle_adr[j]->vel[0] += vrel[0] * kappa_inv_m_one;
2910  particle_adr[j]->vel[1] += vrel[1] * kappa_inv_m_one;
2911  particle_adr[j]->vel[2] += vrel[2] * kappa_inv_m_one;
2912  }
2913  }
2914  }
2915 #endif
2916 
2917 #ifdef AR_SLOWDOWN_TREE
2918 
2920 
2926  template <class Tptcl>
2927  void writeBackSlowDownParticlesIter(const Tptcl& _particle_cm, const Float* _vel_sd_up, const Float& _inv_nest_sd_up, AR::BinaryTree<Tparticle>& _bin) {
2928  Float inv_nest_sd = _inv_nest_sd_up/_bin.slowdown.getSlowDownFactor();
2929  Float* vel_cm = _bin.getVel();
2930  for (int k=0; k<2; k++) {
2931  if (_bin.isMemberTree(k)) {
2932  auto* bink = _bin.getMemberAsTree(k);
2933  Float* vel = bink->getVel();
2934  Float vel_sd[3] = {(vel[0] - vel_cm[0]) * inv_nest_sd + _vel_sd_up[0],
2935  (vel[1] - vel_cm[1]) * inv_nest_sd + _vel_sd_up[1],
2936  (vel[2] - vel_cm[2]) * inv_nest_sd + _vel_sd_up[2]};
2937  writeBackSlowDownParticlesIter(_particle_cm, vel_sd, inv_nest_sd, *bink);
2938  }
2939  else {
2940  int i = _bin.getMemberIndex(k);
2941  auto& pk = particles[i];
2942  auto* pk_adr = particles.getMemberOriginAddress(i);
2943  Float* vel = pk.getVel();
2944  Float vel_sd[3] = {(vel[0] - vel_cm[0]) * inv_nest_sd + _vel_sd_up[0],
2945  (vel[1] - vel_cm[1]) * inv_nest_sd + _vel_sd_up[1],
2946  (vel[2] - vel_cm[2]) * inv_nest_sd + _vel_sd_up[2]};
2947  pk_adr->mass = pk.mass;
2948 
2949  pk_adr->pos[0] = pk.pos[0] + _particle_cm.pos[0];
2950  pk_adr->pos[1] = pk.pos[1] + _particle_cm.pos[1];
2951  pk_adr->pos[2] = pk.pos[2] + _particle_cm.pos[2];
2952 
2953  pk_adr->vel[0] = vel_sd[0] + _particle_cm.vel[0];
2954  pk_adr->vel[1] = vel_sd[1] + _particle_cm.vel[1];
2955  pk_adr->vel[2] = vel_sd[2] + _particle_cm.vel[2];
2956  }
2957  }
2958  }
2959 
2961 
2964  template <class Tptcl>
2965  void writeBackSlowDownParticles(const Tptcl& _particle_cm) {
2967  auto& bin_root=info.getBinaryTreeRoot();
2968  Float vel_cm[3] = {0.0,0.0,0.0};
2969  Float sd_factor=1.0;
2970  writeBackSlowDownParticlesIter(_particle_cm, vel_cm, sd_factor, bin_root);
2971  }
2972 
2973 #endif
2974 
2977  template <class Tptcl>
2978  void writeBackParticlesOriginFrame() {
2979  ASSERT(particles.getMode()==COMM::ListMode::copy);
2980  auto* particle_adr = particles.getOriginAddressArray();
2981  auto* particle_data= particles.getDataAddress();
2982  if (particles.isOriginFrame()) {
2983  for (int i=0; i<particles.getSize(); i++) {
2984  *(Tptcl*)particle_adr[i] = particle_data[i];
2985  }
2986  }
2987  else {
2988  for (int i=0; i<particles.getSize(); i++) {
2989  Tptcl pc = particle_data[i];
2990 
2991  pc.pos[0] = particle_data[i].pos[0] + particles.cm.pos[0];
2992  pc.pos[1] = particle_data[i].pos[1] + particles.cm.pos[1];
2993  pc.pos[2] = particle_data[i].pos[2] + particles.cm.pos[2];
2994 
2995  pc.vel[0] = particle_data[i].vel[0] + particles.cm.vel[0];
2996  pc.vel[1] = particle_data[i].vel[1] + particles.cm.vel[1];
2997  pc.vel[2] = particle_data[i].vel[2] + particles.cm.vel[2];
2998 
2999  *(Tptcl*)particle_adr[i] = pc;
3000  }
3001  }
3002  }
3003 
3005 
3007  Float getTime() const {
3008  return time_;
3009  }
3010 
3012 
3014  Float getEkin() const {
3015  return ekin_;
3016  }
3017 
3019 
3021  Float getEpot() const {
3022  return epot_;
3023  }
3024 
3026 
3028  Float getEtotRef() const {
3029  return etot_ref_;
3030  }
3031 
3033 
3035  Float getEtot() const {
3036  return ekin_ + epot_;
3037  }
3038 
3040 
3042  Float getEpert() const {
3043  Float epert=0.0;
3044  int n_particle = particles.getSize();
3045  for (int i=0; i<n_particle; i++) {
3046  epert += force_[i].pot_pert*particles[i].mass;
3047  }
3048  return epert;
3049  }
3050 
3052 
3054  Float getEnergyError() const {
3055  return ekin_ + epot_ - etot_ref_;
3056  }
3057 
3059  Float getEnergyErrorFromBackup(Float* _bk) const {
3060  return -_bk[1] + _bk[2] + _bk[3];
3061  }
3062 
3064  Float getEtotRefFromBackup(Float* _bk) const {
3065  return _bk[1];
3066  }
3067 
3069  Float getEtotFromBackup(Float* _bk) const {
3070  return _bk[2] + _bk[3];
3071  }
3072 
3074  void resetDEChangeBinaryInterrupt() {
3075  de_change_interrupt_ = 0.0;
3076  dH_change_interrupt_ = 0.0;
3077  }
3078 
3080  Float getDEChangeBinaryInterrupt() const {
3081  return de_change_interrupt_;
3082  }
3083 
3085  Float getDHChangeBinaryInterrupt() const {
3086  return dH_change_interrupt_;
3087  }
3088 
3090  Float getH() const {
3091 #ifdef AR_TTL
3092  //return (ekin_ - etot_ref_)/gt_drift_inv_ + epot_/gt_kick_inv_;
3093  return (ekin_ + epot_ - etot_ref_)/gt_kick_inv_;
3094 #else
3095  return manager->interaction.calcH(ekin_ - etot_ref_, epot_);
3096 #endif
3097  }
3098 
3100  Float getHFromBackup(Float* _bk) const {
3101  Float& etot_ref =_bk[1];
3102  Float& ekin = _bk[2];
3103  Float& epot = _bk[3];
3104 #ifdef AR_TTL
3105 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3106  //Float& gt_drift_inv = _bk[13];
3107  Float& gt_kick_inv = _bk[14];
3108 #else
3109  //Float& gt_drift_inv = _bk[6];
3110  Float& gt_kick_inv = _bk[7];
3111 #endif
3112  return (ekin + epot - etot_ref)/gt_kick_inv;
3113  //return (ekin - etot_ref)/gt_drift_inv + epot/gt_kick_inv;
3114 #else
3115  return manager->interaction.calcH(ekin - etot_ref, epot);
3116 #endif
3117  }
3118 
3119 
3120 
3121 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3122  void resetDESlowDownChangeCum() {
3124  de_sd_change_cum_ = 0.0;
3125  dH_sd_change_cum_ = 0.0;
3126  }
3127 
3129  Float getDESlowDownChangeCum() const {
3130  return de_sd_change_cum_;
3131  }
3132 
3134  Float getDHSlowDownChangeCum() const {
3135  return dH_sd_change_cum_;
3136  }
3137 
3139  void resetDESlowDownChangeBinaryInterrupt() {
3140  de_sd_change_interrupt_ = 0.0;
3141  dH_sd_change_interrupt_ = 0.0;
3142  }
3143 
3145  Float getDESlowDownChangeBinaryInterrupt() const {
3146  return de_sd_change_interrupt_;
3147  }
3148 
3150  Float getDHSlowDownChangeBinaryInterrupt() const {
3151  return dH_sd_change_interrupt_;
3152  }
3153 
3155 
3157  Float getEkinSlowDown() const {
3158  return ekin_sd_;
3159  }
3160 
3162 
3164  Float getEpotSlowDown() const {
3165  return epot_sd_;
3166  }
3167 
3169 
3171  Float getEtotSlowDownRef() const {
3172  return etot_sd_ref_;
3173  }
3174 
3176 
3178  Float getEtotSlowDown() const {
3179  return ekin_sd_ + epot_sd_;
3180  }
3181 
3183 
3185  Float getEnergyErrorSlowDown() const {
3186  return ekin_sd_ + epot_sd_ - etot_sd_ref_;
3187  }
3188 
3190  Float getEnergyErrorSlowDownFromBackup(Float* _bk) const {
3191  return -_bk[6] + _bk[7] + _bk[8];
3192  }
3193 
3195  Float getHSlowDown() const {
3196 #ifdef AR_TTL
3197  //return (ekin_sd_ - etot_sd_ref_)/gt_drift_inv_ + epot_sd_/gt_kick_inv_;
3198  return (ekin_sd_ + epot_sd_ - etot_sd_ref_)/gt_kick_inv_;
3199 #else
3200  return manager->interaction.calcH(ekin_sd_ - etot_sd_ref_, epot_sd_);
3201 #endif
3202  }
3203 
3205  Float getHSlowDownFromBackup(Float* _bk) const {
3206  Float& etot_sd_ref =_bk[6];
3207  Float& ekin_sd = _bk[7];
3208  Float& epot_sd = _bk[8];
3209 #ifdef AR_TTL
3210  //Float& gt_drift_inv = _bk[13];
3211  Float& gt_kick_inv = _bk[14];
3212  //return (ekin_sd - etot_sd_ref)/gt_drift_inv + epot_sd/gt_kick_inv;
3213  return (ekin_sd + epot_sd - etot_sd_ref)/gt_kick_inv;
3214 #else
3215  return manager->interaction.calcH(ekin_sd - etot_sd_ref, epot_sd);
3216 #endif
3217  }
3218 
3220  Float getEtotSlowDownRefFromBackup(Float* _bk) const {
3221  return _bk[6];
3222  }
3223 
3225  Float getEtotSlowDownFromBackup(Float* _bk) const {
3226  return _bk[7] + _bk[8];
3227  }
3228 
3229 #endif
3230 
3232  int getBackupDataSize() const {
3233  int bk_size = 6;
3234 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3235  bk_size += 7;
3236  //bk_size += SlowDown::getBackupDataSize();
3237 #endif
3238 #ifdef AR_TTL
3239  bk_size += 2;
3240 #endif
3241  bk_size += particles.getBackupDataSize();
3242  return bk_size;
3243  }
3244 
3245 
3247 
3250  int backupIntData(Float* _bk) {
3251  int bk_size=0;
3252  _bk[bk_size++] = time_; //0
3253  _bk[bk_size++] = etot_ref_; //1
3254  _bk[bk_size++] = ekin_; //2
3255  _bk[bk_size++] = epot_; //3
3256  _bk[bk_size++] = de_change_interrupt_; //4
3257  _bk[bk_size++] = dH_change_interrupt_; //5
3258 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3259  _bk[bk_size++] = etot_sd_ref_; //6
3260  _bk[bk_size++] = ekin_sd_; //7
3261  _bk[bk_size++] = epot_sd_; //8
3262  _bk[bk_size++] = de_sd_change_cum_; //9
3263  _bk[bk_size++] = dH_sd_change_cum_; //10
3264  _bk[bk_size++] = de_sd_change_interrupt_; //11
3265  _bk[bk_size++] = dH_sd_change_interrupt_; //12
3266 #endif
3267 
3268 #ifdef AR_TTL
3269  _bk[bk_size++] = gt_drift_inv_; //13 / 6
3270  _bk[bk_size++] = gt_kick_inv_; //14 / 7
3271 #endif
3272 
3273  bk_size += particles.backupParticlePosVel(&_bk[bk_size]);
3274 //#if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3275 // bk_size += info.getBinaryTreeRoot().slowdown.backup(&_bk[bk_size]); // slowdownfactor
3276 //#endif
3277  return bk_size;
3278  }
3279 
3281 
3284  int restoreIntData(Float* _bk) {
3285  int bk_size = 0;
3286  time_ = _bk[bk_size++];
3287  etot_ref_ = _bk[bk_size++];
3288  ekin_ = _bk[bk_size++];
3289  epot_ = _bk[bk_size++];
3290  de_change_interrupt_= _bk[bk_size++];
3291  dH_change_interrupt_= _bk[bk_size++];
3292 
3293 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3294  etot_sd_ref_ = _bk[bk_size++];
3295  ekin_sd_ = _bk[bk_size++];
3296  epot_sd_ = _bk[bk_size++];
3297 
3298  de_sd_change_cum_= _bk[bk_size++];
3299  dH_sd_change_cum_= _bk[bk_size++];
3300  de_sd_change_interrupt_= _bk[bk_size++];
3301  dH_sd_change_interrupt_= _bk[bk_size++];
3302 #endif
3303 #ifdef AR_TTL
3304  gt_drift_inv_ = _bk[bk_size++];
3305  gt_kick_inv_ = _bk[bk_size++];
3306 #endif
3307  bk_size += particles.restoreParticlePosVel(&_bk[bk_size]);
3308 //#if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3309 // bk_size += info.getBinaryTreeRoot().slowdown.restore(&_bk[bk_size]);
3310 //#endif
3311  return bk_size;
3312  }
3313 
3314 #ifdef AR_TTL
3315 
3320  Float getGTDriftInv() const {
3321  return gt_drift_inv_;
3322  }
3323 #endif
3324 
3326 
3332  template<class Tptcl>
3333  void printGroupInfo(const int _type, std::ostream& _fout, const int _width, const Tptcl* _pcm=NULL) {
3334  auto& bin_root = info.getBinaryTreeRoot();
3335  //auto* p1 = bin_root.getLeftMember();
3336  //auto* p2 = bin_root.getRightMember();
3337 
3338  bool reset_flag = (_type==1 && bin_root.semi<0 && bin_root.ecca>0);
3339 
3340  if (info.checkAndSetBinaryPairIDIter(bin_root, reset_flag)) {
3341  if (_type==0) return; // if it is new but already existed binary, do not print
3342  else if (!reset_flag) return; // in the end case, if the system is still bound, do not print
3343  }
3344 
3345  Float pos_cm[3], vel_cm[3];
3346  auto& pcm_loc = particles.cm;
3347  if (_pcm!=NULL) {
3348  pos_cm[0] = pcm_loc.pos[0] + _pcm->pos[0];
3349  pos_cm[1] = pcm_loc.pos[1] + _pcm->pos[1];
3350  pos_cm[2] = pcm_loc.pos[2] + _pcm->pos[2];
3351  vel_cm[0] = pcm_loc.vel[0] + _pcm->vel[0];
3352  vel_cm[1] = pcm_loc.vel[1] + _pcm->vel[1];
3353  vel_cm[2] = pcm_loc.vel[2] + _pcm->vel[2];
3354  }
3355  else {
3356  pos_cm[0] = pcm_loc.pos[0];
3357  pos_cm[1] = pcm_loc.pos[1];
3358  pos_cm[2] = pcm_loc.pos[2];
3359  vel_cm[0] = pcm_loc.vel[0];
3360  vel_cm[1] = pcm_loc.vel[1];
3361  vel_cm[2] = pcm_loc.vel[2];
3362  }
3363 #pragma omp critical
3364  {
3365  _fout<<std::setw(_width)<<_type
3366  <<std::setw(_width)<<bin_root.getMemberN()
3367  <<std::setw(_width)<<time_ + info.time_offset;
3368  _fout<<std::setw(_width)<<pos_cm[0]
3369  <<std::setw(_width)<<pos_cm[1]
3370  <<std::setw(_width)<<pos_cm[2]
3371  <<std::setw(_width)<<vel_cm[0]
3372  <<std::setw(_width)<<vel_cm[1]
3373  <<std::setw(_width)<<vel_cm[2];
3374  bin_root.printBinaryTreeIter(_fout, _width);
3375  _fout<<std::endl;
3376  }
3377  //if (_type==0) { // register pair id to avoid repeating printing
3378  // p1->setBinaryPairID(p2->id);
3379  // p2->setBinaryPairID(p1->id);
3380  //}
3381  //else { // break case reset pair id
3382  // p1->setBinaryPairID(0);
3383  // p2->setBinaryPairID(0);
3384  //}
3385  }
3386 
3388 
3393  void printColumnTitle(std::ostream & _fout, const int _width=20, const int _n_sd=0) {
3394  _fout<<std::setw(_width)<<"Time"
3395  <<std::setw(_width)<<"dE"
3396  <<std::setw(_width)<<"Etot"
3397  <<std::setw(_width)<<"Ekin"
3398  <<std::setw(_width)<<"Epot"
3399  <<std::setw(_width)<<"Gt_drift"
3400  <<std::setw(_width)<<"H"
3401  <<std::setw(_width)<<"dE_intr"
3402  <<std::setw(_width)<<"dH_intr";
3403  perturber.printColumnTitle(_fout, _width);
3404  info.printColumnTitle(_fout, _width);
3405  profile.printColumnTitle(_fout, _width);
3406 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3407  _fout<<std::setw(_width)<<"dE_SD"
3408  <<std::setw(_width)<<"Etot_SD"
3409  <<std::setw(_width)<<"Ekin_SD"
3410  <<std::setw(_width)<<"Epot_SD"
3411  <<std::setw(_width)<<"dE_SDC_cum"
3412  <<std::setw(_width)<<"dH_SDC_cum"
3413  <<std::setw(_width)<<"dE_SDC_intr"
3414  <<std::setw(_width)<<"dH_SDC_intr";
3415  _fout<<std::setw(_width)<<"N_SD";
3416  for (int i=0; i<_n_sd; i++) {
3417  _fout<<std::setw(_width)<<"I1"
3418  <<std::setw(_width)<<"I2";
3419  SlowDown::printColumnTitle(_fout, _width);
3420  }
3421 #endif
3422  particles.printColumnTitle(_fout, _width);
3423  }
3424 
3426 
3431  void printColumn(std::ostream & _fout, const int _width=20, const int _n_sd=0){
3432  _fout<<std::setw(_width)<<getTime()
3433  <<std::setw(_width)<<getEnergyError()
3434  <<std::setw(_width)<<etot_ref_
3435  <<std::setw(_width)<<ekin_
3436  <<std::setw(_width)<<epot_
3437 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3438 #ifdef AR_TTL
3439  <<std::setw(_width)<<1.0/gt_drift_inv_
3440 #else
3441  <<std::setw(_width)<<1.0/manager->interaction.calcGTDriftInv(ekin_sd_-etot_sd_ref_)
3442 #endif
3443  <<std::setw(_width)<<getHSlowDown()
3444 #else
3445 #ifdef AR_TTL
3446  <<std::setw(_width)<<1.0/gt_drift_inv_
3447 #else
3448  <<std::setw(_width)<<1.0/manager->interaction.calcGTDriftInv(ekin_-etot_ref_)
3449 #endif
3450  <<std::setw(_width)<<getH()
3451 #endif
3452  <<std::setw(_width)<<de_change_interrupt_
3453  <<std::setw(_width)<<dH_change_interrupt_;
3454  perturber.printColumn(_fout, _width);
3455  info.printColumn(_fout, _width);
3456  profile.printColumn(_fout, _width);
3457 #if (defined AR_SLOWDOWN_ARRAY) || (defined AR_SLOWDOWN_TREE)
3458  _fout<<std::setw(_width)<<getEnergyErrorSlowDown()
3459  <<std::setw(_width)<<etot_sd_ref_
3460  <<std::setw(_width)<<ekin_sd_
3461  <<std::setw(_width)<<epot_sd_
3462  <<std::setw(_width)<<de_sd_change_cum_
3463  <<std::setw(_width)<<dH_sd_change_cum_
3464  <<std::setw(_width)<<de_sd_change_interrupt_
3465  <<std::setw(_width)<<dH_sd_change_interrupt_;
3466  SlowDown sd_empty;
3467 #ifdef AR_SLOWDOWN_ARRAY
3468  int n_sd_now = binary_slowdown.getSize();
3469  _fout<<std::setw(_width)<<n_sd_now;
3470  for (int i=0; i<_n_sd; i++) {
3471  if (i<n_sd_now) {
3472  _fout<<std::setw(_width)<<binary_slowdown[i]->getMemberIndex(0)
3473  <<std::setw(_width)<<binary_slowdown[i]->getMemberIndex(1);
3474  binary_slowdown[i]->slowdown.printColumn(_fout, _width);
3475  }
3476  else {
3477  _fout<<std::setw(_width)<<-1
3478  <<std::setw(_width)<<-1;
3479  sd_empty.printColumn(_fout, _width);
3480  }
3481  }
3482 #else
3483  int n_sd_now = info.binarytree.getSize();
3484  _fout<<std::setw(_width)<<n_sd_now;
3485  for (int i=0; i<_n_sd; i++) {
3486  if (i<n_sd_now) {
3487  _fout<<std::setw(_width)<<info.binarytree[i].getMemberIndex(0)
3488  <<std::setw(_width)<<info.binarytree[i].getMemberIndex(1);
3489  info.binarytree[i].slowdown.printColumn(_fout, _width);
3490  }
3491  else {
3492  _fout<<std::setw(_width)<<-1
3493  <<std::setw(_width)<<-1;
3494  sd_empty.printColumn(_fout, _width);
3495  }
3496  }
3497 #endif
3498 #endif
3499  particles.printColumn(_fout, _width);
3500  }
3501 
3503 
3505  void writeBinary(FILE *_fout) {
3506  fwrite(&time_, sizeof(Float), 1, _fout);
3507  fwrite(&etot_ref_, sizeof(Float), 1, _fout);
3508  fwrite(&ekin_, sizeof(Float), 1, _fout);
3509  fwrite(&epot_, sizeof(Float), 1, _fout);
3510 #ifdef AR_TTL
3511  fwrite(&gt_drift_inv_, sizeof(Float), 1, _fout);
3512 #endif
3513  int size = force_.getSize();
3514  fwrite(&size, sizeof(int), 1, _fout);
3515  for (int i=0; i<size; i++) force_[i].writeBinary(_fout);
3516 
3517  particles.writeBinary(_fout);
3518  perturber.writeBinary(_fout);
3519  info.writeBinary(_fout);
3520  profile.writeBinary(_fout);
3521  }
3522 
3524 
3526  void readBinary(FILE *_fin) {
3527  size_t rcount = fread(&time_, sizeof(Float), 1, _fin);
3528  rcount += fread(&etot_ref_, sizeof(Float), 1, _fin);
3529  rcount += fread(&ekin_, sizeof(Float), 1, _fin);
3530  rcount += fread(&epot_, sizeof(Float), 1, _fin);
3531  if (rcount<4) {
3532  std::cerr<<"Error: Data reading fails! requiring data number is 4, only obtain "<<rcount<<".\n";
3533  abort();
3534  }
3535 #ifdef AR_TTL
3536  rcount = fread(&gt_drift_inv_, sizeof(Float), 1, _fin);
3537  if (rcount<1) {
3538  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
3539  abort();
3540  }
3541 #endif
3542  int size;
3543  rcount = fread(&size, sizeof(int),1, _fin);
3544  if(rcount<1) {
3545  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
3546  abort();
3547  }
3548  if(size<0) {
3549  std::cerr<<"Error: array size <0 "<<size<<"<=0!\n";
3550  abort();
3551  }
3552  if (size>0) {
3553  force_.setMode(COMM::ListMode::local);
3554  force_.reserveMem(size);
3555  force_.resizeNoInitialize(size);
3556  for (int i=0; i<size; i++) force_[i].readBinary(_fin);
3557  }
3558 
3559  particles.setMode(COMM::ListMode::local);
3560  particles.readBinary(_fin);
3561  perturber.readBinary(_fin);
3562  info.readBinary(_fin);
3563  profile.readBinary(_fin);
3564  }
3565 
3566  };
3567 }
COMM::ParticleGroup::shiftToCenterOfMassFrame
void shiftToCenterOfMassFrame()
shift particle to their c.m. frame
Definition: particle_group.h:205
COMM::List::getDataAddress
Ttype * getDataAddress() const
return member data array address
Definition: list.h:170
AR::Profile::printColumnTitle
void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: AR/profile.h:50
AR::TimeTransformedSymplecticIntegrator::getEnergyErrorFromBackup
Float getEnergyErrorFromBackup(Float *_bk) const
get energy error from backup data
Definition: symplectic_integrator.h:3059
COMM::List< Tparticle >::isModified
bool isModified() const
Get modified status.
Definition: list.h:457
COMM::ParticleGroup::printColumnTitle
void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: particle_group.h:102
AR::TimeTransformedSymplecticIntegrator::integrateTwoOneStep
void integrateTwoOneStep(const Float _ds, Float _time_table[])
integration for two body one step
Definition: symplectic_integrator.h:1720
COMM::ParticleGroup::writeBinary
void writeBinary(FILE *_fout)
write particle data to files (notice original address is lost)
Definition: particle_group.h:121
AR::TimeTransformedSymplecticIntegrator::manager
TimeTransformedSymplecticManager< Tmethod > * manager
integration manager
Definition: symplectic_integrator.h:210
AR::printFeatures
void printFeatures(std::ostream &fout)
print features
Definition: symplectic_integrator.h:20
COMM::List< Tparticle >::getOriginAddressArray
Tparticle ** getOriginAddressArray() const
return member original address array
Definition: list.h:175
AR::TimeTransformedSymplecticManager::ds_scale
Float ds_scale
minimum real time step allown
Definition: symplectic_integrator.h:63
ISNAN
#define ISNAN(x)
Definition: Float.h:50
COMM::ParticleGroup::backupParticlePosVel
int backupParticlePosVel(Float *_bk)
Backup member particle position and velocity.
Definition: particle_group.h:63
COMM::List< Tparticle >::getMemberOriginAddress
Tparticle * getMemberOriginAddress(const int _index) const
return one member original address
Definition: list.h:198
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
AR::TimeTransformedSymplecticIntegrator::clear
void clear()
Clear function.
Definition: symplectic_integrator.h:263
AR
Algorithmic regularization (time transformed explicit symplectic integrator) namespace.
Definition: force.h:5
AR::Profile::step_count
UInt64 step_count
Definition: AR/profile.h:33
COMM::ParticleGroup::isOriginFrame
bool isOriginFrame() const
return true if the system is the in their origin frame
Definition: particle_group.h:281
AR::TimeTransformedSymplecticIntegrator::getDEChangeBinaryInterrupt
Float getDEChangeBinaryInterrupt() const
get cumulative energy change due to interruption
Definition: symplectic_integrator.h:3080
AR::TimeTransformedSymplecticIntegrator::TimeTransformedSymplecticIntegrator
TimeTransformedSymplecticIntegrator()
Constructor.
Definition: symplectic_integrator.h:220
AR::Profile::readBinary
void readBinary(FILE *_fin)
read class data with BINARY format and initial the array
Definition: AR/profile.h:79
AR::TimeTransformedSymplecticIntegrator::profile
Profile profile
profile to measure the performance
Definition: symplectic_integrator.h:217
AR::printDebugFeatures
void printDebugFeatures(std::ostream &fout)
print debug features
Definition: symplectic_integrator.h:41
AR::TimeTransformedSymplecticIntegrator::getHFromBackup
Float getHFromBackup(Float *_bk) const
get Hamiltonian from backup data
Definition: symplectic_integrator.h:3100
COMM::List< Tparticle >::getMode
ListMode getMode() const
Definition: list.h:155
AR::TimeTransformedSymplecticIntegrator::getEtotRefFromBackup
Float getEtotRefFromBackup(Float *_bk) const
get integrated energy from backup data
Definition: symplectic_integrator.h:3064
COMM::List::clear
void clear()
Clear function.
Definition: list.h:76
AR::Profile::step_count_tsyn
UInt64 step_count_tsyn
Definition: AR/profile.h:34
AR::InterruptBinary::time_now
Float time_now
Definition: interrupt.h:14
NUMERIC_FLOAT_MAX
const Float NUMERIC_FLOAT_MAX
Definition: Float.h:29
AR::TimeTransformedSymplecticManager::slowdown_pert_ratio_ref
Float slowdown_pert_ratio_ref
scaling factor to determine ds
Definition: symplectic_integrator.h:64
slow_down.h
AR::TimeTransformedSymplecticIntegrator::integrateToTime
InterruptBinary< Tparticle > integrateToTime(const Float _time_end)
Definition: symplectic_integrator.h:1926
AR::TimeTransformedSymplecticIntegrator::writeBackParticlesOriginFrame
void writeBackParticlesOriginFrame()
write back particles to original address
Definition: symplectic_integrator.h:2978
AR::TimeTransformedSymplecticManager::step
SymplecticStep step
class contain interaction function
Definition: symplectic_integrator.h:73
AR::InterruptBinary::time_end
Float time_end
Definition: interrupt.h:15
COMM::ParticleGroup::getBackupDataSize
int getBackupDataSize() const
get backup data size
Definition: particle_group.h:55
AR::TimeTransformedSymplecticIntegrator::perturber
Tpert perturber
perturber class
Definition: symplectic_integrator.h:215
AR::SymplecticStep::readBinary
void readBinary(FILE *_fin)
read class data with BINARY format and initial the array
Definition: symplectic_step.h:344
COMM::BinaryTree::getLeftMember
Tptcl * getLeftMember() const
get left member
Definition: binary_tree.h:963
AR::TimeTransformedSymplecticIntegrator::printColumnTitle
void printColumnTitle(std::ostream &_fout, const int _width=20, const int _n_sd=0)
print titles of class members using column style
Definition: symplectic_integrator.h:3393
AR::TimeTransformedSymplecticManager::time_step_min
Float time_step_min
maximum energy error requirement
Definition: symplectic_integrator.h:62
AR::TimeTransformedSymplecticIntegrator::info
Tinfo info
information of the system
Definition: symplectic_integrator.h:216
AR::InterruptBinary
Binary interrupt recoreder.
Definition: interrupt.h:12
AR::TimeTransformedSymplecticIntegrator::getTime
Float getTime() const
Get current physical time.
Definition: symplectic_integrator.h:3007
AR::TimeTransformedSymplecticIntegrator::particles
COMM::ParticleGroup< Tparticle, Tpcm > particles
particle group manager
Definition: symplectic_integrator.h:211
AR::TimeTransformedSymplecticManager::checkParams
bool checkParams()
check whether parameters values are correct
Definition: symplectic_integrator.h:86
AR::TimeTransformedSymplecticIntegrator::getEpot
Float getEpot() const
Get current potential energy.
Definition: symplectic_integrator.h:3021
AR::TimeTransformedSymplecticIntegrator::getBackupDataSize
int getBackupDataSize() const
get backup data size
Definition: symplectic_integrator.h:3232
AR::Profile::step_count_sum
UInt64 step_count_sum
Definition: AR/profile.h:31
AR::TimeTransformedSymplecticIntegrator
Time Transformed Symplectic integrator class for a group of particles.
Definition: symplectic_integrator.h:175
AR::TimeTransformedSymplecticIntegrator::getH
Float getH() const
get Hamiltonian
Definition: symplectic_integrator.h:3090
AR::TimeTransformedSymplecticIntegrator::readBinary
void readBinary(FILE *_fin)
read class data with BINARY format and initial the array
Definition: symplectic_integrator.h:3526
ROUND_OFF_ERROR_LIMIT
const Float ROUND_OFF_ERROR_LIMIT
Definition: Float.h:28
AR::Profile
profiling class for AR integrator
Definition: AR/profile.h:28
AR::SymplecticStep
class to manager symplectic KD steps
Definition: symplectic_step.h:11
interrupt.h
information.h
COMM::ParticleGroup::calcCenterOfMass
void calcCenterOfMass()
calculate center-of-mass
Definition: particle_group.h:252
AR::InterruptBinary::adr
AR::BinaryTree< Tparticle > * adr
Definition: interrupt.h:13
particle_group.h
COMM::ParticleGroup
Particle group class to store and manage a group of particle.
Definition: particle_group.h:13
AR::TimeTransformedSymplecticIntegrator::integrateOneStep
void integrateOneStep(const Float _ds, Float _time_table[])
integration for one step
Definition: symplectic_integrator.h:1653
COMM::List< Tparticle >::increaseSizeNoInitialize
void increaseSizeNoInitialize(const int _n)
increase size without initialization
Definition: list.h:246
COMM::BinaryTree::getMemberN
int getMemberN() const
get total number of members
Definition: binary_tree.h:801
AR::TimeTransformedSymplecticManager::energy_error_relative_max
Float energy_error_relative_max
maximum time error (absolute), should be positive and larger than round-off error
Definition: symplectic_integrator.h:61
AR::TimeTransformedSymplecticIntegrator::writeBinary
void writeBinary(FILE *_fout)
write class data with BINARY format
Definition: symplectic_integrator.h:3505
COMM::ParticleGroup::cm
Tpcm cm
Definition: particle_group.h:19
AR::Force
force class of one particle
Definition: force.h:10
COMM::List< Tparticle >::setModifiedFalse
void setModifiedFalse()
Reset modified status to false.
Definition: list.h:462
AR::Force::acc_in
Float acc_in[3]
total acceleration from all inner particle members of the group
Definition: force.h:12
COMM::List::getSize
int getSize() const
get current member number
Definition: list.h:151
Float
double Float
Definition: Float.h:25
AR::TimeTransformedSymplecticIntegrator::operator=
TimeTransformedSymplecticIntegrator & operator=(const TimeTransformedSymplecticIntegrator &_sym)
operator =
Definition: symplectic_integrator.h:301
AR::TimeTransformedSymplecticIntegrator::resetDEChangeBinaryInterrupt
void resetDEChangeBinaryInterrupt()
reset cumulative energy/hamiltonian change due to interruption
Definition: symplectic_integrator.h:3074
AR::TimeTransformedSymplecticIntegrator::restoreIntData
int restoreIntData(Float *_bk)
Restore integration data.
Definition: symplectic_integrator.h:3284
COMM::BinaryTree::calcSemiEccPeriod
void calcSemiEccPeriod(const Float &_G)
calculate Kepler orbit semi, ecc and period only
Definition: binary_tree.h:772
AR::TimeTransformedSymplecticManager::interaction
Tmethod interaction
1: detect interruption; 0: no detection
Definition: symplectic_integrator.h:72
AR::TimeTransformedSymplecticIntegrator::getEpert
Float getEpert() const
get perturbation potential energy
Definition: symplectic_integrator.h:3042
COMM::BinaryTree::isMemberTree
bool isMemberTree(const size_t i) const
Definition: binary_tree.h:957
AR::TimeTransformedSymplecticManager::slowdown_timescale_max
Float slowdown_timescale_max
slowdown perturbation /inner ratio reference factor
Definition: symplectic_integrator.h:68
AR::TimeTransformedSymplecticIntegrator::printColumn
void printColumn(std::ostream &_fout, const int _width=20, const int _n_sd=0)
print data of class members using column style
Definition: symplectic_integrator.h:3431
AR::TimeTransformedSymplecticIntegrator::getDHChangeBinaryInterrupt
Float getDHChangeBinaryInterrupt() const
get cumulative hamiltonian change due to interruption
Definition: symplectic_integrator.h:3085
AR::SymplecticStep::print
void print(std::ostream &_fout) const
print parameters
Definition: symplectic_step.h:360
COMM::BinaryTree::getMemberAsTree
BinaryTreeLocal * getMemberAsTree(const size_t i) const
get member as tree
Definition: binary_tree.h:946
AR::InterruptBinary::checkParams
bool checkParams()
Definition: interrupt.h:36
AR::TimeTransformedSymplecticManager::print
void print(std::ostream &_fout) const
print parameters
Definition: symplectic_integrator.h:150
AR::TimeTransformedSymplecticIntegrator::~TimeTransformedSymplecticIntegrator
~TimeTransformedSymplecticIntegrator()
destructor
Definition: symplectic_integrator.h:294
AR::SymplecticStep::writeBinary
void writeBinary(FILE *_fout)
write class data with BINARY format
Definition: symplectic_step.h:337
AR::TimeTransformedSymplecticIntegrator::reserveIntegratorMem
void reserveIntegratorMem()
reserve memory for force
Definition: symplectic_integrator.h:248
AR::TimeTransformedSymplecticIntegrator::correctCenterOfMassDrift
void correctCenterOfMassDrift()
correct CM drift
Definition: symplectic_integrator.h:2829
AR::Force::acc_pert
Float acc_pert[3]
perturbation from outside of the group
Definition: force.h:13
AR::Profile::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: AR/profile.h:62
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
AR::TimeTransformedSymplecticManager::writeBinary
void writeBinary(FILE *_fout)
write class data with BINARY format
Definition: symplectic_integrator.h:107
AR::Profile::writeBinary
void writeBinary(FILE *_fout)
write class data with BINARY format
Definition: AR/profile.h:72
AR::TimeTransformedSymplecticIntegrator::getEkin
Float getEkin() const
Get current kinetic energy.
Definition: symplectic_integrator.h:3014
COMM::ParticleGroup::restoreParticlePosVel
int restoreParticlePosVel(Float *_bk)
restore member particle position and velocity
Definition: particle_group.h:82
COMM::List
list class to store and manage a group of member
Definition: list.h:19
symplectic_step.h
AR::printReference
void printReference(std::ostream &fout, const int offset=4)
print reference to cite
Definition: symplectic_integrator.h:48
AR::InterruptBinary::status
InterruptStatus status
Definition: interrupt.h:16
to_int
#define to_int(x)
Definition: Float.h:26
COMM::List::resizeNoInitialize
void resizeNoInitialize(const int _n)
increase size without initialization
Definition: list.h:268
AR::TimeTransformedSymplecticManager::readBinary
void readBinary(FILE *_fin, int _version=0)
read class data with BINARY format and initial the array
Definition: symplectic_integrator.h:118
COMM::BinaryTree
Binary tree cell.
Definition: binary_tree.h:492
AR::TimeTransformedSymplecticManager::step_count_max
long long unsigned int step_count_max
slowdown maximum timescale to calculate maximum slowdown factor
Definition: symplectic_integrator.h:69
AR::TimeTransformedSymplecticIntegrator::backupIntData
int backupIntData(Float *_bk)
Backup integration data.
Definition: symplectic_integrator.h:3250
COMM::ListMode::local
@ local
AR::TimeTransformedSymplecticIntegrator::printGroupInfo
void printGroupInfo(const int _type, std::ostream &_fout, const int _width, const Tptcl *_pcm=NULL)
print group information
Definition: symplectic_integrator.h:3333
profile.h
AR::TimeTransformedSymplecticIntegrator::getEtot
Float getEtot() const
Get current total energy from ekin and epot.
Definition: symplectic_integrator.h:3035
COMM::List::setMode
void setMode(const ListMode _mode)
set mode
Definition: list.h:39
force.h
AR::TimeTransformedSymplecticManager::TimeTransformedSymplecticManager
TimeTransformedSymplecticManager()
class to manager kick drift step
Definition: symplectic_integrator.h:76
list.h
AR::TimeTransformedSymplecticIntegrator::getEtotFromBackup
Float getEtotFromBackup(Float *_bk) const
get total energy from backup data (ekin+epot)
Definition: symplectic_integrator.h:3069
AR::SlowDown::initialSlowDownReference
void initialSlowDownReference(const Float _kappa_ref, const Float _timescale_max)
initialize slow-down parameters
Definition: slow_down.h:42
COMM::ParticleGroup::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: particle_group.h:112
AR::TimeTransformedSymplecticIntegrator::getEtotRef
Float getEtotRef() const
Get current total integrated energy.
Definition: symplectic_integrator.h:3028
AR::SlowDown
Slow-down parameter control class.
Definition: slow_down.h:11
AR::TimeTransformedSymplecticIntegrator::getEnergyError
Float getEnergyError() const
get energy error
Definition: symplectic_integrator.h:3054
AR::Profile::step_count_tsyn_sum
UInt64 step_count_tsyn_sum
Definition: AR/profile.h:32
AR::TimeTransformedSymplecticIntegrator::checkParams
bool checkParams()
check whether parameters values are correct
Definition: symplectic_integrator.h:237
AR::TimeTransformedSymplecticIntegrator::initialIntegration
void initialIntegration(const Float _time)
initialization for integration
Definition: symplectic_integrator.h:1546
AR::Profile::clear
void clear()
Definition: AR/profile.h:40
COMM::List::reserveMem
void reserveMem(const int _nmax)
Memory allocation for storing members.
Definition: list.h:48
AR::InterruptBinary::clear
void clear()
Definition: interrupt.h:29
COMM::ListMode::copy
@ copy
COMM::BinaryTree::getMember
Tptcl * getMember(const size_t i) const
get member
Definition: binary_tree.h:940
COMM::ParticleGroup::readBinary
void readBinary(FILE *_fin)
Read particle data from file.
Definition: particle_group.h:152
AR::TimeTransformedSymplecticManager::time_error_max
Float time_error_max
Definition: symplectic_integrator.h:60
ISINF
#define ISINF(x)
Definition: Float.h:51
AR::SlowDown::printColumn
void printColumn(std::ostream &_fout, const int _width=20)
print data of class members using column style
Definition: slow_down.h:204
AR::TimeTransformedSymplecticManager::interrupt_detection_option
int interrupt_detection_option
maximum step counts
Definition: symplectic_integrator.h:70
AR::TimeTransformedSymplecticManager
Time Transformed Symplectic integrator manager.
Definition: symplectic_integrator.h:58
AR::SymplecticStep::getOrder
int getOrder() const
get Symplectic order
Definition: symplectic_step.h:330