From b6087c9b145156e32da571f7ee751aa702ce8071 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Tue, 18 Aug 2020 21:59:37 -0400 Subject: [PATCH] Pluto integrator: fixed a small math mistake. I didn't implement the math I had derived on paper exactly right. There was one place where I meant to have (dt^3 / 6), but I had (dt^2 / 6). This has only a tiny effect on the error. I know this can work better, so I'm still searching for the real problem. C PlutoCheck: 2049-12-19T12:00:00.000Z = 18250.000000 UT = 18250.001076 TT C PlutoCheck: calc pos = [ 37.3682426267669996, -10.0216432448322834, -14.3724661626442582] C PlutoCheck: ref pos = [ 37.4377303529113306, -10.2466292445590774, -14.4773101309873091] C PlutoCheck: del pos = [ -0.0694877261443310, 0.2249859997267940, 0.1048439683430509] C PlutoCheck: error = 2.577586e-01 --- generate/template/astronomy.c | 8 +++----- source/c/astronomy.c | 8 +++----- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 5bb0e096..af11d8fa 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -2098,15 +2098,13 @@ body_grav_calc_t GravSim( /* out: [pos, vel, acc] of the simulated bod /* Calculate where the major bodies (Sun, Jupiter...Neptune) will be at the next time step. */ MajorBodyBary(bary2, tt2); - calc2.r = approx_pos; - /* Calculate acceleration experienced by small body at approximate next location. */ - next_acc = SmallBodyAcceleration(calc2.r, bary2); + next_acc = SmallBodyAcceleration(approx_pos, bary2); /* Assume each component of the acceleration vector ramps linearly over the interval. */ /* Integrating over the interval [tt1, tt2], we get expressions for r2, v2. */ delta_acc = VecSub(next_acc, calc1->a); - calc2.r = VecAdd(approx_pos, VecMul(dt*dt/6, delta_acc)); + calc2.r = VecAdd(approx_pos, VecMul(dt*dt*dt/6, delta_acc)); calc2.v = VecAdd(calc1->v, VecMul(dt, calc1->a)); VecIncr(&calc2.v, VecMul(dt/2, delta_acc)); calc2.a = SmallBodyAcceleration(calc2.r, bary2); @@ -2116,7 +2114,7 @@ body_grav_calc_t GravSim( /* out: [pos, vel, acc] of the simulated bod } -double CalcPlutoDeltaTime = 10.0; +double CalcPlutoDeltaTime = 1.0; static astro_vector_t CalcPluto(astro_time_t time) diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 0724871c..73a9fedb 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -3149,15 +3149,13 @@ body_grav_calc_t GravSim( /* out: [pos, vel, acc] of the simulated bod /* Calculate where the major bodies (Sun, Jupiter...Neptune) will be at the next time step. */ MajorBodyBary(bary2, tt2); - calc2.r = approx_pos; - /* Calculate acceleration experienced by small body at approximate next location. */ - next_acc = SmallBodyAcceleration(calc2.r, bary2); + next_acc = SmallBodyAcceleration(approx_pos, bary2); /* Assume each component of the acceleration vector ramps linearly over the interval. */ /* Integrating over the interval [tt1, tt2], we get expressions for r2, v2. */ delta_acc = VecSub(next_acc, calc1->a); - calc2.r = VecAdd(approx_pos, VecMul(dt*dt/6, delta_acc)); + calc2.r = VecAdd(approx_pos, VecMul(dt*dt*dt/6, delta_acc)); calc2.v = VecAdd(calc1->v, VecMul(dt, calc1->a)); VecIncr(&calc2.v, VecMul(dt/2, delta_acc)); calc2.a = SmallBodyAcceleration(calc2.r, bary2); @@ -3167,7 +3165,7 @@ body_grav_calc_t GravSim( /* out: [pos, vel, acc] of the simulated bod } -double CalcPlutoDeltaTime = 10.0; +double CalcPlutoDeltaTime = 1.0; static astro_vector_t CalcPluto(astro_time_t time)