<p>I'm trying to make it more general code so I can trace even a parabola or hyperbola for <span class="math-container">$e=1, E=0$</span> and <span class="math-container">$e&gt;1, E&gt;0$</span> 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 (<span class="math-container">$-A/r^3$</span>) in the energy term. <em>Feel free to change my entire code and energy term.</em></p>

<pre><code>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)
radius = zeros([n])
theta = zeros([n])
x = zeros([n])
y = zeros([n])

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

for i in range(0,n):
    x[i] = radius[i]*cos(phi[i])

for i in range(0,n):
    y[i] = radius[i]*sin(phi[i])

print('r =',radius)
print('x =',x)
print('y =',y)
plt.plot(x,y)
plt.show()
</code></pre>

CC BY-SA 4.0
2019-05-27T10:24:38.800Z
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. https://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](https://space.stackexchange.com/a/27604/12102) for some equations for bound and unbound orbits, and [this answer](https://space.stackexchange.com/a/23257/12102) for a numerical integration example. Feel free to ping me with questions.
uhoh
CC BY-SA 4.0
2019-05-27T10:25:48.307Z
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
CC BY-SA 4.0
2019-05-27T10:32:11.927Z
also see [this answer](https://space.stackexchange.com/a/22170/12102)
uhoh
<orbit><planet><astrophysics><python>
To show perihelion precession motion of Mercury in Python with matplotlib
zircon
# To show perihelion precession motion of Mercury in Python with matplotlib

(Asked by: username_0 on 2019-05-27)

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)
radius = zeros([n])
theta = zeros([n])
x = zeros([n])
y = zeros([n])

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

for i in range(0,n):
    x[i] = radius[i]*cos(phi[i])

for i in range(0,n):
    y[i] = radius[i]*sin(phi[i])

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

**Tags:** orbit, planet, astrophysics, python

> **Comments:**
> **username_1**: 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. https://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](https://space.stackexchange.com/a/27604/12102) for some equations for bound and unbound orbits, and [this answer](https://space.stackexchange.com/a/23257/12102) for a numerical integration example. Feel free to ping me with questions.  
> **username_1**: 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!*  
> **username_1**: also see [this answer](https://space.stackexchange.com/a/22170/12102)  

---

## 0 answers

