mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-05-24 00:35:19 -04:00
Code Issues Packages Projects Releases Wiki Activity

Pluto integrator: fixed a small math mistake.

Browse Source 
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
This commit is contained in:
Don Cross
2020-08-18 21:59:37 -04:00
parent f279cdd9b5
commit b6087c9b14
2 changed files with 6 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
generate/template/astronomy.c
View File

@@ -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)

8
source/c/astronomy.c
View File

@@ -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)