# Precession of Mercury Python simulation

I was trying to simulate the precession of Mercury based on the perturbed solution: $$\frac{1}{r}=\frac{m}{B^{2}}(1+e\cos\phi+3\frac{m^{2}}{B^{2}}(1+e\phi \sin\phi +e^{2}(\frac{1}{2}-\frac{1}{6}\cos2\phi))),$$ where $e$ is the eccentricity of mercury, B is a constant which is the angular momentum per mass. here is my code in Python: http://sprunge.us/TFhU

In the code I have to multiply a large number(SF) so that the precession is slightly visible(or you can make a bigger SF value then the precession would be much easier to see but it would distort the orbits). I guess it is okay if I use this as a illustration for the mercury precession, however, I want to be able to obtain the value of the precession rate which is known as 43 arcseconds per century, does anyone know how to do that in the code or anyone has python code that does a better job? (I also came through a model based on these two equations: $$r_{n+1} = 2r_{n}-r_{n-1}+\Delta \tau^{2}(-\frac{GM}{{r_{n}}^{2}}+\frac{B^{2}}{{r_{n}}^{3}}-\frac{3GMB^{2}}{{r_{n}}^{4}}),\phi_{n+1}=\phi_{n}+\Delta\tau\frac{B}{[\frac{1}{2}(r_{n+1}+r_{n})]^{2}},$$ where B is same as above, but I couldn't get it to work.)

(Update: I used Vpython and Leapfrog method to do it and the simulation worked out all right. here is the code and modify it as however you want to, http://sprunge.us/ehgS (in linux the sphere doesn't leave trail probably because the vpython in linux is very out of date). However, I still need to find a way to output the precession rate 43" out of the data.)