Skip to content

Commit

Permalink
fixes collision test
Browse files Browse the repository at this point in the history
  • Loading branch information
Tiger Lu authored and Tiger Lu committed Jan 25, 2024
1 parent 29f3563 commit 2f4b841
Showing 1 changed file with 12 additions and 10 deletions.
22 changes: 12 additions & 10 deletions src/integrator_trace.c
Original file line number Diff line number Diff line change
Expand Up @@ -351,26 +351,28 @@ void reb_integrator_trace_bs_step(struct reb_simulation* const r, double dt){
// run
const double old_dt = r->dt;
const double old_t = r->t;
double t = r->t;
//double t = r->t;
const double t_needed = r->t + dt;
reb_integrator_bs_reset(r);

struct reb_ode* nbody_ode = reb_ode_create(r, ri_trace->encounter_N*3*2);
nbody_ode->derivatives = reb_integrator_trace_nbody_derivatives;
nbody_ode->needs_nbody = 0;
double* const y = nbody_ode->y;
double* y;

while(t < t_needed && fabs(dt/old_dt)>1e-14 ){
while(r->t < t_needed && fabs(dt/old_dt)>1e-14 ){
y = nbody_ode->y;

// In case of overshoot
if (t + dt > t_needed){
dt = t_needed - t;
if (r->t + dt > t_needed){
dt = t_needed - r->t;
}

struct reb_particle star = r->particles[0]; // backup velocity
r->particles[0].vx = 0; // star does not move in dh
r->particles[0].vy = 0;
r->particles[0].vz = 0;

for (unsigned int i=0; i<ri_trace->encounter_N; i++){
const int mi = ri_trace->encounter_map[i];
const struct reb_particle p = r->particles[mi];
Expand All @@ -384,11 +386,11 @@ void reb_integrator_trace_bs_step(struct reb_simulation* const r, double dt){

int success = reb_integrator_bs_step(r, dt);
if (success){
t += dt;
r->t += dt;
}
dt = r->ri_bs.dt_proposed;
reb_integrator_trace_update_particles(r, nbody_ode->y);

r->particles[0].vx = star.vx; // restore every timestep for collisions
r->particles[0].vy = star.vy;
r->particles[0].vz = star.vz;
Expand Down Expand Up @@ -430,7 +432,7 @@ void reb_integrator_trace_bs_step(struct reb_simulation* const r, double dt){
r->particles[mi] = ri_trace->particles_backup_kepler[mi];
}
}

reb_ode_free(nbody_ode);

r->t = old_t;
Expand Down Expand Up @@ -630,8 +632,8 @@ void reb_integrator_trace_part2(struct reb_simulation* const r){

reb_integrator_trace_interaction_step(r, r->dt/2.);
reb_integrator_trace_jump_step(r, r->dt/2.);
reb_integrator_trace_kepler_step(r, r->dt); // always accept this
reb_integrator_trace_com_step(r,r->dt);
reb_integrator_trace_kepler_step(r, r->dt); // always accept this
reb_integrator_trace_jump_step(r, r->dt/2.);
reb_integrator_trace_interaction_step(r, r->dt/2.);

Expand Down

0 comments on commit 2f4b841

Please sign in to comment.