rebound-bind 5.0.0

/**
 * @file    integrator_eos.c
 * @brief   Embedded Operator Splitting (EOS) method
 * @author  Hanno Rein
 * 
 * @section LICENSE
 * Copyright (c) 2019 Hanno Rein
 *
 * 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>
#include "gravity.h"
#include "integrator_eos.h"
#include "integrator_leapfrog.h"
#include "tools.h"
#include "binarydata.h"
#define MIN(a, b) ((a) > (b) ? (b) : (a))    ///< Returns the minimum of a and b
#define MAX(a, b) ((a) > (b) ? (a) : (b))    ///< Returns the maximum of a and b

void reb_integrator_eos_step(struct reb_simulation* r, void* state);
void reb_integrator_eos_synchronize(struct reb_simulation* r, void* state);
void* reb_integrator_eos_create();
const struct reb_binarydata_field_descriptor reb_integrator_eos_field_descriptor_list[];

const struct reb_integrator reb_integrator_eos = {
    .documentation = 
    "This is the Embedded Operator Splitting (EOS) methods described in [Rein (2019)]."
    "\n\n"
    "[Rein (2019)]: https://ui.adsabs.harvard.edu/abs/2020MNRAS.492.5413R/abstract\n"
    ,
    .step = reb_integrator_eos_step,
    .create = reb_integrator_eos_create,
    .synchronize = reb_integrator_eos_synchronize,
    .field_descriptor_list = reb_integrator_eos_field_descriptor_list,
};

const struct reb_binarydata_field_descriptor reb_integrator_eos_field_descriptor_list[] = {
    { "Outer operator splitting scheme",
        REB_INT,         "phi0",                  offsetof(struct reb_integrator_eos_state, phi0), 0, 0, REB_GENERATE_ENUM_DESCRIPTORS(REB_INTEGRATOR_EOS_TYPE)},
    { "Inner operator splitting scheme", 
        REB_INT,         "phi1",                  offsetof(struct reb_integrator_eos_state, phi1), 0, 0, REB_GENERATE_ENUM_DESCRIPTORS(REB_INTEGRATOR_EOS_TYPE)},
    { "Number of sub-timesteps. Default: 2.",
        REB_UINT,        "n",                     offsetof(struct reb_integrator_eos_state, n), 0, 0, 0},
    { "If set to 0, always combine drift steps at the beginning and end of `phi0`. If set to 1, `n` needs to be bigger than 1.", 
        REB_UINT,        "safe_mode",             offsetof(struct reb_integrator_eos_state, safe_mode), 0, 0, 0},
    { 0 }, // Null terminated list
};

void* reb_integrator_eos_create(){
    struct reb_integrator_eos_state* eos = calloc(sizeof(struct reb_integrator_eos_state),1);
    eos->n = 2;
    eos->phi0 = REB_INTEGRATOR_EOS_TYPE_LF;
    eos->phi1 = REB_INTEGRATOR_EOS_TYPE_LF;
    eos->safe_mode = 1;
    return eos;
}


static const double lf4_2_a = 0.211324865405187117745425609749;

static const double lf8_6_4_a[4] = {0.0711334264982231177779387300061549964174, 0.241153427956640098736487795326289649618, 0.521411761772814789212136078067994229991, -0.333698616227678005726562603400438876027};
static const double lf8_6_4_b[4] = {0.183083687472197221961703757166430291072, 0.310782859898574869507522291054262796375, -0.0265646185119588006972121379164987592663, 0.0653961422823734184559721793911134363710};

static const double pmlf6_a[2] = {-0.0682610383918630,0.568261038391863038121699}; 
static const double pmlf6_b[2] = {0.2621129352517028, 0.475774129496594366806050}; 
static const double pmlf6_c[2] = {0., 0.0164011128160783}; 
static const double pmlf6_z[6] = { 0.07943288242455420, 0.02974829169467665, -0.7057074964815896, 0.3190423451260838, -0.2869147334299646, 0.564398710666239478150885};
static const double pmlf6_y[6] = {1.3599424487455264, -0.6505973747535132, -0.033542814598338416, -0.040129915275115030, 0.044579729809902803, -0.680252073928462652752103};
static const double pmlf6_v[6] = {-0.034841228074994859, 0.031675672097525204, -0.005661054677711889, 0.004262222269023640, 0.005, -0.005};

static const double pmlf4_y[3] = {0.1859353996846055, 0.0731969797858114, -0.1576624269298081};
static const double pmlf4_z[3] = {0.8749306155955435, -0.237106680151022, -0.5363539829039128};

static const double plf7_6_4_a[2] = {0.5600879810924619,-0.060087981092461900000};
static const double plf7_6_4_b[2] = {1.5171479707207228, -2.0342959414414456000};
static const double plf7_6_4_z[6] = {-0.3346222298730800, 1.0975679907321640, -1.0380887460967830, 0.6234776317921379, -1.1027532063031910, -0.0141183222088869};
static const double plf7_6_4_y[6] = {-1.6218101180868010, 0.0061709468110142, 0.8348493592472594, -0.0511253369989315, 0.5633782670698199, -0.5};

static inline void reb_integrator_eos_interaction_shell0(struct reb_simulation* r, double y, double v){
    // Calculate gravity using standard gravity routine
    r->gravity_ignore_terms = REB_GRAVITY_IGNORE_TERMS_INVOLVING_0;
    r->gravity = REB_GRAVITY_BASIC;
    reb_simulation_update_acceleration(r);
    if (v!=0.){
        reb_gravity_basic_calculate_and_apply_jerk(r,v);
    }
    // Apply acceleration (jerk already applied)
    struct reb_particle* restrict const particles = r->particles;
    const size_t N = r->N;
    for (size_t i=0;i<N;i++){
        particles[i].vx += y*particles[i].ax;
        particles[i].vy += y*particles[i].ay;
        particles[i].vz += y*particles[i].az;
    }
}

static inline void reb_integrator_eos_interaction_shell1(struct reb_simulation* r, double y, double v){
    const size_t N = r->N;
    const size_t N_active = r->N_active==SIZE_MAX?N:r->N_active;
    const int testparticle_type   = r->testparticle_type;
    struct reb_particle* restrict const particles = r->particles;

    const double G = r->G;

    if (v!=0.){ // is jerk even used?
                // Normal force calculation 
        particles[0].ax = 0;
        particles[0].ay = 0;
        particles[0].az = 0;
        // Interactions between central object and all other active particles
        for (size_t j=1; j<N_active; j++){
            const double dx = particles[0].x - particles[j].x;
            const double dy = particles[0].y - particles[j].y;
            const double dz = particles[0].z - particles[j].z;
            const double dr = sqrt(dx*dx + dy*dy + dz*dz);

            const double prefact = G/(dr*dr*dr);
            const double prefactj = -prefact*particles[j].m;
            particles[0].ax    += prefactj*dx;
            particles[0].ay    += prefactj*dy;
            particles[0].az    += prefactj*dz;
            const double prefacti = prefact*particles[0].m;
            particles[j].ax    = prefacti*dx;
            particles[j].ay    = prefacti*dy;
            particles[j].az    = prefacti*dz;
        }
        // Interactions between central object and all test particles
        for (size_t j=N_active; j<N; j++){
            const double dx = particles[0].x - particles[j].x;
            const double dy = particles[0].y - particles[j].y;
            const double dz = particles[0].z - particles[j].z;
            const double dr = sqrt(dx*dx + dy*dy + dz*dz);

            const double prefact = G/(dr*dr*dr);
            const double prefacti = prefact*particles[0].m;
            particles[j].ax    = prefacti*dx;
            particles[j].ay    = prefacti*dy;
            particles[j].az    = prefacti*dz;
            if (testparticle_type){
                const double prefactj = -prefact*particles[j].m;
                particles[0].ax    += prefactj*dx;
                particles[0].ay    += prefactj*dy;
                particles[0].az    += prefactj*dz;
            }
        }
        // Jerk calculation
        // Interactions between central object and all other active particles
        for (size_t i=1; i<N_active; i++){
            const double dx = particles[0].x - particles[i].x; 
            const double dy = particles[0].y - particles[i].y; 
            const double dz = particles[0].z - particles[i].z; 

            const double dax = particles[0].ax - particles[i].ax; 
            const double day = particles[0].ay - particles[i].ay; 
            const double daz = particles[0].az - particles[i].az; 

            const double dr = sqrt(dx*dx + dy*dy + dz*dz);
            const double alphasum = dax*dx+day*dy+daz*dz;
            const double prefact2 = 2.*v*G /(dr*dr*dr);
            const double prefact2i = prefact2*particles[i].m;
            const double prefact2j = prefact2*particles[0].m;
            const double prefact1 = alphasum*prefact2/dr *3./dr;
            const double prefact1i = prefact1*particles[i].m;
            const double prefact1j = prefact1*particles[0].m;
            particles[0].vx    += -dax*prefact2i + dx*prefact1i;
            particles[0].vy    += -day*prefact2i + dy*prefact1i;
            particles[0].vz    += -daz*prefact2i + dz*prefact1i;
            particles[i].vx    += y*particles[i].ax + dax*prefact2j - dx*prefact1j;
            particles[i].vy    += y*particles[i].ay + day*prefact2j - dy*prefact1j;
            particles[i].vz    += y*particles[i].az + daz*prefact2j - dz*prefact1j;
        }
        // Interactions between central object and all test particles
        for (size_t i=N_active; i<N; i++){
            const double dx = particles[0].x - particles[i].x; 
            const double dy = particles[0].y - particles[i].y; 
            const double dz = particles[0].z - particles[i].z; 

            const double dax = particles[0].ax - particles[i].ax; 
            const double day = particles[0].ay - particles[i].ay; 
            const double daz = particles[0].az - particles[i].az; 

            const double dr = sqrt(dx*dx + dy*dy + dz*dz);
            const double alphasum = dax*dx+day*dy+daz*dz;
            const double prefact2 = 2.*v*G /(dr*dr*dr);
            const double prefact2j = prefact2*particles[0].m;
            const double prefact1 = alphasum*prefact2/dr *3./dr;
            const double prefact1j = prefact1*particles[0].m;
            if (testparticle_type){
                const double prefact2i = prefact2*particles[i].m;
                const double prefact1i = prefact1*particles[i].m;
                particles[0].vx    += -dax*prefact2i + dx*prefact1i;
                particles[0].vy    += -day*prefact2i + dy*prefact1i;
                particles[0].vz    += -daz*prefact2i + dz*prefact1i;
            }
            particles[i].vx    += y*particles[i].ax + dax*prefact2j - dx*prefact1j;
            particles[i].vy    += y*particles[i].ay + day*prefact2j - dy*prefact1j;
            particles[i].vz    += y*particles[i].az + daz*prefact2j - dz*prefact1j;
        }
        particles[0].vx += y*particles[0].ax;
        particles[0].vy += y*particles[0].ay;
        particles[0].vz += y*particles[0].az;
    }else{
        // Normal force calculation 
        // Interactions between central object and all other active particles
        for (size_t j=1; j<N_active; j++){
            const double dx = particles[0].x - particles[j].x;
            const double dy = particles[0].y - particles[j].y;
            const double dz = particles[0].z - particles[j].z;
            const double dr = sqrt(dx*dx + dy*dy + dz*dz);

            const double prefact = y*G/(dr*dr*dr);
            const double prefactj = -prefact*particles[j].m;
            particles[0].vx    += prefactj*dx;
            particles[0].vy    += prefactj*dy;
            particles[0].vz    += prefactj*dz;
            const double prefacti = prefact*particles[0].m;
            particles[j].vx    += prefacti*dx;
            particles[j].vy    += prefacti*dy;
            particles[j].vz    += prefacti*dz;
        }
        // Interactions between central object and all test particles
        for (size_t j=N_active; j<N; j++){
            const double dx = particles[0].x - particles[j].x;
            const double dy = particles[0].y - particles[j].y;
            const double dz = particles[0].z - particles[j].z;
            const double dr = sqrt(dx*dx + dy*dy + dz*dz);

            const double prefact = y*G/(dr*dr*dr);
            const double prefacti = prefact*particles[0].m;
            particles[j].vx    += prefacti*dx;
            particles[j].vy    += prefacti*dy;
            particles[j].vz    += prefacti*dz;
            if (testparticle_type){
                const double prefactj = -prefact*particles[j].m;
                particles[0].vx    += prefactj*dx;
                particles[0].vy    += prefactj*dy;
                particles[0].vz    += prefactj*dz;
            }
        }
    }
}

static inline void reb_integrator_eos_preprocessor(struct reb_simulation* const r, double dt, enum REB_INTEGRATOR_EOS_TYPE type, void (*drift_step)(struct reb_simulation* const r, double a), void (*interaction_step)(struct reb_simulation* const r, double y, double v)){
    switch(type){
        case REB_INTEGRATOR_EOS_TYPE_PMLF6:
            for (int i=0;i<6;i++){
                drift_step(r, dt*pmlf6_z[i]);
                interaction_step(r, dt*pmlf6_y[i], dt*dt*dt*pmlf6_v[i]);
            }
            break;
        case REB_INTEGRATOR_EOS_TYPE_PMLF4:
            for (int i=0;i<3;i++){
                interaction_step(r, dt*pmlf4_y[i], 0.);
                drift_step(r, dt*pmlf4_z[i]);
            }
            break;
        case REB_INTEGRATOR_EOS_TYPE_PLF7_6_4:
            for (int i=0;i<6;i++){
                drift_step(r, dt*plf7_6_4_z[i]);
                interaction_step(r, dt*plf7_6_4_y[i], 0.);
            }
            break;
        default:
            break;
    }
}

static inline void reb_integrator_eos_postprocessor(struct reb_simulation* const r, double dt, enum REB_INTEGRATOR_EOS_TYPE type, void (*drift_step)(struct reb_simulation* const r, double a), void (*interaction_step)(struct reb_simulation* const r, double y, double v)){
    switch(type){
        case REB_INTEGRATOR_EOS_TYPE_PMLF6:
            for (int i=5;i>=0;i--){
                interaction_step(r, -dt*pmlf6_y[i], -dt*dt*dt*pmlf6_v[i]); 
                drift_step(r, -dt*pmlf6_z[i]);
            }
            break;
        case REB_INTEGRATOR_EOS_TYPE_PMLF4:
            for (int i=2;i>=0;i--){
                drift_step(r, -dt*pmlf4_z[i]);
                interaction_step(r, -dt*pmlf4_y[i], 0.); 
            }
            break;
        case REB_INTEGRATOR_EOS_TYPE_PLF7_6_4:
            for (int i=5;i>=0;i--){
                interaction_step(r, -dt*plf7_6_4_y[i], 0.);
                drift_step(r, -dt*plf7_6_4_z[i]);
            }
            break;
        default:
            break;
    }
}
static void reb_integrator_eos_drift_shell1(struct reb_simulation* const r, double dt){
    struct reb_particle* restrict const particles = r->particles;
    size_t N = r->N;
    for (size_t i=0;i<N;i++){  
        particles[i].x += dt*particles[i].vx;
        particles[i].y += dt*particles[i].vy;
        particles[i].z += dt*particles[i].vz;
    } 
}

static void reb_integrator_eos_drift_shell0(struct reb_simulation* const r, double _dt){
    struct reb_integrator_eos_state* const eos = r->integrator.state;
    const size_t n = eos->n;
    const double dt = _dt/n;
    reb_integrator_eos_preprocessor(r, dt, eos->phi1, reb_integrator_eos_drift_shell1, reb_integrator_eos_interaction_shell1);
    switch(eos->phi1){
        case REB_INTEGRATOR_EOS_TYPE_LF:
            reb_integrator_eos_drift_shell1(r, dt*0.5); 
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, dt, 0.);
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, dt);
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*0.5);
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF4:
            reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf4_a);
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, dt*2.*reb_integrator_leapfrog_lf4_a, 0.);
                reb_integrator_eos_drift_shell1(r, dt*(0.5-reb_integrator_leapfrog_lf4_a));
                reb_integrator_eos_interaction_shell1(r, dt*(1.-4.*reb_integrator_leapfrog_lf4_a), 0.);
                reb_integrator_eos_drift_shell1(r, dt*(0.5-reb_integrator_leapfrog_lf4_a));
                reb_integrator_eos_interaction_shell1(r, dt*2.*reb_integrator_leapfrog_lf4_a, 0.);
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, dt*2.*reb_integrator_leapfrog_lf4_a);
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf4_a);
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF6:
            reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf6_a[0]*0.5);
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[0], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[0]+reb_integrator_leapfrog_lf6_a[1])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[1], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[1]+reb_integrator_leapfrog_lf6_a[2])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[2], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[2]+reb_integrator_leapfrog_lf6_a[3])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[3], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[3]+reb_integrator_leapfrog_lf6_a[4])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[4], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[3]+reb_integrator_leapfrog_lf6_a[4])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[3], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[2]+reb_integrator_leapfrog_lf6_a[3])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[2], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[1]+reb_integrator_leapfrog_lf6_a[2])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[1], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf6_a[0]+reb_integrator_leapfrog_lf6_a[1])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf6_a[0], 0.);
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf6_a[0]);
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf6_a[0]*0.5);
            break; 
        case REB_INTEGRATOR_EOS_TYPE_LF8: 
            reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf8_a[0]*0.5);
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[0], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[0]+reb_integrator_leapfrog_lf8_a[1])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[1], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[1]+reb_integrator_leapfrog_lf8_a[2])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[2], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[2]+reb_integrator_leapfrog_lf8_a[3])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[3], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[3]+reb_integrator_leapfrog_lf8_a[4])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[4], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[4]+reb_integrator_leapfrog_lf8_a[5])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[5], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[5]+reb_integrator_leapfrog_lf8_a[6])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[6], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[6]+reb_integrator_leapfrog_lf8_a[7])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[7], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[7]+reb_integrator_leapfrog_lf8_a[8])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[8], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[7]+reb_integrator_leapfrog_lf8_a[8])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[7], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[6]+reb_integrator_leapfrog_lf8_a[7])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[6], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[5]+reb_integrator_leapfrog_lf8_a[6])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[5], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[4]+reb_integrator_leapfrog_lf8_a[5])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[4], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[3]+reb_integrator_leapfrog_lf8_a[4])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[3], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[2]+reb_integrator_leapfrog_lf8_a[3])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[2], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[1]+reb_integrator_leapfrog_lf8_a[2])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[1], 0.);
                reb_integrator_eos_drift_shell1(r, dt*(reb_integrator_leapfrog_lf8_a[0]+reb_integrator_leapfrog_lf8_a[1])*0.5);
                reb_integrator_eos_interaction_shell1(r, dt*reb_integrator_leapfrog_lf8_a[0], 0.);
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf8_a[0]);
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*reb_integrator_leapfrog_lf8_a[0]*0.5);
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF4_2: 
            reb_integrator_eos_drift_shell1(r, dt*lf4_2_a); 
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, dt*0.5, 0.); 
                reb_integrator_eos_drift_shell1(r, dt*(1.-2.*lf4_2_a));
                reb_integrator_eos_interaction_shell1(r, dt*0.5, 0.); 
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, 2.*dt*lf4_2_a); 
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*lf4_2_a); 
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF8_6_4:
            reb_integrator_eos_drift_shell1(r, dt*lf8_6_4_a[0]);   
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, lf8_6_4_b[0]*dt,0);
                reb_integrator_eos_drift_shell1(r, lf8_6_4_a[1]*dt);   
                reb_integrator_eos_interaction_shell1(r, lf8_6_4_b[1]*dt,0);
                reb_integrator_eos_drift_shell1(r, lf8_6_4_a[2]*dt);   
                reb_integrator_eos_interaction_shell1(r, lf8_6_4_b[2]*dt,0);
                reb_integrator_eos_drift_shell1(r, lf8_6_4_a[3]*dt);   
                reb_integrator_eos_interaction_shell1(r, lf8_6_4_b[3]*dt,0);
                reb_integrator_eos_drift_shell1(r, lf8_6_4_a[3]*dt);   
                reb_integrator_eos_interaction_shell1(r, lf8_6_4_b[2]*dt,0);
                reb_integrator_eos_drift_shell1(r, lf8_6_4_a[2]*dt);   
                reb_integrator_eos_interaction_shell1(r, lf8_6_4_b[1]*dt,0);
                reb_integrator_eos_drift_shell1(r, lf8_6_4_a[1]*dt);   
                reb_integrator_eos_interaction_shell1(r, lf8_6_4_b[0]*dt,0);
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, 2.*dt*lf8_6_4_a[0]);   
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*lf8_6_4_a[0]);   
            break;
        case REB_INTEGRATOR_EOS_TYPE_PMLF4:
            reb_integrator_eos_drift_shell1(r, dt*0.5); 
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, dt, dt*dt*dt/24.); 
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, dt);
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*0.5);
            break;
        case REB_INTEGRATOR_EOS_TYPE_PMLF6:
            reb_integrator_eos_drift_shell1(r, dt*pmlf6_a[0]); 
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, dt*pmlf6_b[0], dt*dt*dt*pmlf6_c[0]); 
                reb_integrator_eos_drift_shell1(r, dt*pmlf6_a[1]);
                reb_integrator_eos_interaction_shell1(r, dt*pmlf6_b[1], dt*dt*dt*pmlf6_c[1]); 
                reb_integrator_eos_drift_shell1(r, dt*pmlf6_a[1]);
                reb_integrator_eos_interaction_shell1(r, dt*pmlf6_b[0], dt*dt*dt*pmlf6_c[0]);
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, 2.*dt*pmlf6_a[0]);
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*pmlf6_a[0]);
            break;
        case REB_INTEGRATOR_EOS_TYPE_PLF7_6_4:
            reb_integrator_eos_drift_shell1(r, dt*plf7_6_4_a[0]);   
            for (size_t i=0;i<n;i++){
                reb_integrator_eos_interaction_shell1(r, plf7_6_4_b[0]*dt,0);
                reb_integrator_eos_drift_shell1(r, plf7_6_4_a[1]*dt);   
                reb_integrator_eos_interaction_shell1(r, plf7_6_4_b[1]*dt,0);
                reb_integrator_eos_drift_shell1(r, plf7_6_4_a[1]*dt);   
                reb_integrator_eos_interaction_shell1(r, plf7_6_4_b[0]*dt,0);
                if (i<n-1){
                    reb_integrator_eos_drift_shell1(r, 2.*dt*plf7_6_4_a[0]);   
                }
            }
            reb_integrator_eos_drift_shell1(r, dt*plf7_6_4_a[0]);   
            break;
    }
    reb_integrator_eos_postprocessor(r, dt, eos->phi1, reb_integrator_eos_drift_shell1, reb_integrator_eos_interaction_shell1);
}

void reb_integrator_eos_step(struct reb_simulation* r, void* state){
    if (r->gravity != REB_GRAVITY_BASIC){
        reb_simulation_warning(r,"EOS only supports the BASIC gravity routine.");
    }
    if (r->N_var){
        // TODO Variational equations no longer working.
        reb_simulation_warning(r, "Variational particles/MEGNO in EOS no longer supported since REBOUND version 5.");
    }
    r->gravity = REB_GRAVITY_NONE;

    struct reb_integrator_eos_state* const eos = state;
    const double dt = r->dt;

    double dtfac = 1.;
    if (r->is_synchronized){
        reb_integrator_eos_preprocessor(r, r->dt, eos->phi0, reb_integrator_eos_drift_shell0, reb_integrator_eos_interaction_shell0);
    }else{
        dtfac = 2.;
    }
    switch(eos->phi0){
        case REB_INTEGRATOR_EOS_TYPE_LF:
            reb_integrator_eos_drift_shell0(r, dt*0.5*dtfac);
            reb_integrator_eos_interaction_shell0(r, dt, 0.);
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF4:
            reb_integrator_eos_drift_shell0(r, dt*reb_integrator_leapfrog_lf4_a*dtfac);
            reb_integrator_eos_interaction_shell0(r, dt*2.*reb_integrator_leapfrog_lf4_a, 0.);
            reb_integrator_eos_drift_shell0(r, dt*(0.5-reb_integrator_leapfrog_lf4_a));
            reb_integrator_eos_interaction_shell0(r, dt*(1.-4.*reb_integrator_leapfrog_lf4_a), 0.);
            reb_integrator_eos_drift_shell0(r, dt*(0.5-reb_integrator_leapfrog_lf4_a));
            reb_integrator_eos_interaction_shell0(r, dt*2.*reb_integrator_leapfrog_lf4_a, 0.);
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF6:
            reb_integrator_eos_drift_shell0(r, dt*reb_integrator_leapfrog_lf6_a[0]*0.5*dtfac);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[0], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[0]+reb_integrator_leapfrog_lf6_a[1])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[1], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[1]+reb_integrator_leapfrog_lf6_a[2])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[2], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[2]+reb_integrator_leapfrog_lf6_a[3])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[3], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[3]+reb_integrator_leapfrog_lf6_a[4])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[4], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[3]+reb_integrator_leapfrog_lf6_a[4])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[3], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[2]+reb_integrator_leapfrog_lf6_a[3])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[2], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[1]+reb_integrator_leapfrog_lf6_a[2])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[1], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf6_a[0]+reb_integrator_leapfrog_lf6_a[1])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf6_a[0], 0.);
            break; 
        case REB_INTEGRATOR_EOS_TYPE_LF8: 
            reb_integrator_eos_drift_shell0(r, dt*reb_integrator_leapfrog_lf8_a[0]*0.5*dtfac);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[0], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[0]+reb_integrator_leapfrog_lf8_a[1])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[1], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[1]+reb_integrator_leapfrog_lf8_a[2])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[2], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[2]+reb_integrator_leapfrog_lf8_a[3])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[3], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[3]+reb_integrator_leapfrog_lf8_a[4])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[4], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[4]+reb_integrator_leapfrog_lf8_a[5])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[5], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[5]+reb_integrator_leapfrog_lf8_a[6])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[6], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[6]+reb_integrator_leapfrog_lf8_a[7])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[7], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[7]+reb_integrator_leapfrog_lf8_a[8])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[8], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[7]+reb_integrator_leapfrog_lf8_a[8])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[7], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[6]+reb_integrator_leapfrog_lf8_a[7])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[6], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[5]+reb_integrator_leapfrog_lf8_a[6])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[5], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[4]+reb_integrator_leapfrog_lf8_a[5])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[4], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[3]+reb_integrator_leapfrog_lf8_a[4])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[3], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[2]+reb_integrator_leapfrog_lf8_a[3])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[2], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[1]+reb_integrator_leapfrog_lf8_a[2])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[1], 0.);
            reb_integrator_eos_drift_shell0(r, dt*(reb_integrator_leapfrog_lf8_a[0]+reb_integrator_leapfrog_lf8_a[1])*0.5);
            reb_integrator_eos_interaction_shell0(r, dt*reb_integrator_leapfrog_lf8_a[0], 0.);
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF4_2: 
            reb_integrator_eos_drift_shell0(r, dt*lf4_2_a*dtfac); 
            reb_integrator_eos_interaction_shell0(r, dt*0.5, 0.); 
            reb_integrator_eos_drift_shell0(r, dt*(1.-2.*lf4_2_a));
            reb_integrator_eos_interaction_shell0(r, dt*0.5, 0.); 
            break;
        case REB_INTEGRATOR_EOS_TYPE_LF8_6_4:
            reb_integrator_eos_drift_shell0(r, dt*lf8_6_4_a[0]*dtfac);   
            reb_integrator_eos_interaction_shell0(r, lf8_6_4_b[0]*dt,0);
            reb_integrator_eos_drift_shell0(r, lf8_6_4_a[1]*dt);   
            reb_integrator_eos_interaction_shell0(r, lf8_6_4_b[1]*dt,0);
            reb_integrator_eos_drift_shell0(r, lf8_6_4_a[2]*dt);   
            reb_integrator_eos_interaction_shell0(r, lf8_6_4_b[2]*dt,0);
            reb_integrator_eos_drift_shell0(r, lf8_6_4_a[3]*dt);   
            reb_integrator_eos_interaction_shell0(r, lf8_6_4_b[3]*dt,0);
            reb_integrator_eos_drift_shell0(r, lf8_6_4_a[3]*dt);   
            reb_integrator_eos_interaction_shell0(r, lf8_6_4_b[2]*dt,0);
            reb_integrator_eos_drift_shell0(r, lf8_6_4_a[2]*dt);   
            reb_integrator_eos_interaction_shell0(r, lf8_6_4_b[1]*dt,0);
            reb_integrator_eos_drift_shell0(r, lf8_6_4_a[1]*dt);   
            reb_integrator_eos_interaction_shell0(r, lf8_6_4_b[0]*dt,0);
            break;
        case REB_INTEGRATOR_EOS_TYPE_PMLF4:
            reb_integrator_eos_drift_shell0(r, dt*0.5*dtfac); 
            reb_integrator_eos_interaction_shell0(r, dt, dt*dt*dt/24.); 
            break;
        case REB_INTEGRATOR_EOS_TYPE_PMLF6:
            reb_integrator_eos_drift_shell0(r, dt*pmlf6_a[0]*dtfac); 
            reb_integrator_eos_interaction_shell0(r, dt*pmlf6_b[0], dt*dt*dt*pmlf6_c[0]); 
            reb_integrator_eos_drift_shell0(r, dt*pmlf6_a[1]);
            reb_integrator_eos_interaction_shell0(r, dt*pmlf6_b[1], dt*dt*dt*pmlf6_c[1]); 
            reb_integrator_eos_drift_shell0(r, dt*pmlf6_a[1]);
            reb_integrator_eos_interaction_shell0(r, dt*pmlf6_b[0], dt*dt*dt*pmlf6_c[0]);
            break;
        case REB_INTEGRATOR_EOS_TYPE_PLF7_6_4:
            reb_integrator_eos_drift_shell0(r, dt*plf7_6_4_a[0]*dtfac);   
            reb_integrator_eos_interaction_shell0(r, plf7_6_4_b[0]*dt,0);
            reb_integrator_eos_drift_shell0(r, plf7_6_4_a[1]*dt);   
            reb_integrator_eos_interaction_shell0(r, plf7_6_4_b[1]*dt,0);
            reb_integrator_eos_drift_shell0(r, plf7_6_4_a[1]*dt);   
            reb_integrator_eos_interaction_shell0(r, plf7_6_4_b[0]*dt,0);
            break;
    }

    r->is_synchronized = 0;
    if (eos->safe_mode){
        reb_integrator_eos_synchronize(r, state);
    }

    r->t+=r->dt;
    r->dt_last_done = r->dt;
}

void reb_integrator_eos_synchronize(struct reb_simulation* r, void* state){
    struct reb_integrator_eos_state* const eos = state;
    const double dt = r->dt;
    if (r->is_synchronized == 0){
        switch(eos->phi0){
            case REB_INTEGRATOR_EOS_TYPE_PMLF4:
            case REB_INTEGRATOR_EOS_TYPE_LF:
                reb_integrator_eos_drift_shell0(r, dt*0.5); 
                break;
            case REB_INTEGRATOR_EOS_TYPE_PMLF6:
                reb_integrator_eos_drift_shell0(r, dt*pmlf6_a[0]); 
                break;
            case REB_INTEGRATOR_EOS_TYPE_LF4:
                reb_integrator_eos_drift_shell0(r, dt*reb_integrator_leapfrog_lf4_a);
                break;
            case REB_INTEGRATOR_EOS_TYPE_LF4_2:
                reb_integrator_eos_drift_shell0(r, dt*lf4_2_a); 
                break;
            case REB_INTEGRATOR_EOS_TYPE_PLF7_6_4:
                reb_integrator_eos_drift_shell0(r, dt*plf7_6_4_a[0]);   
                break;
            case REB_INTEGRATOR_EOS_TYPE_LF8_6_4:
                reb_integrator_eos_drift_shell0(r, dt*lf8_6_4_a[0]);   
                break;
            case REB_INTEGRATOR_EOS_TYPE_LF6:
                reb_integrator_eos_drift_shell0(r, dt*reb_integrator_leapfrog_lf6_a[0]*0.5);
                break;
            case REB_INTEGRATOR_EOS_TYPE_LF8: 
                reb_integrator_eos_drift_shell0(r, dt*reb_integrator_leapfrog_lf8_a[0]*0.5);
                break;
        }
        reb_integrator_eos_postprocessor(r, r->dt, eos->phi0, reb_integrator_eos_drift_shell0, reb_integrator_eos_interaction_shell0);
        r->is_synchronized = 1;
    }
}