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

What should `x`, `y`, `v_x` and `v_y` be? Are they coordinates of points in the orbit and the components of the speed of the Earth? — Pier Paolo, Feb 09 '16 at 22:38
Basically the x and y are the coordinates of the possible positions of the earth. The position of the sun is at zero/zero. v_x and v_y represents the linear change of the position in x and y direction. — Heisenberg, Feb 09 '16 at 22:44

Pier Paolo · Answer 1 · 2016-02-10T02:32:58.923

I think that there are two mistakes in your code:

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

the difference is: you find a parametrization which solves the problem: the task is to find a solution of the physical problem: i just deduce what v_x and v_y shall do: in x direction we have: F_x= G * Me*Ms*x F_y=G*Me*Ms*y Me mass of earth, Ms mass of sun. We dont respect that the sun also drains the moon. With Newtons equation (F=a*m) we get: 1.) dv_x/dt = - G*Ms*x/r³ and 2.) dv_y/dt = - G*Ms*y/r³ if we taylor now the for small t: v_x/y(t+dt)=v_x/y(t)-4*(pi²/r_i³)*x_i/y_i*dt — Heisenberg, Feb 10 '16 at 00:38

Python Kepler´s law Plotting

1 Answers1