2
$\begingroup$

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.

10 year plot

40 year plot

100 year plot

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

Improve this question
$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$
    – uhoh
    Commented Apr 4, 2018 at 16:05
  • $\begingroup$ 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. $\endgroup$
    – uhoh
    Commented Apr 4, 2018 at 16:14
  • 1
    $\begingroup$ 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? $\endgroup$
    – Tarius
    Commented Apr 4, 2018 at 16:29
  • 1
    $\begingroup$ 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. $\endgroup$
    – uhoh
    Commented Apr 4, 2018 at 16:53
  • 1
    $\begingroup$ @uhoh I've just updated the azimuth as per your instruction. Thanks ever so much for the help, works a treat. $\endgroup$
    – Tarius
    Commented Apr 5, 2018 at 19:35

1 Answer 1

Reset to default
1
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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")? $\endgroup$
    – uhoh
    Commented Apr 6, 2018 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.