@@ -208,18 +208,18 @@ IAS15 integrator
208
208
void wez_ias15 (void (* func )(double t , double * q , double * a , int nargs , struct potentialArg * potentialArgs ),
209
209
int dim ,
210
210
double * yo ,
211
- int nt ,
212
- double dt ,
211
+ int nt ,
212
+ double dt ,
213
213
double * t ,
214
- int nargs ,
214
+ int nargs ,
215
215
struct potentialArg * potentialArgs ,
216
- double rtol ,
216
+ double rtol ,
217
217
double atol ,
218
- double * result ,
218
+ double * result ,
219
219
int * err ){
220
220
//Declare and initialize
221
221
double * x = (double * ) malloc ( dim * sizeof (double ) );
222
- double * v = (double * ) malloc ( dim * sizeof (double ) );
222
+ double * v = (double * ) malloc ( dim * sizeof (double ) );
223
223
double * a = (double * ) malloc ( dim * sizeof (double ) );
224
224
double * xs = (double * ) malloc ( dim * sizeof (double ) ); //x substep
225
225
@@ -238,7 +238,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
238
238
239
239
x -= dim ;
240
240
v -= dim ;
241
-
241
+
242
242
double diff_G ;
243
243
244
244
save_ias15 (dim , x , v , result );
@@ -248,9 +248,9 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
248
248
//Estimate necessary stepsize, use the returned time interval if the user does not provide
249
249
double init_dt = (* (t + 1 ))- (* t );
250
250
if ( dt == -9999.99 ) {
251
- dt = init_dt ;
251
+ dt = init_dt ;
252
252
}
253
-
253
+
254
254
long ndt = (long ) (init_dt /dt );
255
255
int timestep_sign ;
256
256
if (init_dt > 0 ){
@@ -287,7 +287,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
287
287
288
288
while (time_remaining > 0 ) {
289
289
double to_temp ;
290
- double dt_temp ;
290
+ double dt_temp ;
291
291
if (time_remaining < fabs (dt )){
292
292
dt_temp = timestep_sign * time_remaining ;
293
293
} else {
@@ -304,7 +304,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
304
304
305
305
int iterations = 0 ;
306
306
double integrator_error = integrator_error_threshold + 1 ; //init value above the threshold
307
- while (true){
307
+ while (true){
308
308
if (iterations == 12 ){
309
309
//TODO: spit the dummy
310
310
//printf("Iterations took over 12, break iteration, break iteration loop\n");
@@ -316,7 +316,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
316
316
317
317
//at start of each step reset xs
318
318
for (int l = 0 ; l < dim ; l ++ ) * (xs + l )= * (x + l );
319
-
319
+
320
320
double max_delta_B6 = 0.0 ; //also = max delta G
321
321
double max_a = 0.0 ;
322
322
for (int k = 1 ; k < (order + 1 ); k ++ ){
@@ -349,7 +349,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
349
349
iterations += 1 ;
350
350
}
351
351
352
- //global error strategy for timestep
352
+ //global error strategy for timestep
353
353
double max_B6 = 0.0 ;
354
354
double max_a = 0.0 ;
355
355
for (int i = 0 ; i < dim ; i ++ ){
@@ -367,7 +367,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
367
367
368
368
if (fabs (dt_temp ) > fabs (dt_required )){
369
369
//rejected, try again with dt required
370
- dt = dt_required ;
370
+ dt = dt_required ;
371
371
} else {
372
372
//accepted, update position/velocity and do next timestep with dt required
373
373
time_remaining -= fabs (dt ); //will eventually get negative as we stepped forward the minimum of dt and time_remaining
@@ -376,9 +376,9 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
376
376
update_position (x , x , v , dim , 1 , dt_temp , Fs , Bs );
377
377
update_velocity (v , v , dim , 1 , dt_temp , Fs , Bs );
378
378
next_sequence_Bs (1 , Bs , Es , BDs , dim );
379
-
379
+
380
380
//finally, update the next step
381
- dt = dt_required ;
381
+ dt = dt_required ;
382
382
}
383
383
}
384
384
@@ -405,36 +405,36 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
405
405
}
406
406
407
407
void update_velocity (double * v , double * v1 , int dim , double h_n , double T , double * Fs , double * Bs ){
408
- for (int ii = 0 ; ii < dim ; ii ++ ){
409
- * (v + ii ) = * (v1 + ii ) +
410
- (h_n * T *
411
- (Fs [ii * (order + 1 ) + 0 ] + h_n *
412
- (Bs [ii * order + 0 ]/2 + h_n *
413
- (Bs [ii * order + 1 ]/3 + h_n *
414
- (Bs [ii * order + 2 ]/4 + h_n *
415
- (Bs [ii * order + 3 ]/5 + h_n *
416
- (Bs [ii * order + 4 ]/6 + h_n *
417
- (Bs [ii * order + 5 ]/7 + h_n *
418
- (Bs [ii * order + 6 ]/8 //+ h_n *
408
+ for (int ii = 0 ; ii < dim ; ii ++ ){
409
+ * (v + ii ) = * (v1 + ii ) +
410
+ (h_n * T *
411
+ (Fs [ii * (order + 1 ) + 0 ] + h_n *
412
+ (Bs [ii * order + 0 ]/2 + h_n *
413
+ (Bs [ii * order + 1 ]/3 + h_n *
414
+ (Bs [ii * order + 2 ]/4 + h_n *
415
+ (Bs [ii * order + 3 ]/5 + h_n *
416
+ (Bs [ii * order + 4 ]/6 + h_n *
417
+ (Bs [ii * order + 5 ]/7 + h_n *
418
+ (Bs [ii * order + 6 ]/8 //+ h_n *
419
419
//(Bs[ii * order + 7]/9
420
420
)))))))));
421
421
}
422
422
}
423
-
423
+
424
424
425
425
void update_position (double * y , double * y1 , double * v , int dim , double h_n , double T , double * Fs , double * Bs ){
426
- for (int ii = 0 ; ii < dim ; ii ++ ){
427
- * (y + ii ) = * (y1 + ii ) +
428
- (* (v + ii ) * h_n * T ) +
429
- ((h_n * T )* (h_n * T )) *
430
- (Fs [ii * (order + 1 ) + 0 ]/2 + h_n *
431
- (Bs [ii * order + 0 ]/6 + h_n *
432
- (Bs [ii * order + 1 ]/12 + h_n *
433
- (Bs [ii * order + 2 ]/20 + h_n *
434
- (Bs [ii * order + 3 ]/30 + h_n *
435
- (Bs [ii * order + 4 ]/42 + h_n *
436
- (Bs [ii * order + 5 ]/56 + h_n *
437
- (Bs [ii * order + 6 ]/72 //+ h_n *
426
+ for (int ii = 0 ; ii < dim ; ii ++ ){
427
+ * (y + ii ) = * (y1 + ii ) +
428
+ (* (v + ii ) * h_n * T ) +
429
+ ((h_n * T )* (h_n * T )) *
430
+ (Fs [ii * (order + 1 ) + 0 ]/2 + h_n *
431
+ (Bs [ii * order + 0 ]/6 + h_n *
432
+ (Bs [ii * order + 1 ]/12 + h_n *
433
+ (Bs [ii * order + 2 ]/20 + h_n *
434
+ (Bs [ii * order + 3 ]/30 + h_n *
435
+ (Bs [ii * order + 4 ]/42 + h_n *
436
+ (Bs [ii * order + 5 ]/56 + h_n *
437
+ (Bs [ii * order + 6 ]/72 //+ h_n *
438
438
//(Bs[ii * order + 7]/90)
439
439
))))))));
440
440
}
@@ -532,7 +532,7 @@ void update_Bs_from_Gs(int current_truncation_order, int i, double * Bs, double
532
532
Bs [j + 2 ] += diff_G * c_4_2 ;
533
533
Bs [j + 3 ] += diff_G * c_4_3 ;
534
534
Bs [j + 4 ] += diff_G ;
535
- }
535
+ }
536
536
if (current_truncation_order == 6 ){
537
537
Bs [j ] += diff_G * c_5_0 ;
538
538
Bs [j + 1 ] += diff_G * c_5_1 ;
@@ -550,4 +550,4 @@ void update_Bs_from_Gs(int current_truncation_order, int i, double * Bs, double
550
550
Bs [j + 5 ] += diff_G * c_6_5 ;
551
551
Bs [j + 6 ] += diff_G ;
552
552
}
553
- }
553
+ }
0 commit comments