/**
* @file integrator_ias15.c
* @brief IAS15 integrator.
* @author Hanno Rein <hanno@hanno-rein.de>
* @details This file implements the IAS15 integration scheme.
* IAS stands for Integrator with Adaptive Step-size control, 15th
* order. This scheme is a fifteenth order integrator well suited for
* high accuracy orbit integration with non-conservative forces.
* For more details see Rein & Spiegel 2014. Also see Everhart, 1985,
* ASSL Vol. 115, IAU Colloq. 83, Dynamics of Comets, Their Origin
* and Evolution, 185 for the original implementation by Everhart.
* Part of this code is based a function from the ORSE package.
* See orsa.sourceforge.net for more details on their implementation.
*
*
* @section LICENSE
* Copyright (c) 2011-2012 Hanno Rein, Dave Spiegel.
* Copyright (c) 2002-2004 Pasquale Tricarico.
*
* This file is part of rebound.
*
* rebound is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* rebound is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with rebound. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include "rebound.h"
#include "rebound_internal.h"
#include <string.h>
#include <math.h>
// Uncomment the following line to generate numerical constants with extended precision.
//#define GENERATE_CONSTANTS
#include "particle.h"
#include "gravity.h"
#include "tools.h"
#include "integrator_ias15.h"
#include "binarydata.h"
void* reb_integrator_ias15_create();
void reb_integrator_ias15_free(void* p);
void reb_integrator_ias15_step(struct reb_simulation* const r, void* state);
const struct reb_binarydata_field_descriptor reb_integrator_ias15_field_descriptor_list[];
const struct reb_integrator reb_integrator_ias15 = {
.documentation =
"IAS15 is the default integrator of REBOUND."
"IAS15 stands for Integrator with Adaptive Step-size control, 15th order. "
"It is a very high order, non-symplectic integrator which can handle arbitrary "
"forces (including those who are velocity dependent). "
"It is in most cases accurate down to machine precision (16 significant decimal digits). "
"The IAS15 implementation in REBOUND can integrate variational equations. "
"The algorithm is described in detail in [Rein & Spiegel (2015)] and also in the "
"original paper by [Everhart (1985)]. "
"The timestepping choice is described in [Pham et al. (2024)]. "
"A short [YouTube video] also provides a basic introduction to this integrator."
"\n\n"
"[Rein & Spiegel (2015)]: https://ui.adsabs.harvard.edu/abs/2015MNRAS.446.1424R/abstract\n"
"[Pham et al. (2024)]: https://ui.adsabs.harvard.edu/abs/2024OJAp....7E...1P/abstract\n"
"[Everhart (1985)]: https://ui.adsabs.harvard.edu/abs/1985ASSL..115..185E/abstract\n"
"[Aarseth (1985)]: https://ui.adsabs.harvard.edu/abs/1985IAUS..113..251A/abstract\n"
"[YouTube video]: https://www.youtube.com/watch?v=UILEgdZt-fw\n"
,
.step = reb_integrator_ias15_step,
.create = reb_integrator_ias15_create,
.free = reb_integrator_ias15_free,
.field_descriptor_list = reb_integrator_ias15_field_descriptor_list,
};
const struct reb_binarydata_field_descriptor reb_integrator_ias15_field_descriptor_list[] = {
{ "IAS15 is an adaptive integrator. It chooses its timesteps automatically. "
"This parameter controls the accuracy of the integrator. The default "
"value is 1e-9. Setting this parameter to 0 turns off adaptive timestepping "
"and a constant timestep will is used. Turning off adaptive time-stepping is "
"rarely useful. "
"\n\n"
"It is tempting to change this parameter to achieve a speedup at the loss of "
"some accuracy. However, that makes rarely sense. The reason is that IAS15 "
"is a very high (15th!) order integrator. Suppose we increase the timestep "
"by a factor of 10. This will increase the error by a factor of 1e15. "
"In other words, a simulation that previously was converged to machine "
"precision will now have an error of order unity. ",
REB_DOUBLE, "epsilon", offsetof(struct reb_integrator_ias15_state, epsilon), 0, 0, 0},
{ "This sets the minimum allowed timestep. The default value is 0. Set this "
"to a finite value if the adaptive timestep becomes excessively small, for "
"example during close encounters or because of finite floating point precision. "
"Use with caution and make sure the simulation results still make physical sense "
"as you might be in danger of ignoring small timescales in the problem. ",
REB_DOUBLE, "min_dt", offsetof(struct reb_integrator_ias15_state, min_dt), 0, 0, 0},
{ "This flag determines how the adaptive timestep is chosen. The previous name "
"of this flag was `epsilon_global`. The default is PRS23 which corresponds to "
"the timestep criterion described in [Pham et al. (2024)]. This default "
"should be optimal in almost all cases. "
"If set to INDIVIDUAL, the fractional error is estimated via `max(acceleration_error/acceleration)` "
"and the timestep criterion of [Rein & Spiegel (2015)] is used. "
"If set to GLOBAL, IAS15 estimates the fractional error via `max(acceleration_error)/max(acceleration)` "
"where the maximum is taken over all particles and the timestep criterion of "
"[Rein & Spiegel (2015)] is used. This was the default until January 2024. "
"If set to AARSETH85, then the criterion of [Aarseth (1985)]) is used.",
REB_UINT, "adaptive_mode", offsetof(struct reb_integrator_ias15_state, adaptive_mode), 0, 0, REB_GENERATE_ENUM_DESCRIPTORS(REB_IAS15_ADAPTIVEMODE)},
{ "", REB_UINT64, "iterations_max_exceeded", offsetof(struct reb_integrator_ias15_state, iterations_max_exceeded), 0, 0, 0},
{ "", REB_POINTER, "at", offsetof(struct reb_integrator_ias15_state, at), offsetof(struct reb_integrator_ias15_state, N_allocated), sizeof(double), 0},
{ "", REB_POINTER, "x0", offsetof(struct reb_integrator_ias15_state, x0), offsetof(struct reb_integrator_ias15_state, N_allocated), sizeof(double), 0},
{ "", REB_POINTER, "v0", offsetof(struct reb_integrator_ias15_state, v0), offsetof(struct reb_integrator_ias15_state, N_allocated), sizeof(double), 0},
{ "", REB_POINTER, "a0", offsetof(struct reb_integrator_ias15_state, a0), offsetof(struct reb_integrator_ias15_state, N_allocated), sizeof(double), 0},
{ "", REB_POINTER, "csx", offsetof(struct reb_integrator_ias15_state, csx), offsetof(struct reb_integrator_ias15_state, N_allocated), sizeof(double), 0},
{ "", REB_POINTER, "csv", offsetof(struct reb_integrator_ias15_state, csv), offsetof(struct reb_integrator_ias15_state, N_allocated), sizeof(double), 0},
{ "", REB_POINTER, "csa0", offsetof(struct reb_integrator_ias15_state, csa0), offsetof(struct reb_integrator_ias15_state, N_allocated), sizeof(double), 0},
{ "", REB_POINTER, "g", offsetof(struct reb_integrator_ias15_state, g), offsetof(struct reb_integrator_ias15_state, N_allocated), 7*sizeof(double), 0},
{ "", REB_POINTER, "b", offsetof(struct reb_integrator_ias15_state, b), offsetof(struct reb_integrator_ias15_state, N_allocated), 7*sizeof(double), 0},
{ "", REB_POINTER, "csb", offsetof(struct reb_integrator_ias15_state, csb), offsetof(struct reb_integrator_ias15_state, N_allocated), 7*sizeof(double), 0},
{ "", REB_POINTER, "e", offsetof(struct reb_integrator_ias15_state, e), offsetof(struct reb_integrator_ias15_state, N_allocated), 7*sizeof(double), 0},
{ "", REB_POINTER, "br", offsetof(struct reb_integrator_ias15_state, br), offsetof(struct reb_integrator_ias15_state, N_allocated), 7*sizeof(double), 0},
{ "", REB_POINTER, "er", offsetof(struct reb_integrator_ias15_state, er), offsetof(struct reb_integrator_ias15_state, N_allocated), 7*sizeof(double), 0},
{ 0 }, // Null terminated list
};
void* reb_integrator_ias15_create(){
struct reb_integrator_ias15_state* ias15 = calloc(sizeof(struct reb_integrator_ias15_state),1);
ias15->epsilon = 1e-9;
ias15->min_dt = 0.0;
ias15->adaptive_mode = REB_IAS15_ADAPTIVEMODE_PRS23; // new default since January 2024
return ias15;
}
void reb_integrator_ias15_free(void* p){
struct reb_integrator_ias15_state* ias15 = p;
free(ias15->g);
free(ias15->e);
free(ias15->b);
free(ias15->csb);
free(ias15->er);
free(ias15->br);
free(ias15->at);
free(ias15->x0);
free(ias15->v0);
free(ias15->a0);
free(ias15->csx);
free(ias15->csv);
free(ias15->csa0);
free(ias15);
}
// Generic pointer with 7 elements, for internal use only (IAS15).
struct reb_dp7 {
double* REB_RESTRICT p0;
double* REB_RESTRICT p1;
double* REB_RESTRICT p2;
double* REB_RESTRICT p3;
double* REB_RESTRICT p4;
double* REB_RESTRICT p5;
double* REB_RESTRICT p6;
};
struct reb_dpconst7 {
double* const REB_RESTRICT p0;
double* const REB_RESTRICT p1;
double* const REB_RESTRICT p2;
double* const REB_RESTRICT p3;
double* const REB_RESTRICT p4;
double* const REB_RESTRICT p5;
double* const REB_RESTRICT p6;
};
// Helper functions for resetting the b and e coefficients
static void copybuffers(const struct reb_dpconst7 _a, const struct reb_dpconst7 _b, size_t N3);
static void predict_next_step(double ratio, size_t N3, const struct reb_dpconst7 _e, const struct reb_dpconst7 _b, const struct reb_dpconst7 e, const struct reb_dpconst7 b);
/////////////////////////
// Constants
static const double safety_factor = 0.25; /**< Maximum increase/deacrease of consecutve timesteps. */
// Gauss Radau spacings
static const double h[8] = { 0.0, 0.0562625605369221464656521910318, 0.180240691736892364987579942780, 0.352624717113169637373907769648, 0.547153626330555383001448554766, 0.734210177215410531523210605558, 0.885320946839095768090359771030, 0.977520613561287501891174488626};
// Other constants
static const double rr[28] = {0.0562625605369221464656522, 0.1802406917368923649875799, 0.1239781311999702185219278, 0.3526247171131696373739078, 0.2963621565762474909082556, 0.1723840253762772723863278, 0.5471536263305553830014486, 0.4908910657936332365357964, 0.3669129345936630180138686, 0.1945289092173857456275408, 0.7342101772154105315232106, 0.6779476166784883850575584, 0.5539694854785181665356307, 0.3815854601022408941493028, 0.1870565508848551485217621, 0.8853209468390957680903598, 0.8290583863021736216247076, 0.7050802551022034031027798, 0.5326962297259261307164520, 0.3381673205085403850889112, 0.1511107696236852365671492, 0.9775206135612875018911745, 0.9212580530243653554255223, 0.7972799218243951369035945, 0.6248958964481178645172667, 0.4303669872307321188897259, 0.2433104363458769703679639, 0.0921996667221917338008147};
static const double c[21] = {-0.0562625605369221464656522, 0.0101408028300636299864818, -0.2365032522738145114532321, -0.0035758977292516175949345, 0.0935376952594620658957485, -0.5891279693869841488271399, 0.0019565654099472210769006, -0.0547553868890686864408084, 0.4158812000823068616886219, -1.1362815957175395318285885, -0.0014365302363708915424460, 0.0421585277212687077072973, -0.3600995965020568122897665, 1.2501507118406910258505441, -1.8704917729329500633517991, 0.0012717903090268677492943, -0.0387603579159067703699046, 0.3609622434528459832253398, -1.4668842084004269643701553, 2.9061362593084293014237913, -2.7558127197720458314421588};
static const double d[21] = {0.0562625605369221464656522, 0.0031654757181708292499905, 0.2365032522738145114532321, 0.0001780977692217433881125, 0.0457929855060279188954539, 0.5891279693869841488271399, 0.0000100202365223291272096, 0.0084318571535257015445000, 0.2535340690545692665214616, 1.1362815957175395318285885, 0.0000005637641639318207610, 0.0015297840025004658189490, 0.0978342365324440053653648, 0.8752546646840910912297246, 1.8704917729329500633517991, 0.0000000317188154017613665, 0.0002762930909826476593130, 0.0360285539837364596003871, 0.5767330002770787313544596, 2.2485887607691597933926895, 2.7558127197720458314421588};
// Weights for integration of a first order differential equation (Note: interval length = 2)
static const double w[8] = {0.03125, 0.185358154802979278540728972807180754479812609, 0.304130620646785128975743291458180383736715043, 0.376517545389118556572129261157225608762708603, 0.391572167452493593082499533303669362149363727, 0.347014795634501068709955597003528601733139176, 0.249647901329864963257869294715235590174262844, 0.114508814744257199342353731044292225247093225};
// Machine independent implementation of pow(*,1./7.)
static double sqrt7(double a){
// Without scaling, this is only accurate for arguments in [1e-7, 1e2]
// With scaling: [1e-14, 1e8]
double scale = 1;
while(a<1e-7 && isnormal(a)){
scale *= 0.1;
a *= 1e7;
}
while(a>1e2 && isnormal(a)){
scale *= 10;
a *= 1e-7;
}
double x = 1.;
for (size_t k=0; k<20;k++){ // A smaller number should be ok too.
double x6 = x*x*x*x*x*x;
x += (a/x6-x)/7.;
}
return x*scale;
}
// Casting to reb_dp7 datatype for more convenient access
static struct reb_dpconst7 dpcast(double* dp, size_t N3){
struct reb_dpconst7 dpc = {
.p0 = dp,
.p1 = dp+1*N3,
.p2 = dp+2*N3,
.p3 = dp+3*N3,
.p4 = dp+4*N3,
.p5 = dp+5*N3,
.p6 = dp+6*N3,
};
return dpc;
}
static void realloc_dp7(double* REB_RESTRICT * p, const size_t N3){
*p = realloc(*p, sizeof(double)*N3*7);
for (size_t i=0;i<N3*7;i++){
(*p)[i] = 0.0;
}
}
static inline void add_cs(double* p, double* csp, double inp){
const double y = inp - *csp;
const double t = *p + y;
*csp = (t - *p) - y;
*p = t;
}
static void reb_integrator_ias15_alloc(struct reb_simulation* r, struct reb_integrator_ias15_state* ias15){
size_t N;
if (r->map){
N = r->N_map;
}else{
N = r->N + r->N_var;
}
size_t N3 = 3*N;
if (N3 > ias15->N_allocated) {
realloc_dp7(&ias15->g, N3);
realloc_dp7(&ias15->e, N3);
realloc_dp7(&ias15->b, N3);
realloc_dp7(&ias15->csb, N3);
realloc_dp7(&ias15->er, N3);
realloc_dp7(&ias15->br, N3);
ias15->at = realloc(ias15->at, sizeof(double)*N3);
ias15->x0 = realloc(ias15->x0, sizeof(double)*N3);
ias15->v0 = realloc(ias15->v0, sizeof(double)*N3);
ias15->a0 = realloc(ias15->a0, sizeof(double)*N3);
ias15->csx= realloc(ias15->csx, sizeof(double)*N3);
ias15->csv= realloc(ias15->csv, sizeof(double)*N3);
ias15->csa0 = realloc(ias15->csa0,sizeof(double)*N3);
double* REB_RESTRICT const csx = ias15->csx;
double* REB_RESTRICT const csv = ias15->csv;
for (size_t i=0;i<N3;i++){
// Kill compensated summation coefficients
csx[i] = 0;
csv[i] = 0;
}
ias15->N_allocated = N3;
}
}
// Tries to do one actual timestep.
// Returns 1 if successful.
// Returns 0 if step is rejected.
static int reb_integrator_ias15_step_try(struct reb_simulation* r, struct reb_integrator_ias15_state* ias15) {
reb_integrator_ias15_alloc(r, ias15);
struct reb_particle* const particles = r->particles;
struct reb_particle* const particles_var = r->particles_var;
size_t N; // Number of real particles
size_t N_var; // Number of variational particles
size_t* map = r->map; // this map allow for integrating only a selection of particles
if (map){
N = r->N_map;
N_var = 0;
}else{
N = r->N;
N_var = r->N_var;
}
const size_t N3 = 3*(N+N_var);
reb_simulation_update_acceleration(r);
double* REB_RESTRICT const csx = ias15->csx;
double* REB_RESTRICT const csv = ias15->csv;
double* REB_RESTRICT const csa0 = ias15->csa0;
double* REB_RESTRICT const at = ias15->at;
double* REB_RESTRICT const x0 = ias15->x0;
double* REB_RESTRICT const v0 = ias15->v0;
double* REB_RESTRICT const a0 = ias15->a0;
struct reb_vec3d* gravity_cs = r->gravity_cs;
const struct reb_dpconst7 g = dpcast(ias15->g, N3);
const struct reb_dpconst7 e = dpcast(ias15->e, N3);
const struct reb_dpconst7 b = dpcast(ias15->b, N3);
const struct reb_dpconst7 csb= dpcast(ias15->csb, N3);
const struct reb_dpconst7 er = dpcast(ias15->er, N3);
const struct reb_dpconst7 br = dpcast(ias15->br, N3);
for(size_t k=0;k<N;k++) {
size_t mk = map ? map[k] : k;
x0[3*k] = particles[mk].x;
x0[3*k+1] = particles[mk].y;
x0[3*k+2] = particles[mk].z;
v0[3*k] = particles[mk].vx;
v0[3*k+1] = particles[mk].vy;
v0[3*k+2] = particles[mk].vz;
a0[3*k] = particles[mk].ax;
a0[3*k+1] = particles[mk].ay;
a0[3*k+2] = particles[mk].az;
}
for(size_t mk=0;mk<N_var;mk++) {
size_t k = mk + N;
x0[3*k] = particles_var[mk].x;
x0[3*k+1] = particles_var[mk].y;
x0[3*k+2] = particles_var[mk].z;
v0[3*k] = particles_var[mk].vx;
v0[3*k+1] = particles_var[mk].vy;
v0[3*k+2] = particles_var[mk].vz;
a0[3*k] = particles_var[mk].ax;
a0[3*k+1] = particles_var[mk].ay;
a0[3*k+2] = particles_var[mk].az;
}
if (r->gravity==REB_GRAVITY_COMPENSATED){
for(size_t k=0;k<N+N_var;k++) {
size_t mk = map ? map[k] : k;
csa0[3*k] = gravity_cs[mk].x;
csa0[3*k+1] = gravity_cs[mk].y;
csa0[3*k+2] = gravity_cs[mk].z;
}
}else{
gravity_cs = (struct reb_vec3d*)csa0; // Always 0.
for(size_t k=0;k<N3;k++) {
csa0[k] = 0;
}
}
for (size_t k=0;k<N3;k++){
// Memset might be faster!
csb.p0[k] = 0.;
csb.p1[k] = 0.;
csb.p2[k] = 0.;
csb.p3[k] = 0.;
csb.p4[k] = 0.;
csb.p5[k] = 0.;
csb.p6[k] = 0.;
}
for(size_t k=0;k<N3;k++) {
g.p0[k] = b.p6[k]*d[15] + b.p5[k]*d[10] + b.p4[k]*d[6] + b.p3[k]*d[3] + b.p2[k]*d[1] + b.p1[k]*d[0] + b.p0[k];
g.p1[k] = b.p6[k]*d[16] + b.p5[k]*d[11] + b.p4[k]*d[7] + b.p3[k]*d[4] + b.p2[k]*d[2] + b.p1[k];
g.p2[k] = b.p6[k]*d[17] + b.p5[k]*d[12] + b.p4[k]*d[8] + b.p3[k]*d[5] + b.p2[k];
g.p3[k] = b.p6[k]*d[18] + b.p5[k]*d[13] + b.p4[k]*d[9] + b.p3[k];
g.p4[k] = b.p6[k]*d[19] + b.p5[k]*d[14] + b.p4[k];
g.p5[k] = b.p6[k]*d[20] + b.p5[k];
g.p6[k] = b.p6[k];
}
double integrator_megno_thisdt = 0.;
double integrator_megno_thisdt_init = 0.;
if (r->calculate_megno){
integrator_megno_thisdt_init = w[0]* (r->t-r->megno_initial_t) * reb_tools_megno_deltad_delta(r);
}
double t_beginning = r->t;
double predictor_corrector_error = 1e300;
double predictor_corrector_error_last = 2;
size_t iterations = 0;
// Predictor corrector loop
// Stops if one of the following conditions is satisfied:
// 1) predictor_corrector_error better than 1e-16
// 2) predictor_corrector_error starts to oscillate
// 3) more than 12 iterations
while(1){
if(predictor_corrector_error<1e-16){
break;
}
if(iterations > 2 && predictor_corrector_error_last <= predictor_corrector_error){
break;
}
if (iterations>=12){
ias15->iterations_max_exceeded++;
const size_t integrator_iterations_warning = 10;
if (ias15->iterations_max_exceeded==integrator_iterations_warning ){
reb_simulation_warning(r, "At least 10 predictor corrector loops in IAS15 did not converge. This is typically an indication of the timestep being too large.");
}
break; // Quit predictor corrector loop
}
predictor_corrector_error_last = predictor_corrector_error;
predictor_corrector_error = 0;
iterations++;
integrator_megno_thisdt = integrator_megno_thisdt_init;
for(size_t n=1;n<8;n++) { // Loop over interval using Gauss-Radau spacings
r->t = t_beginning + r->dt * h[n];
// Prepare particles arrays for force calculation
for(size_t i=0;i<N+N_var;i++) { // Predict positions at interval n using b values
const size_t k0 = 3*i+0;
const size_t k1 = 3*i+1;
const size_t k2 = 3*i+2;
double xk0;
double xk1;
double xk2;
xk0 = -csx[k0] + ((((((((b.p6[k0]*7.*h[n]/9. + b.p5[k0])*3.*h[n]/4. + b.p4[k0])*5.*h[n]/7. + b.p3[k0])*2.*h[n]/3. + b.p2[k0])*3.*h[n]/5. + b.p1[k0])*h[n]/2. + b.p0[k0])*h[n]/3. + a0[k0])*r->dt*h[n]/2. + v0[k0])*r->dt*h[n];
xk1 = -csx[k1] + ((((((((b.p6[k1]*7.*h[n]/9. + b.p5[k1])*3.*h[n]/4. + b.p4[k1])*5.*h[n]/7. + b.p3[k1])*2.*h[n]/3. + b.p2[k1])*3.*h[n]/5. + b.p1[k1])*h[n]/2. + b.p0[k1])*h[n]/3. + a0[k1])*r->dt*h[n]/2. + v0[k1])*r->dt*h[n];
xk2 = -csx[k2] + ((((((((b.p6[k2]*7.*h[n]/9. + b.p5[k2])*3.*h[n]/4. + b.p4[k2])*5.*h[n]/7. + b.p3[k2])*2.*h[n]/3. + b.p2[k2])*3.*h[n]/5. + b.p1[k2])*h[n]/2. + b.p0[k2])*h[n]/3. + a0[k2])*r->dt*h[n]/2. + v0[k2])*r->dt*h[n];
if (i<N){
size_t mi = map ? map[i] : i;
particles[mi].x = xk0 + x0[k0];
particles[mi].y = xk1 + x0[k1];
particles[mi].z = xk2 + x0[k2];
}else{
size_t mi = i-N;
particles_var[mi].x = xk0 + x0[k0];
particles_var[mi].y = xk1 + x0[k1];
particles_var[mi].z = xk2 + x0[k2];
}
}
if (r->calculate_megno || (r->additional_forces && r->force_is_velocity_dependent)){
for(size_t i=0;i<N+N_var;i++) { // Predict velocities at interval n using b values
const size_t k0 = 3*i+0;
const size_t k1 = 3*i+1;
const size_t k2 = 3*i+2;
double vk0;
double vk1;
double vk2;
vk0 = -csv[k0] + (((((((b.p6[k0]*7.*h[n]/8. + b.p5[k0])*6.*h[n]/7. + b.p4[k0])*5.*h[n]/6. + b.p3[k0])*4.*h[n]/5. + b.p2[k0])*3.*h[n]/4. + b.p1[k0])*2.*h[n]/3. + b.p0[k0])*h[n]/2. + a0[k0])*r->dt*h[n];
vk1 = -csv[k1] + (((((((b.p6[k1]*7.*h[n]/8. + b.p5[k1])*6.*h[n]/7. + b.p4[k1])*5.*h[n]/6. + b.p3[k1])*4.*h[n]/5. + b.p2[k1])*3.*h[n]/4. + b.p1[k1])*2.*h[n]/3. + b.p0[k1])*h[n]/2. + a0[k1])*r->dt*h[n];
vk2 = -csv[k2] + (((((((b.p6[k2]*7.*h[n]/8. + b.p5[k2])*6.*h[n]/7. + b.p4[k2])*5.*h[n]/6. + b.p3[k2])*4.*h[n]/5. + b.p2[k2])*3.*h[n]/4. + b.p1[k2])*2.*h[n]/3. + b.p0[k2])*h[n]/2. + a0[k2])*r->dt*h[n];
if (i<N){
size_t mi = map ? map[i] : i;
particles[mi].vx = vk0 + v0[k0];
particles[mi].vy = vk1 + v0[k1];
particles[mi].vz = vk2 + v0[k2];
}else{
size_t mi = i-N;
particles_var[mi].vx = vk0 + v0[k0];
particles_var[mi].vy = vk1 + v0[k1];
particles_var[mi].vz = vk2 + v0[k2];
}
}
}
reb_simulation_update_acceleration(r); // Calculate forces at interval n
if (r->calculate_megno){
integrator_megno_thisdt += w[n] * (r->t-r->megno_initial_t) * reb_tools_megno_deltad_delta(r);
}
for(size_t k=0;k<N;++k) {
size_t mk = map ? map[k] : k;
at[3*k] = particles[mk].ax;
at[3*k+1] = particles[mk].ay;
at[3*k+2] = particles[mk].az;
}
for(size_t mk=0;mk<N_var;++mk) {
size_t k = mk+N;;
at[3*k] = particles_var[mk].ax;
at[3*k+1] = particles_var[mk].ay;
at[3*k+2] = particles_var[mk].az;
}
switch (n) { // Improve b and g values
case 1:
for(size_t k=0;k<N3;++k) {
double tmp = g.p0[k];
double gk = at[k];
double gk_cs = ((double*)(gravity_cs))[k];
add_cs(&gk, &gk_cs, -a0[k]);
add_cs(&gk, &gk_cs, csa0[k]);
g.p0[k] = gk/rr[0];
add_cs(&(b.p0[k]), &(csb.p0[k]), g.p0[k]-tmp);
} break;
case 2:
for(size_t k=0;k<N3;++k) {
double tmp = g.p1[k];
double gk = at[k];
double gk_cs = ((double*)(gravity_cs))[k];
add_cs(&gk, &gk_cs, -a0[k]);
add_cs(&gk, &gk_cs, csa0[k]);
g.p1[k] = (gk/rr[1] - g.p0[k])/rr[2];
tmp = g.p1[k] - tmp;
add_cs(&(b.p0[k]), &(csb.p0[k]), tmp * c[0]);
add_cs(&(b.p1[k]), &(csb.p1[k]), tmp);
} break;
case 3:
for(size_t k=0;k<N3;++k) {
double tmp = g.p2[k];
double gk = at[k];
double gk_cs = ((double*)(gravity_cs))[k];
add_cs(&gk, &gk_cs, -a0[k]);
add_cs(&gk, &gk_cs, csa0[k]);
g.p2[k] = ((gk/rr[3] - g.p0[k])/rr[4] - g.p1[k])/rr[5];
tmp = g.p2[k] - tmp;
add_cs(&(b.p0[k]), &(csb.p0[k]), tmp * c[1]);
add_cs(&(b.p1[k]), &(csb.p1[k]), tmp * c[2]);
add_cs(&(b.p2[k]), &(csb.p2[k]), tmp);
} break;
case 4:
for(size_t k=0;k<N3;++k) {
double tmp = g.p3[k];
double gk = at[k];
double gk_cs = ((double*)(gravity_cs))[k];
add_cs(&gk, &gk_cs, -a0[k]);
add_cs(&gk, &gk_cs, csa0[k]);
g.p3[k] = (((gk/rr[6] - g.p0[k])/rr[7] - g.p1[k])/rr[8] - g.p2[k])/rr[9];
tmp = g.p3[k] - tmp;
add_cs(&(b.p0[k]), &(csb.p0[k]), tmp * c[3]);
add_cs(&(b.p1[k]), &(csb.p1[k]), tmp * c[4]);
add_cs(&(b.p2[k]), &(csb.p2[k]), tmp * c[5]);
add_cs(&(b.p3[k]), &(csb.p3[k]), tmp);
} break;
case 5:
for(size_t k=0;k<N3;++k) {
double tmp = g.p4[k];
double gk = at[k];
double gk_cs = ((double*)(gravity_cs))[k];
add_cs(&gk, &gk_cs, -a0[k]);
add_cs(&gk, &gk_cs, csa0[k]);
g.p4[k] = ((((gk/rr[10] - g.p0[k])/rr[11] - g.p1[k])/rr[12] - g.p2[k])/rr[13] - g.p3[k])/rr[14];
tmp = g.p4[k] - tmp;
add_cs(&(b.p0[k]), &(csb.p0[k]), tmp * c[6]);
add_cs(&(b.p1[k]), &(csb.p1[k]), tmp * c[7]);
add_cs(&(b.p2[k]), &(csb.p2[k]), tmp * c[8]);
add_cs(&(b.p3[k]), &(csb.p3[k]), tmp * c[9]);
add_cs(&(b.p4[k]), &(csb.p4[k]), tmp);
} break;
case 6:
for(size_t k=0;k<N3;++k) {
double tmp = g.p5[k];
double gk = at[k];
double gk_cs = ((double*)(gravity_cs))[k];
add_cs(&gk, &gk_cs, -a0[k]);
add_cs(&gk, &gk_cs, csa0[k]);
g.p5[k] = (((((gk/rr[15] - g.p0[k])/rr[16] - g.p1[k])/rr[17] - g.p2[k])/rr[18] - g.p3[k])/rr[19] - g.p4[k])/rr[20];
tmp = g.p5[k] - tmp;
add_cs(&(b.p0[k]), &(csb.p0[k]), tmp * c[10]);
add_cs(&(b.p1[k]), &(csb.p1[k]), tmp * c[11]);
add_cs(&(b.p2[k]), &(csb.p2[k]), tmp * c[12]);
add_cs(&(b.p3[k]), &(csb.p3[k]), tmp * c[13]);
add_cs(&(b.p4[k]), &(csb.p4[k]), tmp * c[14]);
add_cs(&(b.p5[k]), &(csb.p5[k]), tmp);
} break;
case 7:
{
double maxak = 0.0;
double maxb6ktmp = 0.0;
for(size_t k=0;k<N3;++k) {
double tmp = g.p6[k];
double gk = at[k];
double gk_cs = ((double*)(gravity_cs))[k];
add_cs(&gk, &gk_cs, -a0[k]);
add_cs(&gk, &gk_cs, csa0[k]);
g.p6[k] = ((((((gk/rr[21] - g.p0[k])/rr[22] - g.p1[k])/rr[23] - g.p2[k])/rr[24] - g.p3[k])/rr[25] - g.p4[k])/rr[26] - g.p5[k])/rr[27];
tmp = g.p6[k] - tmp;
add_cs(&(b.p0[k]), &(csb.p0[k]), tmp * c[15]);
add_cs(&(b.p1[k]), &(csb.p1[k]), tmp * c[16]);
add_cs(&(b.p2[k]), &(csb.p2[k]), tmp * c[17]);
add_cs(&(b.p3[k]), &(csb.p3[k]), tmp * c[18]);
add_cs(&(b.p4[k]), &(csb.p4[k]), tmp * c[19]);
add_cs(&(b.p5[k]), &(csb.p5[k]), tmp * c[20]);
add_cs(&(b.p6[k]), &(csb.p6[k]), tmp);
// Monitor change in b.p6[k] relative to at[k]. The predictor corrector scheme is converged if it is close to 0.
if (ias15->adaptive_mode!=REB_IAS15_ADAPTIVEMODE_INDIVIDUAL){
const double ak = fabs(at[k]);
if (isnormal(ak) && ak>maxak){
maxak = ak;
}
const double b6ktmp = fabs(tmp); // change of b6ktmp coefficient
if (isnormal(b6ktmp) && b6ktmp>maxb6ktmp){
maxb6ktmp = b6ktmp;
}
}else{
const double ak = at[k];
const double b6ktmp = tmp;
const double errork = fabs(b6ktmp/ak);
if (isnormal(errork) && errork>predictor_corrector_error){
predictor_corrector_error = errork;
}
}
}
if (ias15->adaptive_mode!=REB_IAS15_ADAPTIVEMODE_INDIVIDUAL){
predictor_corrector_error = maxb6ktmp/maxak;
}
break;
}
}
}
}
// Set time back to initial value (will be updated below)
r->t = t_beginning;
// Find new timestep
const double dt_done = r->dt;
double dt_new;
if (ias15->epsilon>0){
// Estimate error (given by last term in series expansion)
// There are two options:
// ias15->adaptive_mode==REB_IAS15_ADAPTIVEMODE_GLOBAL (used to be default until January 2024)
// First, we determine the maximum acceleration and the maximum of the last term in the series.
// Then, the two are divided.
// ias15->adaptive_mode==REB_IAS15_ADAPTIVEMODE_INDIVIDUAL
// Here, the fractional error is calculated for each particle individually and we use the maximum of the fractional error.
// This might fail in cases where a particle does not experience any (physical) acceleration besides roundoff errors.
// ias15->adaptive_mode==REB_IAS15_ADAPTIVEMODE_PRS23
// Here, the acceleration, jerk and snap are used to estimate the new timestep.
// The method is described in detail in Pham, Rein, Spiegel 2023
if (ias15->adaptive_mode==REB_IAS15_ADAPTIVEMODE_INDIVIDUAL || ias15->adaptive_mode==REB_IAS15_ADAPTIVEMODE_GLOBAL){ // Old adaptive timestepping methods
double integrator_error = 0.0; // Try to estimate integrator error based on last polynomial
if (ias15->adaptive_mode==REB_IAS15_ADAPTIVEMODE_GLOBAL){
double maxa = 0.0;
double maxj = 0.0;
for(size_t i=0;i<N;i++){ // Looping over all particles and all 3 components of the acceleration.
size_t mi = map ? map[i] : i;
const double v2 = particles[mi].vx*particles[mi].vx+particles[mi].vy*particles[mi].vy+particles[mi].vz*particles[mi].vz;
const double x2 = particles[mi].x*particles[mi].x+particles[mi].y*particles[mi].y+particles[mi].z*particles[mi].z;
// Skip slowly varying accelerations
if (fabs(v2*r->dt*r->dt/x2) < 1e-16) continue;
for(size_t k=3*i;k<3*(i+1);k++) {
const double ak = fabs(at[k]);
if (isnormal(ak) && ak>maxa){
maxa = ak;
}
const double b6k = fabs(b.p6[k]);
if (isnormal(b6k) && b6k>maxj){
maxj = b6k;
}
}
integrator_error = maxj/maxa;
}
}else{ // adaptive_mode == REB_IAS15_ADAPTIVEMODE_INDIVIDUAL
for(size_t k=0;k<N3;k++) {
const double ak = at[k];
const double bk = b.p6[k];
const double errork = fabs(bk/ak);
if (isnormal(errork) && errork>integrator_error){
integrator_error = errork;
}
}
}
// Use error estimate to predict new timestep
if (isnormal(integrator_error)){
dt_new = sqrt7(ias15->epsilon/integrator_error)*dt_done;
}else{ // In the rare case that the error estimate doesn't give a finite number (e.g. when all forces accidentally cancel up to machine precision).
dt_new = dt_done/safety_factor; // by default, increase timestep a little
};
}else{ // adaptive_mode >= 2 (New adaptive timestepping method, default since January 2024)
double min_timescale2 = INFINITY; // note factor of dt_done**2 not included
for(size_t i=0;i<N;i++){
double a0i = 0; // (acceleration at beginning of timestep)^2
double y2 = 0; // (accaleration at end of timestep)^2
double y3 = 0; // (jerk * dt_done)^2
double y4 = 0; // (snap * dt_done^2)^2
double y5 = 0; // (crackle * dt_done^3)^2
for(size_t k=3*i;k<3*(i+1);k++) {
a0i += a0[k]*a0[k];
double tmp = a0[k] + b.p0[k] + b.p1[k] + b.p2[k] + b.p3[k] + b.p4[k] + b.p5[k] + b.p6[k];
y2 += tmp*tmp;
tmp = b.p0[k] + 2.* b.p1[k] + 3.* b.p2[k] + 4.* b.p3[k] + 5.* b.p4[k] + 6.* b.p5[k] + 7.* b.p6[k];
y3 += tmp*tmp;
tmp = 2.* b.p1[k] + 6.* b.p2[k] + 12.* b.p3[k] + 20.* b.p4[k] + 30.* b.p5[k] + 42.* b.p6[k];
y4 += tmp*tmp;
tmp = 6.* b.p2[k] + 24.* b.p3[k] + 60.* b.p4[k] + 120.* b.p5[k] + 210.* b.p6[k];
y5 += tmp*tmp;
}
if (!isnormal(a0i)){
// Skip particles which do not experience any acceleration or
// have acceleration which is inf or Nan.
continue;
}
double timescale2 = 0;
if (ias15->adaptive_mode == REB_IAS15_ADAPTIVEMODE_PRS23){
timescale2 = 2.*y2/(y3+sqrt(y4*y2)); // PRS23
}else if (ias15->adaptive_mode == REB_IAS15_ADAPTIVEMODE_AARSETH85){
timescale2 = (sqrt(y2*y4)+y3) / (sqrt(y3*y5)+y4); // A85
}
if (isnormal(timescale2) && timescale2<min_timescale2){
min_timescale2 = timescale2;
}
}
if (isnormal(min_timescale2)){
// Numerical factor below is there to match timestep to that of adaptive_mode==REB_IAS15_ADAPTIVEMODE_GLOBAL and default epsilon
dt_new = sqrt(min_timescale2) * dt_done * sqrt7(ias15->epsilon*5040.0);
}else{
dt_new = dt_done/safety_factor; // by default, increase timestep a little
}
}
if (fabs(dt_new)<ias15->min_dt) dt_new = copysign(ias15->min_dt,dt_new);
if (fabs(dt_new/dt_done) < safety_factor) { // New timestep is significantly smaller.
// Reset particles
for(size_t k=0;k<N;++k) {
size_t mk = map ? map[k] : k;
particles[mk].x = x0[3*k+0]; // Set initial position
particles[mk].y = x0[3*k+1];
particles[mk].z = x0[3*k+2];
particles[mk].vx = v0[3*k+0]; // Set initial velocity
particles[mk].vy = v0[3*k+1];
particles[mk].vz = v0[3*k+2];
particles[mk].ax = a0[3*k+0]; // Set initial acceleration
particles[mk].ay = a0[3*k+1];
particles[mk].az = a0[3*k+2];
}
for(size_t mk=0;mk<N_var;++mk) {
size_t k = mk + N;
particles_var[mk].x = x0[3*k+0]; // Set initial position
particles_var[mk].y = x0[3*k+1];
particles_var[mk].z = x0[3*k+2];
particles_var[mk].vx = v0[3*k+0]; // Set initial velocity
particles_var[mk].vy = v0[3*k+1];
particles_var[mk].vz = v0[3*k+2];
particles_var[mk].ax = a0[3*k+0]; // Set initial acceleration
particles_var[mk].ay = a0[3*k+1];
particles_var[mk].az = a0[3*k+2];
}
r->dt = dt_new;
if (r->dt_last_done!=0.){ // Do not predict next e/b values if this is the first time step.
double ratio = r->dt/r->dt_last_done;
predict_next_step(ratio, N3, er, br, e, b);
}
return 0; // Step rejected. Do again.
}
if (fabs(dt_new/dt_done) > 1.0) { // New timestep is larger.
if (dt_new/dt_done > 1./safety_factor){
dt_new = dt_done /safety_factor; // Don't increase the timestep by too much compared to the last one.
}
}
r->dt = dt_new;
}
// Find new position and velocity values at end of the sequence
for(size_t k=0;k<N3;++k) {
// Note: dt_done*dt_done is not precalculated to avoid
// biased round-off errors when a fixed timestep is used.
add_cs(&(x0[k]), &(csx[k]), b.p6[k]/72.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), b.p5[k]/56.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), b.p4[k]/42.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), b.p3[k]/30.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), b.p2[k]/20.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), b.p1[k]/12.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), b.p0[k]/6.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), a0[k]/2.*dt_done*dt_done);
add_cs(&(x0[k]), &(csx[k]), v0[k]*dt_done);
add_cs(&(v0[k]), &(csv[k]), b.p6[k]/8.*dt_done);
add_cs(&(v0[k]), &(csv[k]), b.p5[k]/7.*dt_done);
add_cs(&(v0[k]), &(csv[k]), b.p4[k]/6.*dt_done);
add_cs(&(v0[k]), &(csv[k]), b.p3[k]/5.*dt_done);
add_cs(&(v0[k]), &(csv[k]), b.p2[k]/4.*dt_done);
add_cs(&(v0[k]), &(csv[k]), b.p1[k]/3.*dt_done);
add_cs(&(v0[k]), &(csv[k]), b.p0[k]/2.*dt_done);
add_cs(&(v0[k]), &(csv[k]), a0[k]*dt_done);
}
r->t += dt_done;
r->dt_last_done = dt_done;
if (r->calculate_megno){
double dY = dt_done*integrator_megno_thisdt;
reb_tools_megno_update(r, dY, dt_done);
}
// Swap particle buffers
for(size_t k=0;k<N;++k) {
size_t mk = map ? map[k] : k;
particles[mk].x = x0[3*k+0]; // Set final position
particles[mk].y = x0[3*k+1];
particles[mk].z = x0[3*k+2];
particles[mk].vx = v0[3*k+0]; // Set final velocity
particles[mk].vy = v0[3*k+1];
particles[mk].vz = v0[3*k+2];
}
for(size_t mk=0;mk<N_var;++mk) {
size_t k = mk + N;
particles_var[mk].x = x0[3*k+0]; // Set final position
particles_var[mk].y = x0[3*k+1];
particles_var[mk].z = x0[3*k+2];
particles_var[mk].vx = v0[3*k+0]; // Set final velocity
particles_var[mk].vy = v0[3*k+1];
particles_var[mk].vz = v0[3*k+2];
}
copybuffers(e,er,N3);
copybuffers(b,br,N3);
double ratio = r->dt/dt_done;
predict_next_step(ratio, N3, e, b, e, b);
return 1; // Success.
}
static void predict_next_step(double ratio, size_t N3, const struct reb_dpconst7 _e, const struct reb_dpconst7 _b, const struct reb_dpconst7 e, const struct reb_dpconst7 b){
if (ratio>20.){
// Do not predict if stepsize increase is very large.
for(size_t k=0;k<N3;++k) {
e.p0[k] = 0.; e.p1[k] = 0.; e.p2[k] = 0.; e.p3[k] = 0.; e.p4[k] = 0.; e.p5[k] = 0.; e.p6[k] = 0.;
b.p0[k] = 0.; b.p1[k] = 0.; b.p2[k] = 0.; b.p3[k] = 0.; b.p4[k] = 0.; b.p5[k] = 0.; b.p6[k] = 0.;
}
}else{
// Predict new B values to use at the start of the next sequence. The predicted
// values from the last call are saved as E. The correction, BD, between the
// actual and predicted values of B is applied in advance as a correction.
//
const double q1 = ratio;
const double q2 = q1 * q1;
const double q3 = q1 * q2;
const double q4 = q2 * q2;
const double q5 = q2 * q3;
const double q6 = q3 * q3;
const double q7 = q3 * q4;
for(size_t k=0;k<N3;++k) {
double be0 = _b.p0[k] - _e.p0[k];
double be1 = _b.p1[k] - _e.p1[k];
double be2 = _b.p2[k] - _e.p2[k];
double be3 = _b.p3[k] - _e.p3[k];
double be4 = _b.p4[k] - _e.p4[k];
double be5 = _b.p5[k] - _e.p5[k];
double be6 = _b.p6[k] - _e.p6[k];
e.p0[k] = q1*(_b.p6[k]* 7.0 + _b.p5[k]* 6.0 + _b.p4[k]* 5.0 + _b.p3[k]* 4.0 + _b.p2[k]* 3.0 + _b.p1[k]*2.0 + _b.p0[k]);
e.p1[k] = q2*(_b.p6[k]*21.0 + _b.p5[k]*15.0 + _b.p4[k]*10.0 + _b.p3[k]* 6.0 + _b.p2[k]* 3.0 + _b.p1[k]);
e.p2[k] = q3*(_b.p6[k]*35.0 + _b.p5[k]*20.0 + _b.p4[k]*10.0 + _b.p3[k]* 4.0 + _b.p2[k]);
e.p3[k] = q4*(_b.p6[k]*35.0 + _b.p5[k]*15.0 + _b.p4[k]* 5.0 + _b.p3[k]);
e.p4[k] = q5*(_b.p6[k]*21.0 + _b.p5[k]* 6.0 + _b.p4[k]);
e.p5[k] = q6*(_b.p6[k]* 7.0 + _b.p5[k]);
e.p6[k] = q7* _b.p6[k];
b.p0[k] = e.p0[k] + be0;
b.p1[k] = e.p1[k] + be1;
b.p2[k] = e.p2[k] + be2;
b.p3[k] = e.p3[k] + be3;
b.p4[k] = e.p4[k] + be4;
b.p5[k] = e.p5[k] + be5;
b.p6[k] = e.p6[k] + be6;
}
}
}
static void copybuffers(const struct reb_dpconst7 _a, const struct reb_dpconst7 _b, size_t N3){
for (size_t i=0;i<N3*7;i++){
_b.p0[i] = _a.p0[i];
}
// The above code seems faster than the code below, probably due to some compiler optimizations.
// for (size_t i=0;i<7;i++){
// memcpy(_b[i],_a[i], sizeof(double)*N3);
// }
}
// Do nothing here. This is only used in a leapfrog-like DKD integrator. IAS15 performs one complete timestep.
void reb_integrator_ias15_step(struct reb_simulation* r, void* state){
r->gravity_ignore_terms = REB_GRAVITY_IGNORE_TERMS_NONE;
if (r->N){
// Try until a step was successful.
while(!reb_integrator_ias15_step_try(r, state));
}else{
r->t += r->dt;
r->dt_last_done = r->dt;
}
}
double reb_integrator_ias15_timescale(struct reb_simulation* r){
// Returns a timescale according to Pham, Rein, Spiegel 2023 (PRS23)
reb_simulation_update_acceleration(r);
size_t N; // Number of real particles
size_t* map = r->map; // this map allow for integrating only a selection of particles
if (map){
N = r->N_map;
}else{
N = r->N;
}
double min_timescale2 = INFINITY;
for (size_t i=0; i<N; i++) {
size_t mi = map ? map[i] : i;
struct reb_particle* p_i = &(r->particles[mi]);
double y2 = p_i->ax*p_i->ax + p_i->ay*p_i->ay + p_i->az*p_i->az;
struct reb_vec3d vec_y3 = {0};
struct reb_vec3d vec_y4 = {0};
if (!isnormal(y2)){
// Skip particles which do not experience any acceleration or
// have acceleration which is inf or Nan.
continue;
}
for (size_t j=0; j<N; j++) { // N^2 loop could be optimized to N^2/2
size_t mj = map ? map[j] : j;
if (mi==mj) continue;
struct reb_particle* p_j = &(r->particles[mj]);
double rij_x = p_j->x - p_i->x;
double rij_y = p_j->y - p_i->y;
double rij_z = p_j->z - p_i->z;
double vij_x = p_j->vx - p_i->vx;
double vij_y = p_j->vy - p_i->vy;
double vij_z = p_j->vz - p_i->vz;
double aij_x = p_j->ax - p_i->ax;
double aij_y = p_j->ay - p_i->ay;
double aij_z = p_j->az - p_i->az;
double r_sq = rij_x * rij_x + rij_y * rij_y + rij_z * rij_z;
double r_mag = sqrt(r_sq); // |r_ij|
double r_cubed = r_sq * r_mag; // |r_ij|^3
double r_fifth = r_cubed * r_sq; // |r_ij|^5
double r_seventh = r_fifth * r_sq;// |r_ij|^7
// Dot products
double r_dot_v = rij_x * vij_x + rij_y * vij_y + rij_z * vij_z; // (r_ij . v_ij)
double r_dot_a = rij_x * aij_x + rij_y * aij_y + rij_z * aij_z; // (r_ij . a_ij)
double v_sq = vij_x * vij_x + vij_y * vij_y + vij_z * vij_z; // v_ij^2
// --- Jerk Calculation for particle i due to particle j ---
// Term 1: v_ij / |r_ij|^3
// Term 2: -3 * r_ij * (r_ij . v_ij) / |r_ij|^5
double jerk_factor1 = p_j->m / r_cubed;
double jerk_factor2 = -3.0 * p_j->m * r_dot_v / r_fifth;
vec_y3.x += jerk_factor1 * vij_x + jerk_factor2 * rij_x;
vec_y3.y += jerk_factor1 * vij_y + jerk_factor2 * rij_y;
vec_y3.z += jerk_factor1 * vij_z + jerk_factor2 * rij_z;
// --- Snap Calculation for particle i due to particle j ---
// Term 1: a_ij / |r_ij|^3
// Term 2: -6 * v_ij * (r_ij . v_ij) / |r_ij|^5
// Term 3: -3 * r_ij * v_ij^2 / |r_ij|^5
// Term 4: -3 * r_ij * (r_ij . a_ij) / |r_ij|^5
// Term 5: +15 * r_ij * (r_ij . v_ij)^2 / |r_ij|^7
double snap_c1 = p_j->m / r_cubed; // for a_ij term
double snap_c2 = -6.0 * p_j->m * r_dot_v / r_fifth; // for v_ij term
double snap_c3_rij = -3.0 * p_j->m * v_sq / r_fifth; // for r_ij term (from v_ij^2)
double snap_c4_rij = -3.0 * p_j->m * r_dot_a / r_fifth; // for r_ij term (from r_ij . a_ij)
double snap_c5_rij = 15.0 * p_j->m * r_dot_v * r_dot_v / r_seventh; // for r_ij term (from (r_ij . v_ij)^2)
vec_y4.x += snap_c1 * aij_x + snap_c2 * vij_x + (snap_c3_rij + snap_c4_rij + snap_c5_rij) * rij_x;
vec_y4.y += snap_c1 * aij_y + snap_c2 * vij_y + (snap_c3_rij + snap_c4_rij + snap_c5_rij) * rij_y;
vec_y4.z += snap_c1 * aij_z + snap_c2 * vij_z + (snap_c3_rij + snap_c4_rij + snap_c5_rij) * rij_z;
}
vec_y3.x *= r->G;
vec_y3.y *= r->G;
vec_y3.z *= r->G;
vec_y4.x *= r->G;
vec_y4.y *= r->G;
vec_y4.z *= r->G;
double y3 = vec_y3.x*vec_y3.x + vec_y3.y*vec_y3.y + vec_y3.z*vec_y3.z;
double y4 = vec_y4.x*vec_y4.x + vec_y4.y*vec_y4.y + vec_y4.z*vec_y4.z;
double timescale2 = 2.*y2/(y3+sqrt(y4*y2)); // PRS23
if (isnormal(timescale2) && timescale2<min_timescale2){
min_timescale2 = timescale2;
}
}
return sqrt(min_timescale2); // Multiply by sqrt7(ias15->epsilon*5040.0) to get IAS default timestep.
}
#ifdef GENERATE_CONSTANTS
#include <gmp.h>
void integrator_generate_constants(void){
printf("Generating constants.\n\n");
mpf_set_default_prec(512);
mpf_t* _h = malloc(sizeof(mpf_t)*8);
for (size_t i=0;i<8;i++){
mpf_init(_h[i]);
}
mpf_t* _r = malloc(sizeof(mpf_t)*28);
for (size_t i=0;i<28;i++){
mpf_init(_r[i]);
}
mpf_t* _c = malloc(sizeof(mpf_t)*21);
mpf_t* _d = malloc(sizeof(mpf_t)*21);
for (size_t i=0;i<21;i++){
mpf_init(_c[i]);
mpf_init(_d[i]);
}
mpf_set_str(_h[0],"0.0",10);
mpf_set_str(_h[1],"0.0562625605369221464656521910318",10);
mpf_set_str(_h[2],"0.180240691736892364987579942780",10);
mpf_set_str(_h[3],"0.352624717113169637373907769648",10);
mpf_set_str(_h[4],"0.547153626330555383001448554766",10);
mpf_set_str(_h[5],"0.734210177215410531523210605558", 10);
mpf_set_str(_h[6],"0.885320946839095768090359771030",10);
mpf_set_str(_h[7],"0.977520613561287501891174488626",10);
size_t l=0;
for (size_t j=1;j<8;++j) {
for(size_t k=0;k<j;++k) {
// rr[l] = h[j] - h[k];
mpf_sub(_r[l],_h[j],_h[k]);
++l;
}
}
//c[0] = -h[1];
mpf_neg(_c[0],_h[1]);
//d[0] = h[1];
mpf_set(_d[0],_h[1]);
l=0;
for (size_t j=2;j<7;++j) {
++l;
// c[l] = -h[j] * c[l-j+1];
mpf_mul(_c[l], _h[j], _c[l-j+1]);
mpf_neg(_c[l], _c[l]);
//d[l] = h[1] * d[l-j+1];
mpf_mul(_d[l], _h[1], _d[l-j+1]);
for(size_t k=2;k<j;++k) {
++l;
//c[l] = c[l-j] - h[j] * c[l-j+1];
mpf_mul(_c[l], _h[j], _c[l-j+1]);
mpf_sub(_c[l], _c[l-j], _c[l]);
//d[l] = d[l-j] + h[k] * d[l-j+1];
mpf_mul(_d[l], _h[k], _d[l-j+1]);
mpf_add(_d[l], _d[l-j], _d[l]);
}
++l;
//c[l] = c[l-j] - h[j];
mpf_sub(_c[l], _c[l-j], _h[j]);
//d[l] = d[l-j] + h[j];
mpf_add(_d[l], _d[l-j], _h[j]);
}
// Output
printf("double rr[28] = {");
for (size_t i=0;i<28;i++){
gmp_printf ("%.*Ff", 25, _r[i]);
if (i!=27) printf(", ");
}
printf("};\n");
printf("double c[21] = {");
for (size_t i=0;i<21;i++){
gmp_printf ("%.*Ff", 25, _c[i]);
if (i!=20) printf(", ");
}
printf("};\n");
printf("double d[21] = {");
for (size_t i=0;i<21;i++){
gmp_printf ("%.*Ff", 25, _d[i]);
if (i!=20) printf(", ");
}
printf("};\n");
exit(0);
}
#endif // GENERATE_CONSTANTS