# To show perihelion precession motion of Mercury in Python with matplotlib

I'm trying to make it more general code so I can trace even a parabola or hyperbola for $$e=1, E=0$$ and $$e>1, E>0$$ respectively. And after achieving that general code I want to make the same code working for perihelion precession motion of Mercury just by adding an additional term ($$-A/r^3$$) in the energy term. Feel free to change my entire code and energy term.

from numpy import *
import matplotlib.pyplot as plt

L = 9.11*10**38     #L = angular momentum
m = 3.28*10**23     #m = mass of mercury
M = 1.99*10**30     #M = mass of sun
a = 5.8*10**10       #a = semi-major axis
G = 6.674*10**-11   #G = Gravitationl constant
k = G*M*m
E = -k/(2*a)        #E = energy
##E = 0
p = L**2/(m*k)
c = 1 + (2*E*L**2)/(m*k**2)
e = sqrt(c)       #e = eccentricity
print(e,c)
def fx(x):
r = p/(1 + e*cos(x))
return float(r)

n = 1000
phi =linspace(0,2*pi,n)
theta = zeros([n])
x = zeros([n])
y = zeros([n])

for i in range(0,n):
theta[i] = 180*phi[i]/pi

for i in range(0,n):

for i in range(0,n):

print('x =',x)
print('y =',y)
plt.plot(x,y)
plt.show()

• It looks like you are moving from a compiled language to python. Numpy arrays make life easier! Here's a better way to use them. pastebin.com/LY7f5t72 btw never use import * it can cause unexpected and hard to find problems. Right now your script crashes because $c<0$ and $e=\sqrt{c}$. Have a look at this answer for some equations for bound and unbound orbits, and this answer for a numerical integration example. Feel free to ping me with questions.
– uhoh
May 27, 2019 at 10:24
• But right now I recommend 1) you update your question to have better python (er my first link or something similar) and then 2) figure out how you are going to handle the square root of a negative number. Good luck and Welcome to Stack Exchange!
– uhoh
May 27, 2019 at 10:25