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
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)
  • 1
    $\begingroup$ 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. $\endgroup$ – uhoh May 27 '19 at 10:24
  • $\begingroup$ 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! $\endgroup$ – uhoh May 27 '19 at 10:25
  • $\begingroup$ also see this answer $\endgroup$ – uhoh May 27 '19 at 10:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.