http://www.usm.unimuenchen.de/ISMSPP/index.htm
https://github.com/sotos82/SolarSystemSimulatorGame
http://www.willbell.com/MATH/mc1.htm
https://en.wikibooks.org/wiki/Celestia/Trajectories
http://www.stjarnhimlen.se/comp/ppcomp.html
14. The position of Pluto
No analytical theory has ever been constructed for the planet Pluto. Our most accurate representation of the motion of this planet is from numerical integrations. Yet, a “curve fit” may be performed to these numerical integrations, and the result will be the formulae below, valid from ca 1800 to ca 2100. Compute d, our day number, as usual (section 3). Then compute these angles:
S = 50.03 + 0.033459652 * d P = 238.95 + 0.003968789 * d
Next compute the heliocentric ecliptic longitude and latitude (degrees), and distance (a.u.):
lonecl = 238.9508 + 0.00400703 * d  19.799 * sin(P) + 19.848 * cos(P) + 0.897 * sin(2*P)  4.956 * cos(2*P) + 0.610 * sin(3*P) + 1.211 * cos(3*P)  0.341 * sin(4*P)  0.190 * cos(4*P) + 0.128 * sin(5*P)  0.034 * cos(5*P)  0.038 * sin(6*P) + 0.031 * cos(6*P) + 0.020 * sin(SP)  0.010 * cos(SP) latecl = 3.9082  5.453 * sin(P)  14.975 * cos(P) + 3.527 * sin(2*P) + 1.673 * cos(2*P)  1.051 * sin(3*P) + 0.328 * cos(3*P) + 0.179 * sin(4*P)  0.292 * cos(4*P) + 0.019 * sin(5*P) + 0.100 * cos(5*P)  0.031 * sin(6*P)  0.026 * cos(6*P) + 0.011 * cos(SP) r = 40.72 + 6.68 * sin(P) + 6.90 * cos(P)  1.18 * sin(2*P)  0.03 * cos(2*P) + 0.15 * sin(3*P)  0.14 * cos(3*P)
Now we know the heliocentric distance and ecliptic longitude/latitude for Pluto. To convert to geocentric coordinates, do as for the other planets.
How Orbital Motion is Calculated
Index
Kepler’s Laws Newtonian Mechanics 
As stated earlier, the motion of a satellite (or of a planet)
However, the actual position of the satellite is given by the true anomalyφ. In polar coordinates (r,f) describing the satellite’s motion in its orbital plane, f is the polar angle. The equation of the orbit is r = a(1 – e^{2})/(1 + e cos φ) The angle φ also grows by 360^{o} each full orbit, but not at all uniformly. By Kepler’s law of areas, it grows rapidly near perigee (point closest to Earth) but slowly near apogee (most distant point). The information needed to derive φ for any time t is contained in the law of areas, but the actual calculation is not easy. The process involves an auxiliary angle, the eccentric anomaly E which like φ and M grows by 360^{o} each orbit. At perigee, all three anomalies equal zero. 
The drawing on the right gives a geometric construction of those angles (no, don’t try to puzzle out the details). The orbital ellipse is enclosed in a circle of radius a, and given a position P of the satellite, a corresponding point Q on the circle can be drawn, sharing the same line perpendicular to the ellipse’s axis. Then E is the angle between the long axis of the ellipse and the line drawn from the center of the circle to Q (“eccentric” might mean here “from the center”).
Kepler’s EquationSuppose the elements a, e and M(0) at time t=0 are given, and we need to find the value of φ at some different time t. With f known, the above equation gives r, and (r ,φ) together pinpoint the satellite’s position in its orbital plane. The first step is to derive M = M(0) + 360°(t/T) We assume the period T is known (this requires the 3rd law and is discussed for circular orbits in sections 20 and 20a). It can then be shown that the angle E satisfies “Kepler’s equation” M = E – (180°/π)e sinE where π = 3.14159256… is the ratio between the circumference of a circle and its diameter. How did that number suddenly crop up, you may ask? The fact is, the division of the circle into 360 degrees may be convenient to use (we inherited it from the ancient Babylonians) but the number 360 has no particular place in mathematics. It is probably related to the number of days in a year. The “natural” division of angles which arises in calculus and other branches of math is into radians, with 360 degrees equal to 2π = 6.2831… radians (making each radian equal to about 57.3 degrees). With angles measured in radians, Kepler’s equation simplifies to M = E – e sinE No matter which form is used, mathematics knows no formula which gives E in terms of M. However, solutions can often be approximated to any degree of accuracy by iteration–by starting with an approximate solution, then improving it again and again by an appropriate procedure (“algorithm“–more about that word, here). If the eccentricity e is not too big–the ellipse not much different from a circle–then M and E are not too different. So an initial guess E’ = M may not be too far off. Putting this guess into the term sinE gives an improved guess E” E” = M + (180°/π)e sinE’ One can now insert E” in the sinE term and get an even closer guess, and so on and so forth… until the first (say) ten decimal digits of the value of Eno longer change, at which point we may decide we have E to sufficient accuracy and stop the process. Computer handle such a process of continuous improvement (“iteration of the solution”–one form of an algorithm) very rapidly, and other methods also exist, with sufficient speed even when e is not very small. Given E, a number of formulas will give the true anomaly φ. For instance, one can first derive r = a(1 – e cos E) And then cos φ can be found from The Orbit in SpaceThe 3 remaining orbital elements are all angles giving the position of the orbit in 3 dimensions. They are described below, but their actual use belongs to a university course in orbital mechanics and will be omitted. The angles:
To orient the orbit in 3 dimensions requires a reference plane and a reference direction. For satellite orbits, the reference plane–the horizontal plane in the drawing–is usually the Earth’s equatorial plane (sometimes it is the plane of the ecliptic). The reference direction in either case is the direction from the center of the Earth to the vernal equinox (which belongs to both above planes). We will call it the x direction, since that is its role in (x,y,z) coordinates used in orbital calculations. Two nonparallel planes always intersect along a line–the way the plane of a door intersects the plane of the wall along the door’s hinge. The orbital plane and the equatorial plane (used for reference) do so too, and their intersection is called the line of nodes N. Let the origin O of our coordinates be the center of the Earth, which is also the focus of the ellipse; this point belongs to both the equatorial plane and the orbital plane, and is therefore also on their intersection line N (drawing). Then…
Suppose you have the orbital elements of some satellite, e.g. the space shuttle (you can often get them off the worldwide web). The first three (a, e, M), with M given at some particular time, enable you to calculate where the satellite will be at any time in its orbit. With (i, ω, Ω) you can then find where it would be in the sky.
