From 2f4b841d38a365500366177c520a53d1df7ebc90 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Thu, 25 Jan 2024 16:54:21 -0500 Subject: [PATCH] fixes collision test --- src/integrator_trace.c | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 239db2d59..d0cd76731 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -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; iencounter_N; i++){ const int mi = ri_trace->encounter_map[i]; const struct reb_particle p = r->particles[mi]; @@ -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; @@ -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; @@ -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.);