You say you want to create a simulation of actual planetary movement, but I think you have to first decide how accurate you want it, and for what time range. I'm going to outline some different options that vary starting from the least accurate to most accurate:
1) Assume elliptical plantary movement with constant angular velocity. Then all you need to do is plug in 1 parametric variable for time into three eqns for x, y, and z.
2) Same as above but using Kepler's law to take into account the angular velocity differences depending on the distance from barycenter.
3) Use an analytical solution such as VSOP87. This would be accurate to about an arcsecond. The VSOP87 solution doesn't give Pluto
4) Attempt "solve" an n-body problem as suggested above. I wouldn't recommend it.
5) Use an integrated solution (basically the n-body problem has already been done). DE200, DE403 are probably your best bet. This would give sub arcsecond resolution.
|