Python Kepler´s law Plotting

Question

I shall plot the earth surrounding the sun. Therefore the task is splitted into 2 subtasks. At the first task I shall approximate that the motion is one a circle.

I used following code to get the solution but somehow the programm will edit one dot instead of several. Can you pls help me with my algorithm.

So my code:

npoints = 360    
x= np.zeros((npoints,1))
y= np.zeros((npoints,1))
v_x=np.zeros((npoints,1))
v_y=np.zeros((npoints,1))
r=1
dt=1
x[0]=1.
y[0]=0.
v_x[0]=-1.
v_y[90]=-1.
v_x[180]=1.
v_y[270]=1.
for step in range(0,npoints-1):
    v_x[step+1]=v_x[step]-4*pi**2*x[step]/(r**3)*dt
    x[step+1]=x[step]+v_x[step+1]*dt
    v_y[step+1]=v_y[step]-4*pi**2*y[step]/(r**3)*dt
    y[step+1]=y[step]+v_y[step+1]*dt


plt.plot(x, y)
plt.axis([-100, 100, -100, 100])
plt.ylabel('y-axis')
plt.xlabel('x-axis')
plt.show()  

Thanks for your help :)

Answer 1

I think that there are two mistakes in your code:

1. You do not use the correct units of measure (or, you do not use them consistently.)
   As far as I can tell from the code you posted, the distance from the Sun should be measured in Astronomical Units and time in fractions of a year, so that r == 1. means a distance of 1AU (~1.49e8 km) and dt == 1. is a year. (As dt == 1. is too large, you can divide it by npoints: dt = 1./npoints. If npoints == 360, the time step will be in the order of a day.)
   Moreover, you will need to express the gravitational parameter mu, too, into a consistent unit of measure. Using the expression for the orbital period T = 2*pi*sqrt(r**3 / mu) and imposing T=1. and r=1., we obtain mu = 4 * pi**2.
2. You are imposing wrong initial conditions for the velocity.
   Let's assume (as you did) that the initial position of the Earth has coordinates (x=1 AU, y=0 AU). The speed is tangent to the orbit, so in the chosen reference frame it has only a vertical component (v_y) and its module is given by the equation for the Circular Velocity. So you impose (v_x=0 AU/yr, v_y=sqrt(mu/r) AU/yr). Note that if you impose this set of initial conditions, you do not have to impose any other condition because the problem is already well posed. (Moreover, conditions like v_y[90]=-1. are overwritten in the for-loop and do not affect your computations at all.)

Here is the complete code:

import numpy as np              # please next time include the relevant
import matplotlib.pyplot as plt # `import` statements and variable
pi = np.pi                      # definitions

npoints = 360
r = 1.          # AU
dt = 1./npoints # fractions of a year
mu = 4 * pi**2  # 
x = np.zeros(npoints)
y = np.zeros(npoints)
v_x = np.zeros(npoints)
v_y = np.zeros(npoints)

# Initial Conditions
x[0] = r               # (x0 = r, y0 = 0) AU
v_y[0] = np.sqrt(mu/r) # (v_x0 = 0, v_y0 = sqrt(mu/r)) AU/yr

for step in range(0,npoints-1):
    v_x[step+1]=v_x[step]-4*pi**2*x[step]/(r**3)*dt
    x[step+1]=x[step]+v_x[step+1]*dt
    v_y[step+1]=v_y[step]-4*pi**2*y[step]/(r**3)*dt
    y[step+1]=y[step]+v_y[step+1]*dt

plt.plot(x, y, 'bo')
plt.axis('equal')
plt.show()