Plotting error values between simulated and JPL data

Question

I'm attempting to compare values between my simulated data for the solar system and the official JPL data, to find errors for the position, velocity, inclination and azimuth.

The simulated data was created by using the initial values from the JPL system for all 10 major bodies (Sun, Planets & Pluto) and then calculating the gravitational perturbations on each body by the other 9 bodies and then combining that with the relativistic effects from the Sun.

The following equations were then used to calculate the position, velocity, inclination and azimuth respectively:

$$r_{p,v} = \sqrt(x^2 +y^2 +z^2) $$ $$\theta = acos(z/r_p)$$ $$\phi = atan2(x,y) $$

The error was then calculated by simply taking the absolute difference between each point from the simulated and JPL data.

I then plotted the points as follows for Earth in 10 day iterations, for 10, 40 and 100 years respectively. Position is in AU, velocity is in AU/day and the time axis is in days.

My question thus is: Do my error plots look reasonable and are there any better methods which I could use to either calculate or represent my data?

Thanks, for all the help in advance.

Edit: Updated graphs with fixed azimuth. Thanks to @uhoh

You may find some useful information in the following questions and their answers: How to calculate the planets and moons beyond Newtons's gravitational force?, and also Calculating the planets and moons based on Newtons's gravitational force Can you mention the units? Position in AU perhaps? velocity in ?? Time in days? — uhoh, Apr 04 '18 at 16:05
The 2π spikes in azimuth are probably not meaningful, just aesthetic. Some modulus housekeeping might fix them. Also, you did not mention the integrator you are using. I've had some helpful answers here: What does “symplectic” mean in reference to numerical integrators, and does SciPy's odeint use them? You didn't mention the Moon! Leaving it out is a huge effect! You can ignore Pluto but don't ignore the Moon. — uhoh, Apr 04 '18 at 16:14
My apologies for not mentioning: the Earth + Moon are counted as one object in this simulation, position is in AU, velocity is in AU/day and time is in days. Could you please explain what housekeeping rules I should use, to fix the 2π spikes? — Tarius, Apr 04 '18 at 16:29
okay thanks for the quick reply, that's going to be very close to doing the Earth and Moon separately, but it's not exact so at some point you'll need to include them as two distinct bodies. For the housekeeping, if $\Delta \phi=\phi_{yours}-\phi_{Horizons}$, then just use $\mod( (\Delta \phi + \pi), \ 2\pi) - \pi$ to fold $\pm 2\pi$ back to 0. I think that's it, or something very similar. Once you do that, you might consider updating the plots, since the spikes are making the interesting part hard to see right now. — uhoh, Apr 04 '18 at 16:53
Also, you might consider adding at least Ceres, and perhaps Pallas and Vesta as well, like I did here. — uhoh, Apr 04 '18 at 17:02
@uhoh I've just updated the azimuth as per your instruction. Thanks ever so much for the help, works a treat. — Tarius, Apr 05 '18 at 19:35

score 1 · Answer 1 · answered Apr 04 '18 at 19:50

Are there any better methods which I could use to either calculate or represent my data?

I'd suggest looking at position errors only, in terms of radial, cross track, and along track errors. Use the JPL position $\boldsymbol r$ and velocity $\boldsymbol v$ to define the equivalent of a local vertical, local horizontal frame: $$\begin{align} \boldsymbol h &= \boldsymbol r \times \boldsymbol v \\ \hat{\boldsymbol r} &= \boldsymbol r\,/\,||\boldsymbol r|| \\ \hat{\boldsymbol h} &= \boldsymbol h\,/\,||\boldsymbol h|| \\ \hat{\boldsymbol u} &= \hat{\boldsymbol r} \times \hat{\boldsymbol h} \end{align}$$ Transform the position error vector to this frame. I suspect you'll find that the along track errors eventually dominate over the others. You can also display your velocity error in the same frame.

Aside:

The simulated data was created by using the initial values from the JPL system for all 10 major bodies.

What you really want is the epoch data JPL used to compute its ephemerides. Those data unfortunately are hard to come by.

The JPL ephemerides are formed by repeatedly propagating those epoch values over time, computing errors against observations, and computing a new set of epoch values so as to minimize the computed errors. Once a good set of epoch values are found, they are once again propagated over time so as to form the basis for the Chebyshev polynomial coefficients that form the kernels. As far as I know, the Chebyshev fit is for position only.

This means that using those Chebyshev coefficients to form an initial position and velocity for your propagator means your propagator would not fair as well as JPL's propagator would against observations — even if your propagator is more accurate than is theirs! (That your propagator is better than theirs is almost certainly is not the case; they have invested several person years of effort in making their propagation extremely accurate.)

This is great stuff t to know! It would make sense that the velocities are just the analytical derivatives of the Chebyshev polynomials for position, rather than independently fitted. Also, I'm getting $\mathbf{\hat{u}}$ pointing roughly retrograde, is that the way it is supposed to be? For example if $\mathbf{r} = \mathbf{\hat{x}}$ and $\mathbf{v} = \mathbf{\hat{y}}$, then $\mathbf{h} = \mathbf{\hat{z}}$ ("up"), but $\mathbf{\hat{u}} = -\mathbf{\hat{y}}$ ("backwards")? — uhoh, Apr 06 '18 at 14:17

Plotting error values between simulated and JPL data

1 Answers1