I'm writing code to simulate planetary and lunar orbits in a heliocentric reference plane. I'm using E.M. Standish's Keplerian Elements for Approximate Positions of the Major Planets as my reference and I think it is a straight path to simulate planetary motions from the equations and data given there.
What I haven't found is a similar reference that has data for the moons of those planets. I'm assuming that I can use those same formulae to compute the planet-centric positions of its moons and then translate them to heliocentric coordinates by adding the planet's heliocentric location.
Additional information: The document referred to presents one (time unvarying) set of Keplerian elements which can be applied to the supplied equations to simulate planetary orbits. The Keplerian elements have been best fit to the actual orbits over a given time range. For some of the planets the elements have a time dependent component which has also been fit to an analytic formula. The entire scheme is thus governed by a small set of constants driving a formula with time as the only input variable.