# 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}}\left(1+e\cos\phi+3\frac{m^{2}}{B^{2}}\left(1+e\phi \sin\phi +e^{2}\left(\frac{1}{2}-\frac{1}{6}\cos2\phi\right)\right)\right),$$ 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:

\begin{align*} r_{n+1} & = 2r_{n}-r_{n-1}+\Delta \tau^{2}\left(-\frac{GM}{{r_{n}}^{2}}+\frac{B^{2}}{{r_{n}}^{3}}-\frac{3GMB^{2}}{{r_{n}}^{4}}\right) \\ \phi_{n+1} & = \phi_{n}+\Delta\tau\frac{B}{\left[\frac{1}{2}(r_{n+1}+r_{n})\right]^{2}}, \end{align*}

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

• External links go stale and disappear as external sites go offline; please consider inlining your Python code stored behind those links into your question so that it is preserved. Apr 18 at 1:40

The Newtonian acceleration vector is: $$$$\vec{g}_n = \frac{\hat{d} G M}{{\lvert\lvert \vec{d} \rvert\rvert}^2}.$$$$
One important value is closely related to the kinematic time dilation: $$$$\label{eq_kinematic} \alpha = 2 - \sqrt{1 - \frac{\lvert\lvert \vec{v}_{o}\rvert\rvert^2}{c^2}}.$$$$ Another important value is the gravitational time dilation: $$$$\beta = \sqrt{1 - \frac{R_{s}}{\lvert \lvert \vec{d} \rvert \rvert}}.$$$$
Finally, the semi-implicit Euler integration is: \begin{align} \vec{v}_{o}(t + \delta_t) &= \vec{v}_{o}(t) + \delta_{t} \alpha \vec{g}_n, \\ \ell_{o}(t + \delta_t) &= \ell_{o}(t) + \delta_{t} \beta \vec{v}_{o}(t + \delta_t). \end{align}