PeTar
N-body code for collisional gravitational systems
hermite_interaction.hpp
Go to the documentation of this file.
1 #pragma once
2 
3 #include "Common/Float.h"
4 #include "changeover.hpp"
5 
7 class HermiteInteraction{
8 public:
9  Float eps_sq; // softening parameter
10  Float gravitational_constant; // gravitational constant
11 
12  // constructor
13  HermiteInteraction(): eps_sq(Float(-1.0)), gravitational_constant(Float(-1.0)) {}
14 
16 
18  bool checkParams() {
19  ASSERT(eps_sq>=0.0);
20  ASSERT(gravitational_constant>0.0);
21  return true;
22  }
23 
25  void print(std::ostream & _fout) const{
26  _fout<<"eps_sq: "<<eps_sq<<std::endl
27  <<"G : "<<gravitational_constant<<std::endl;
28  }
29 
31 
36  template<class Tpi, class Tpj>
37  inline Float calcR2Pair(const Tpi& _pi,
38  const Tpj& _pj) {
39  const Float dr[3] = {_pj.pos[0]-_pi.pos[0],
40  _pj.pos[1]-_pi.pos[1],
41  _pj.pos[2]-_pi.pos[2]};
42  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
43  return dr2;
44  }
45 
47 
53  template<class Tpi, class Tpj>
54  inline Float calcAccJerkPairSingleSingle(H4::ForceH4& _fi,
55  const Tpi& _pi,
56  const Tpj& _pj) {
57  const Float dr[3] = {_pj.pos[0]-_pi.pos[0],
58  _pj.pos[1]-_pi.pos[1],
59  _pj.pos[2]-_pi.pos[2]};
60  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
61  Float dr2_eps = dr2 + eps_sq;
62  const Float dv[3] = {_pj.vel[0] - _pi.vel[0],
63  _pj.vel[1] - _pi.vel[1],
64  _pj.vel[2] - _pi.vel[2]};
65  const Float drdv = dr[0]*dv[0] + dr[1]*dv[1] + dr[2]*dv[2];
66  const Float r = sqrt(dr2_eps);
67  ASSERT(r>0.0);
68  const Float rinv = 1.0/r;
69  const Float drdot = drdv*rinv;
70  const Float kp = ChangeOver::calcPotWTwo(_pi.changeover,_pj.changeover, r);
71  const Float k = ChangeOver::calcAcc0WTwo(_pi.changeover, _pj.changeover, r);
72  const Float kdot = ChangeOver::calcAcc1WTwo(_pi.changeover, _pj.changeover, r, drdot);
73 
74  const Float rinv2 = rinv*rinv;
75 
76  const Float gmor = gravitational_constant*_pj.mass*rinv;
77  const Float gmor3 = gmor*rinv2;
78  const Float gmor3k = gmor3*k;
79  const Float gmor3kd = gmor3*kdot;
80  const Float acc0[3] = {gmor3k*dr[0], gmor3k*dr[1], gmor3k*dr[2]};
81  const Float acc1[3] = {gmor3k*dv[0] - 3.0*drdv*rinv2*acc0[0] + gmor3kd*dr[0],
82  gmor3k*dv[1] - 3.0*drdv*rinv2*acc0[1] + gmor3kd*dr[1],
83  gmor3k*dv[2] - 3.0*drdv*rinv2*acc0[2] + gmor3kd*dr[2]};
84  _fi.acc0[0] += acc0[0];
85  _fi.acc0[1] += acc0[1];
86  _fi.acc0[2] += acc0[2];
87 
88  _fi.acc1[0] += acc1[0];
89  _fi.acc1[1] += acc1[1];
90  _fi.acc1[2] += acc1[2];
91 
92  _fi.pot += - gmor*kp;
93 
94  return dr2;
95  }
96 
98 
104  template<class Tpi, class Tgroup>
105  inline Float calcAccJerkPairSingleGroupMember(H4::ForceH4& _fi,
106  const Tpi& _pi,
107  const Tgroup& _gj) {
108  const int n_member = _gj.particles.getSize();
109  auto* member_adr = _gj.particles.getOriginAddressArray();
110  Float r2_min = NUMERIC_FLOAT_MAX;
111  for (int i=0; i<n_member; i++) {
112  const auto& pj = *member_adr[i];
113  ASSERT(_pi.id!=pj.id);
114  const Float dr[3] = {pj.pos[0]-_pi.pos[0],
115  pj.pos[1]-_pi.pos[1],
116  pj.pos[2]-_pi.pos[2]};
117  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
118  Float dr2_eps = dr2 + eps_sq;
119  const Float dv[3] = {pj.vel[0] - _pi.vel[0],
120  pj.vel[1] - _pi.vel[1],
121  pj.vel[2] - _pi.vel[2]};
122  const Float drdv = dr[0]*dv[0] + dr[1]*dv[1] + dr[2]*dv[2];
123  const Float r = sqrt(dr2_eps);
124  ASSERT(r>0.0);
125  const Float rinv = 1.0/r;
126  const Float drdot = drdv*rinv;
127  const Float kp = ChangeOver::calcPotWTwo(_pi.changeover, pj.changeover, r);
128  const Float k = ChangeOver::calcAcc0WTwo(_pi.changeover, pj.changeover, r);
129  const Float kdot = ChangeOver::calcAcc1WTwo(_pi.changeover, pj.changeover, r, drdot);
130 
131  const Float rinv2 = rinv*rinv;
132 
133  const Float gmor = gravitational_constant*pj.mass*rinv;
134  const Float gmor3 = gmor*rinv2;
135  const Float gmor3k = gmor3*k;
136  const Float gmor3kd = gmor3*kdot;
137  const Float acc0[3] = {gmor3k*dr[0], gmor3k*dr[1], gmor3k*dr[2]};
138  const Float acc1[3] = {gmor3k*dv[0] - 3.0*drdv*rinv2*acc0[0] + gmor3kd*dr[0],
139  gmor3k*dv[1] - 3.0*drdv*rinv2*acc0[1] + gmor3kd*dr[1],
140  gmor3k*dv[2] - 3.0*drdv*rinv2*acc0[2] + gmor3kd*dr[2]};
141  _fi.acc0[0] += acc0[0];
142  _fi.acc0[1] += acc0[1];
143  _fi.acc0[2] += acc0[2];
144 
145  _fi.acc1[0] += acc1[0];
146  _fi.acc1[1] += acc1[1];
147  _fi.acc1[2] += acc1[2];
148 
149  _fi.pot += - gmor*kp;
150 
151  r2_min = std::min(r2_min, dr2);
152  }
153  return r2_min;
154  }
155 
157 
164  template<class Tpi, class Tgroup, class Tpcmj>
165  inline Float calcAccJerkPairSingleGroupCM(H4::ForceH4& _fi,
166  const Tpi& _pi,
167  const Tgroup& _gj,
168  const Tpcmj& _pj) {
169  const Float dr[3] = {_pj.pos[0]-_pi.pos[0],
170  _pj.pos[1]-_pi.pos[1],
171  _pj.pos[2]-_pi.pos[2]};
172  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
173  Float dr2_eps = dr2 + eps_sq;
174  const Float dv[3] = {_pj.vel[0] - _pi.vel[0],
175  _pj.vel[1] - _pi.vel[1],
176  _pj.vel[2] - _pi.vel[2]};
177  const Float drdv = dr[0]*dv[0] + dr[1]*dv[1] + dr[2]*dv[2];
178  const Float r = sqrt(dr2_eps);
179  ASSERT(r>0.0);
180  const Float rinv = 1.0/r;
181  const Float drdot = drdv*rinv;
182 
183  Float mkp = 0.0;
184  Float mk = 0.0;
185  Float mkdot = 0.0;
186  auto* ptcl_mem = _gj.particles.getDataAddress();
187  ASSERT(_gj.particles.cm.mass==_pj.mass);
188  for (int i=0; i<_gj.particles.getSize(); i++) {
189  mkp += ptcl_mem[i].mass * ChangeOver::calcPotWTwo(_pi.changeover, ptcl_mem[i].changeover, r);
190  mk += ptcl_mem[i].mass * ChangeOver::calcAcc0WTwo(_pi.changeover, ptcl_mem[i].changeover, r);
191  mkdot += ptcl_mem[i].mass * ChangeOver::calcAcc1WTwo(_pi.changeover, ptcl_mem[i].changeover, r, drdot);
192  }
193  const Float rinv2 = rinv*rinv;
194  const Float rinv3 = rinv2*rinv;
195 
196  const Float gmor3k = gravitational_constant*mk*rinv3;
197  const Float gmor3kd = gravitational_constant*mkdot*rinv3;
198  const Float acc0[3] = {gmor3k*dr[0], gmor3k*dr[1], gmor3k*dr[2]};
199  const Float acc1[3] = {gmor3k*dv[0] - 3.0*drdv*rinv2*acc0[0] + gmor3kd*dr[0],
200  gmor3k*dv[1] - 3.0*drdv*rinv2*acc0[1] + gmor3kd*dr[1],
201  gmor3k*dv[2] - 3.0*drdv*rinv2*acc0[2] + gmor3kd*dr[2]};
202  _fi.acc0[0] += acc0[0];
203  _fi.acc0[1] += acc0[1];
204  _fi.acc0[2] += acc0[2];
205 
206  _fi.acc1[0] += acc1[0];
207  _fi.acc1[1] += acc1[1];
208  _fi.acc1[2] += acc1[2];
209 
210  _fi.pot += - gravitational_constant*mkp*rinv;
211 
212  return dr2;
213  }
214 
216 
223  template<class Tgroup, class Tpcmi, class Tpj>
224  inline Float calcAccJerkPairGroupCMSingle(H4::ForceH4& _fi,
225  const Tgroup& _gi,
226  const Tpcmi& _pi,
227  const Tpj& _pj) {
228  const Float dr[3] = {_pj.pos[0]-_pi.pos[0],
229  _pj.pos[1]-_pi.pos[1],
230  _pj.pos[2]-_pi.pos[2]};
231  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
232  Float dr2_eps = dr2 + eps_sq;
233  const Float dv[3] = {_pj.vel[0] - _pi.vel[0],
234  _pj.vel[1] - _pi.vel[1],
235  _pj.vel[2] - _pi.vel[2]};
236  const Float drdv = dr[0]*dv[0] + dr[1]*dv[1] + dr[2]*dv[2];
237  const Float r = sqrt(dr2_eps);
238  ASSERT(r>0.0);
239  const Float rinv = 1.0/r;
240  const Float drdot = drdv*rinv;
241 
242  Float kp = 0.0;
243  Float k = 0.0;
244  Float kdot = 0.0;
245  auto* ptcl_mem = _gi.particles.getDataAddress();
246  ASSERT(_gi.particles.cm.mass==_pi.mass);
247  for (int i=0; i<_gi.particles.getSize(); i++) {
248  kp += ptcl_mem[i].mass * ChangeOver::calcPotWTwo(_pi.changeover, ptcl_mem[i].changeover, r);
249  k += ptcl_mem[i].mass * ChangeOver::calcAcc0WTwo(_pj.changeover, ptcl_mem[i].changeover, r);
250  kdot += ptcl_mem[i].mass * ChangeOver::calcAcc1WTwo(_pj.changeover, ptcl_mem[i].changeover, r, drdot);
251  }
252  kp /= _pi.mass;
253  k /= _pi.mass;
254  kdot /= _pi.mass;
255 
256  const Float rinv2 = rinv*rinv;
257 
258  const Float gmor = gravitational_constant*_pj.mass*rinv;
259  const Float gmor3 = gmor*rinv2;
260  const Float gmor3k = gmor3*k;
261  const Float gmor3kd = gmor3*kdot;
262  const Float acc0[3] = {gmor3k*dr[0], gmor3k*dr[1], gmor3k*dr[2]};
263  const Float acc1[3] = {gmor3k*dv[0] - 3.0*drdv*rinv2*acc0[0] + gmor3kd*dr[0],
264  gmor3k*dv[1] - 3.0*drdv*rinv2*acc0[1] + gmor3kd*dr[1],
265  gmor3k*dv[2] - 3.0*drdv*rinv2*acc0[2] + gmor3kd*dr[2]};
266  _fi.acc0[0] += acc0[0];
267  _fi.acc0[1] += acc0[1];
268  _fi.acc0[2] += acc0[2];
269 
270  _fi.acc1[0] += acc1[0];
271  _fi.acc1[1] += acc1[1];
272  _fi.acc1[2] += acc1[2];
273 
274  _fi.pot += - gmor*kp;
275 
276  return dr2;
277  }
278 
279 
281 
288  template<class Tpi, class Tgroup>
289  inline Float calcAccJerkPairGroupCMGroupMember(H4::ForceH4& _fi,
290  const Tgroup& _gi,
291  const Tpi& _pi,
292  const Tgroup& _gj) {
293  const int n_member = _gj.particles.getSize();
294  auto* member_adr = _gj.particles.getOriginAddressArray();
295  Float r2_min = NUMERIC_FLOAT_MAX;
296  ASSERT(_gi.particles.cm.mass==_pi.mass);
297  for (int i=0; i<n_member; i++) {
298  const auto& pj = *member_adr[i];
299  ASSERT(_pi.id!=pj.id);
300  const Float dr[3] = {pj.pos[0]-_pi.pos[0],
301  pj.pos[1]-_pi.pos[1],
302  pj.pos[2]-_pi.pos[2]};
303  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
304  Float dr2_eps = dr2 + eps_sq;
305  const Float dv[3] = {pj.vel[0] - _pi.vel[0],
306  pj.vel[1] - _pi.vel[1],
307  pj.vel[2] - _pi.vel[2]};
308  const Float drdv = dr[0]*dv[0] + dr[1]*dv[1] + dr[2]*dv[2];
309  const Float r = sqrt(dr2_eps);
310  ASSERT(r>0.0);
311  const Float rinv = 1.0/r;
312  const Float drdot = drdv*rinv;
313 
314  Float k = 0.0;
315  Float kdot = 0.0;
316  Float kp = 0.0;
317  auto* ptcl_mem = _gi.particles.getDataAddress();
318  for (int i=0; i<_gi.particles.getSize(); i++) {
319  kp += ptcl_mem[i].mass * ChangeOver::calcPotWTwo(_pi.changeover, ptcl_mem[i].changeover, r);
320  k += ptcl_mem[i].mass * ChangeOver::calcAcc0WTwo(pj.changeover, ptcl_mem[i].changeover, r);
321  kdot += ptcl_mem[i].mass * ChangeOver::calcAcc1WTwo(pj.changeover, ptcl_mem[i].changeover, r, drdot);
322  }
323  kp /= _pi.mass;
324  k /= _pi.mass;
325  kdot /= _pi.mass;
326 
327  const Float rinv2 = rinv*rinv;
328 
329  const Float gmor = gravitational_constant*pj.mass*rinv;
330  const Float gmor3 = gmor*rinv2;
331  const Float gmor3k = gmor3*k;
332  const Float gmor3kd = gmor3*kdot;
333  const Float acc0[3] = {gmor3k*dr[0], gmor3k*dr[1], gmor3k*dr[2]};
334  const Float acc1[3] = {gmor3k*dv[0] - 3.0*drdv*rinv2*acc0[0] + gmor3kd*dr[0],
335  gmor3k*dv[1] - 3.0*drdv*rinv2*acc0[1] + gmor3kd*dr[1],
336  gmor3k*dv[2] - 3.0*drdv*rinv2*acc0[2] + gmor3kd*dr[2]};
337  _fi.acc0[0] += acc0[0];
338  _fi.acc0[1] += acc0[1];
339  _fi.acc0[2] += acc0[2];
340 
341  _fi.acc1[0] += acc1[0];
342  _fi.acc1[1] += acc1[1];
343  _fi.acc1[2] += acc1[2];
344 
345  _fi.pot += - gmor*kp;
346 
347  r2_min = std::min(r2_min, dr2);
348  }
349  return r2_min;
350  }
351 
353 
361  template<class Tpcmi, class Tgroup, class Tpcmj>
362  inline Float calcAccJerkPairGroupCMGroupCM(H4::ForceH4& _fi,
363  const Tgroup& _gi,
364  const Tpcmi& _pi,
365  const Tgroup& _gj,
366  const Tpcmj& _pj) {
367  const Float dr[3] = {_pj.pos[0]-_pi.pos[0],
368  _pj.pos[1]-_pi.pos[1],
369  _pj.pos[2]-_pi.pos[2]};
370  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
371  Float dr2_eps = dr2 + eps_sq;
372  const Float dv[3] = {_pj.vel[0] - _pi.vel[0],
373  _pj.vel[1] - _pi.vel[1],
374  _pj.vel[2] - _pi.vel[2]};
375  const Float drdv = dr[0]*dv[0] + dr[1]*dv[1] + dr[2]*dv[2];
376  const Float r = sqrt(dr2_eps);
377  ASSERT(r>0.0);
378  const Float rinv = 1.0/r;
379  const Float drdot = drdv*rinv;
380 
381  // pj.mass * k
382  Float mk = 0.0;
383  Float mkdot = 0.0;
384  Float mkp = 0.0;
385  auto* ptcl_mem_i = _gi.particles.getDataAddress();
386  auto* ptcl_mem_j = _gj.particles.getDataAddress();
387  for (int i=0; i<_gi.particles.getSize(); i++) {
388  Float mkj = 0.0;
389  Float mkpj = 0.0;
390  Float mkdotj = 0.0;
391  for (int j=0; j<_gj.particles.getSize(); j++) {
392  mkpj += ptcl_mem_j[j].mass * ChangeOver::calcPotWTwo(ptcl_mem_i[i].changeover, ptcl_mem_j[j].changeover, r);
393  mkj += ptcl_mem_j[j].mass * ChangeOver::calcAcc0WTwo(ptcl_mem_i[i].changeover, ptcl_mem_j[j].changeover, r);
394  mkdotj += ptcl_mem_j[j].mass * ChangeOver::calcAcc1WTwo(ptcl_mem_i[i].changeover, ptcl_mem_j[j].changeover, r, drdot);
395  }
396  mkp += ptcl_mem_i[i].mass * mkpj;
397  mk += ptcl_mem_i[i].mass * mkj;
398  mkdot += ptcl_mem_i[i].mass * mkdotj;
399  }
400  mkp /= _pi.mass;
401  mk /= _pi.mass;
402  mkdot /= _pi.mass;
403 
404  const Float rinv2 = rinv*rinv;
405  const Float rinv3 = rinv2*rinv;
406 
407  const Float gmor3k = gravitational_constant*mk*rinv3;
408  const Float gmor3kd = gravitational_constant*mkdot*rinv3;
409  const Float acc0[3] = {gmor3k*dr[0], gmor3k*dr[1], gmor3k*dr[2]};
410  const Float acc1[3] = {gmor3k*dv[0] - 3.0*drdv*rinv2*acc0[0] + gmor3kd*dr[0],
411  gmor3k*dv[1] - 3.0*drdv*rinv2*acc0[1] + gmor3kd*dr[1],
412  gmor3k*dv[2] - 3.0*drdv*rinv2*acc0[2] + gmor3kd*dr[2]};
413  _fi.acc0[0] += acc0[0];
414  _fi.acc0[1] += acc0[1];
415  _fi.acc0[2] += acc0[2];
416 
417  _fi.acc1[0] += acc1[0];
418  _fi.acc1[1] += acc1[1];
419  _fi.acc1[2] += acc1[2];
420 
421  _fi.pot += - gravitational_constant*mkp*rinv;
422 
423  return dr2;
424  }
425 
427  template<class Tpi, class Tpj>
428  Float calcEnergyPotSingleSingle(const Tpi& pi, const Tpj& pj) {
429  const Float dr[3] = {pj.pos[0] - pi.pos[0],
430  pj.pos[1] - pi.pos[1],
431  pj.pos[2] - pi.pos[2]};
432  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
433  Float dr2_eps = dr2 + eps_sq;
434  const Float r = sqrt(dr2_eps);
435  ASSERT(r>0.0);
436  const Float rinv = 1.0/r;
437  const Float k = ChangeOver::calcPotWTwo(pi.changeover, pj.changeover, r);
438 
439  return -gravitational_constant*pi.mass*pj.mass*rinv*k;
440  }
441 
443 
447  template<class Tgroup, class Tpert>
448  Float calcEnergyPertOneGroup(const Tgroup& _group, const Tpert& _perturber) {
449  Float epert = 0.0;
450 #ifdef SOFT_PERT
451  if (_group.perturber.soft_pert!=NULL) {
452  auto* pert = _group.perturber.soft_pert;
453  const int n_member = _group.particles.getSize();
454  for (int j=0; j<n_member; j++) {
455  epert += _group.particles[j].mass*pert->evalPot(_group.particles[j].pos);
456  }
457  }
458 #endif
459  return epert;
460  }
461 
463 
472  template<class Tp, class Tgroup, class Tpert>
473  inline void calcEnergy(H4::HermiteEnergy& _energy, const Tp* _particles, const int _n_particle, const Tgroup* _groups, const int* _group_index, const int _n_group, const Tpert& _perturber) {
474  _energy.ekin = _energy.epot = _energy.epert = 0.0;
475  for (int i=0; i<_n_particle; i++) {
476  auto& pi = _particles[i];
477  if (pi.mass==0.0) continue;
478  _energy.ekin += pi.mass* (pi.vel[0]*pi.vel[0] + pi.vel[1]*pi.vel[1] + pi.vel[2]*pi.vel[2]);
479  Float poti = 0.0;
480  for (int j=0; j<i; j++) {
481  if (i==j) continue;
482  auto& pj = _particles[j];
483  if (pj.mass==0.0) continue;
484  const Float dr[3] = {pj.pos[0] - pi.pos[0],
485  pj.pos[1] - pi.pos[1],
486  pj.pos[2] - pi.pos[2]};
487  Float dr2 = dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2];
488  Float dr2_eps = dr2 + eps_sq;
489  const Float r = sqrt(dr2_eps);
490  ASSERT(r>0.0);
491  const Float rinv = 1.0/r;
492  const Float k = ChangeOver::calcPotWTwo(pi.changeover, pj.changeover, r);
493 
494  poti += -pj.mass*rinv*k;
495  }
496  _energy.epot += gravitational_constant*poti*pi.mass;
497  }
498 
499 #ifdef SOFT_PERT
500  // tidal tensor energy
501  for (int k=0; k<_n_group; k++) {
502  const int i = _group_index[k];
503  _energy.epert += calcEnergyPertOneGroup(_groups[i], _perturber);
504  auto& bink = _groups[i].info.getBinaryTreeRoot();
505  auto& pcm = _groups[i].particles.cm;
506  auto& vcm = pcm.vel;
507  auto& vbin = bink.vel;
508  auto& vbin_bk = _groups[i].info.vcm_record;
509  //Float vbcm[3] = {vcm[0] + vbin_bk[0], vcm[1] + vbin_bk[1], vcm[2] + vbin_bk[2]};
510  Float dvbin[3] = {vbin[0] - vbin_bk[0], vbin[1] - vbin_bk[1], vbin[2] - vbin_bk[2]};
511  Float de_kin = bink.mass*(dvbin[0]*vcm[0]+dvbin[1]*vcm[1]+dvbin[2]*vcm[2]);
512  _energy.epert -= de_kin;
513  }
514 #endif
515  _energy.ekin *= 0.5;
516  //_energy.epert *= 0.5;
517  }
518 
520 
522  void writeBinary(FILE *_fp) const {
523  fwrite(this, sizeof(*this),1,_fp);
524  }
525 
527 
529  void readBinary(FILE *_fin) {
530  size_t rcount = fread(this, sizeof(*this), 1, _fin);
531  if (rcount<1) {
532  std::cerr<<"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<".\n";
533  abort();
534  }
535  }
536 };
HermiteInteraction::calcAccJerkPairGroupCMSingle
Float calcAccJerkPairGroupCMSingle(H4::ForceH4 &_fi, const Tgroup &_gi, const Tpcmi &_pi, const Tpj &_pj)
calculate acceleration and jerk of one pair group cm and single
Definition: hermite_interaction.hpp:224
HermiteInteraction::calcEnergyPotSingleSingle
Float calcEnergyPotSingleSingle(const Tpi &pi, const Tpj &pj)
calculate pair potential energy
Definition: hermite_interaction.hpp:428
HermiteInteraction::checkParams
bool checkParams()
check whether parameters values are correct
Definition: hermite_interaction.hpp:18
HermiteInteraction::calcEnergyPertOneGroup
Float calcEnergyPertOneGroup(const Tgroup &_group, const Tpert &_perturber)
calculate peturbation energy to one group
Definition: hermite_interaction.hpp:448
HermiteInteraction::gravitational_constant
Float gravitational_constant
Definition: hermite_interaction.hpp:10
HermiteInteraction::HermiteInteraction
HermiteInteraction()
Definition: hermite_interaction.hpp:13
HermiteInteraction::print
void print(std::ostream &_fout) const
print parameters
Definition: hermite_interaction.hpp:25
HermiteInteraction::calcAccJerkPairSingleGroupCM
Float calcAccJerkPairSingleGroupCM(H4::ForceH4 &_fi, const Tpi &_pi, const Tgroup &_gj, const Tpcmj &_pj)
calculate acceleration and jerk of one pair single and resolved group
Definition: hermite_interaction.hpp:165
HermiteInteraction::calcEnergy
void calcEnergy(H4::HermiteEnergy &_energy, const Tp *_particles, const int _n_particle, const Tgroup *_groups, const int *_group_index, const int _n_group, const Tpert &_perturber)
calculate kinetic and potential energy of the system
Definition: hermite_interaction.hpp:473
PIKG::min
T min(const Vector3< T > &v)
Definition: pikg_vector.hpp:150
HermiteInteraction::writeBinary
void writeBinary(FILE *_fp) const
write class data to file with binary format
Definition: hermite_interaction.hpp:522
changeover.hpp
ChangeOver::calcPotWTwo
static Float calcPotWTwo(const ChangeOver &_ch1, const ChangeOver &_ch2, const Float &_dr)
calculate changeover function Pot by selecting maximum rout
Definition: changeover.hpp:366
ChangeOver::calcAcc0WTwo
static Float calcAcc0WTwo(const ChangeOver &_ch1, const ChangeOver &_ch2, const Float &_dr)
calculate changeover function Acc0 by selecting maximum rout
Definition: changeover.hpp:372
HermiteInteraction
hermite interaction class
Definition: hermite_interaction.hpp:7
ChangeOver::calcAcc1WTwo
static Float calcAcc1WTwo(const ChangeOver &_ch1, const ChangeOver &_ch2, const Float &_dr, const Float &_drdot)
calculate changeover function Acc1 by selecting maximum rout
Definition: changeover.hpp:378
HermiteInteraction::calcAccJerkPairSingleGroupMember
Float calcAccJerkPairSingleGroupMember(H4::ForceH4 &_fi, const Tpi &_pi, const Tgroup &_gj)
calculate acceleration and jerk of one pair single and resolved group
Definition: hermite_interaction.hpp:105
HermiteInteraction::calcR2Pair
Float calcR2Pair(const Tpi &_pi, const Tpj &_pj)
calculate separation square between i and j particles
Definition: hermite_interaction.hpp:37
HermiteInteraction::calcAccJerkPairGroupCMGroupMember
Float calcAccJerkPairGroupCMGroupMember(H4::ForceH4 &_fi, const Tgroup &_gi, const Tpi &_pi, const Tgroup &_gj)
calculate acceleration and jerk of one pair group cm and resolved group
Definition: hermite_interaction.hpp:289
ASSERT
#define ASSERT(expr)
Definition: hard_assert.hpp:238
HermiteInteraction::eps_sq
Float eps_sq
Definition: hermite_interaction.hpp:9
HermiteInteraction::calcAccJerkPairGroupCMGroupCM
Float calcAccJerkPairGroupCMGroupCM(H4::ForceH4 &_fi, const Tgroup &_gi, const Tpcmi &_pi, const Tgroup &_gj, const Tpcmj &_pj)
calculate acceleration and jerk of one pair single and resolved group
Definition: hermite_interaction.hpp:362
HermiteInteraction::calcAccJerkPairSingleSingle
Float calcAccJerkPairSingleSingle(H4::ForceH4 &_fi, const Tpi &_pi, const Tpj &_pj)
calculate acceleration and jerk of one pair
Definition: hermite_interaction.hpp:54
HermiteInteraction::readBinary
void readBinary(FILE *_fin)
read class data to file with binary format
Definition: hermite_interaction.hpp:529