Skip to content

Commit 096fbdc

Browse files
committed
added fix for backwards timestep
1 parent 0f00a2b commit 096fbdc

File tree

1 file changed

+13
-7
lines changed

1 file changed

+13
-7
lines changed

galpy/util/wez_ias15.c

+13-7
Original file line numberDiff line numberDiff line change
@@ -252,6 +252,12 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
252252
}
253253

254254
long ndt= (long) (init_dt/dt);
255+
int timestep_sign;
256+
if(init_dt > 0){
257+
timestep_sign = 1;
258+
} else {
259+
timestep_sign = -1;
260+
}
255261

256262
double to= *t;
257263
// Handle KeyboardInterrupt gracefully
@@ -277,13 +283,13 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
277283
}
278284

279285
//WHILE THERE IS TIME REMAINING, INTEGRATE A DYNAMIC TIMESTEP FORWARD AND SUBTRACT FROM THE TIME REMAINING
280-
double time_remaining = init_dt;
281-
286+
double time_remaining = fabs(init_dt);
287+
282288
while(time_remaining > 0) {
283289
double to_temp;
284-
double dt_temp; //keeping an extra dt variable in case I want to change my dynamic strategy, I take the min of time remaining and dt, so not sure how taking the time remaining might effect the next dynamic time step
285-
if (time_remaining < dt){
286-
dt_temp = time_remaining;
290+
double dt_temp;
291+
if (time_remaining < fabs(dt)){
292+
dt_temp = timestep_sign * time_remaining;
287293
} else {
288294
dt_temp = dt;
289295
}
@@ -359,12 +365,12 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
359365

360366
double dt_required = dt_temp * pow(precision_parameter / (max_B6/max_a), 1.0/7.0);
361367

362-
if(dt_temp > dt_required){
368+
if(fabs(dt_temp) > fabs(dt_required)){
363369
//rejected, try again with dt required
364370
dt = dt_required;
365371
} else {
366372
//accepted, update position/velocity and do next timestep with dt required
367-
time_remaining -= dt; //will eventually get negative as we stepped forward the minimum of dt and time_remaining
373+
time_remaining -= fabs(dt); //will eventually get negative as we stepped forward the minimum of dt and time_remaining
368374
to = to_temp;
369375

370376
update_position(x, x, v, dim, 1, dt_temp, Fs, Bs);

0 commit comments

Comments
 (0)