Python Plotting Help: Solving Circular Orbit with Matplotlib and Numpy

[yik-boh]

Homework Statement

We need to plot a two points (representation of planets). One is stationary and one is supposed to move on a circular orbit. We are using Python, Matplotlib and Numpy. We're not getting a circular orbit. We're just getting a line.

Homework Equations

a = F / mass1
x = x + (v * t) + (0.5 * a * t**2)
y = y + (v * t) + (0.5 * a * t**2)

The Attempt at a Solution

import matplotlib.pyplot as plt
from math import *

i = 0
t = 0
x1 = []
y1 = []
x = 5
y = 0
bx = 2
by = 0
v = 0
mass1 = input("input mass 1: ")
mass2 = input("input mass 2: ")
G = 6.67e-11
bigx = x
bigy = y
smallx = bx
smally = by

while i < 10000:
    i = i+1
    while t<100:
        r = sqrt((x-bx)**2 + (y-by)**2)
        t = t + 0.001
        F = G * mass1 * mass2 / r**2
        a = F / mass1
        v = v + a * t
        x = x + (v * t) + (0.5 * a * t**2)
        y = y + (v * t) + (0.5 * a * t**2)
        y1.append(y)
        x1.append(x)
       
        
    
plt.plot(x1,y1)
plt.plot(smallx,smally, marker ='o',linestyle = '-', color = 'b')
plt.plot(bigx,bigy, marker ='o', color = 'r')
plt.axis([-20,20,-20,20])
plt.show()

 

[Joffan]

I see no determination of what the x and y components of F, a or v are. Also, once you do that, initial relative velocity of zero will indeed give a straight-line attraction between the two masses, known as "dropping" , as opposed to the orbit I think you want.

 

[yik-boh]

What do you mean no determination? I'm sorry. I'm really confused. I mean how can I set the x and y components of F, a, or v?

 

[Joffan]

The components of F depend on the relative position of the two bodies. The force acts along the line between the two, so there should be an x component and a y component. That translates directly into the x and y components of the acceleration, and that adds appropriately to the x and y components of the velocity, which in turn updates the x and y coordinates of the positions.

You need to determine those components, not simply add one acceleration to everything.

Suppose we start with one big mass at (0,0) and arbitrarily say it's not going to move. Then put the moving mass at say (0,100) with a velocity of (1,0). NOTE that the initial direction of motion is important. Then if the initial acceleration is 0.1, say, it acts on that moving mass in the direction (0,-0.1). So the velocity at the next step is (1,-0.1) and the new position is (1,99.9) and the next force calculation takes place on the new vector.

 

[yik-boh]

But what are the directions of the forces here? Towards the other planet and the other one is where? How about my radius, r = sqrt((x-bx)**2 + (y-by)**2), is this the one I should use?

 

[Joffan]

The force on the moving planet acts along the line joining the two bodies.

Your distance is fine. The problem is that you are treating the force, acceleration and velocity as scalars when they are vectors. If you calculate the offset-x and offset-y between the two planet positions, that gives you an offset vector and the force and acceleration act along the same direction:

        offset_x = (x-bx)
        offset_y = (y-by)
        r = sqrt(offset_x**2 + offset_y**2)
        F = G * mass1 * mass2 / r**2
        F_x = (-offset_x/r)*F
        F_y = (-offset_y/r)*F

You'll also need to have the idea of delta-t when updating the velocity vector using the acceleration vector, and the position using the velocity vector.

 

[yik-boh]

import matplotlib.pyplot as plt
from math import *

i = 0
t = 0
x1 = []
y1 = []
x = 5
y = 0
bx = 2
by = 0
mass1 = input("input mass 1: ")
mass2 = input("input mass 2: ")
G = 6.67e-11
bigx = x
bigy = y
smallx = bx
smally = by

while i < 10000:
    i = i+1
    while t < 100:
        offset_x = (x - bx)
        offset_y = (y - by)
        r = sqrt(offset_x**2 + offset_y**2)
        t = t + 0.001
        F = G * mass1 * mass2 / r**2
        F_x = (-offset_x/r) * F
        F_y = (-offset_y/r) * F
        a_x = F_x / mass1
        a_y = F_y / mass1
        v = 0
        v_x = v + a_x * t
        v_y = v + a_y * t
        x = x + (v_x * t) + (0.5 * a_x * t**2)
        y = y + (v_y * t) + (0.5 * a_y * t**2)
        x1.append(x)
        y1.append(y)
       
        
plt.figure()    
plt.plot(x1,y1)
plt.plot(smallx,smally, marker ='o',linestyle = '-', color = 'b')
plt.plot(bigx,bigy, marker ='o', color = 'r')
plt.axis([-20,20,-20,20])
plt.show()

That's my current code. I tried to to the components thingy but it still doesn't work. :( I only get a horizontal line.

By the way, here's the guideline that my teacher sent us before:1) Use the kinematical equation x_i = x_{i-1} + v_{i-1} t + 1/2 a t^2. Do the same for the y-coordinate. The index i = 0 denotes your initial positions and velocities, you are free to choose that.
2) For every time step (make it small like t += 0.001 or t += 0.0001), you should replace the x_{i-1} and v_{i-1}, which are the previous position and velocity of the particle at the previous time step.
3) For your acceleration a, use the formula a = F / m, where F is the gravitational force between the two masses.
4) Feed all the x_i's and the y_i's into your x = [ ] and y = [ ] arrays (declare these first as empty lists). You will then have a list of values for each of x and y. Use x.append(x_i) to add the x_ith value to the list. After all the iterations, the two lists should have the same number of elements. Mine has 10,000 elements each.
5) Plot using the following codes (be sure you have "import matplotlib.pyplot as plt" in your code header):

plt.figure()
plt.plot(x,y)
plt.savefig("name of the image file.png")

The .png file of the graph will be located in the same directory as your python file. You can also set the axes limits using plt.ylim(-20,20) and plt.xlim(-20,20). This will output a graph whose x- and y- axes are ranging from -20 to 20.

 

[Joffan]

Your v_x should be an update on an existing v_x (not v as you have it), similarly v_y, and of course they should have initial values, preferably to represent movement at right angles to the initial line between the planets as per my earlier comment.

The update velocity/position equations should be using delta_t rather than t, which in your code seems to be 0.001 (and you should set delta_t = 0.001 at the beginning and use that instead of the constant).

I don't know what your outer loop is for. It doesn't seem to do anything.

What masses ae you using as input? Note that G is small, so unless the mass of the fixed planet (mass2) is large, you may not see enough gravitational force to make the moving planet orbit.


Python note: It may be also that your teacher expects you to use tuples rather than separate coordinate variables.

 

[yik-boh]

import matplotlib.pyplot as plt
from math import *

i = 0
t = 0
x1 = []
y1 = []
x = 5
y = 0
bx = 2
by = 0
mass1 = input("input mass 1: ")
mass2 = input("input mass 2: ")
G = 6.67e-11
bigx = x
bigy = y
smallx = bx
smally = by

while i < 10000:
    i = i+1
    while t < 100:
        offset_x = (x - bx)
        offset_y = (y - by)
        r = sqrt(offset_x**2 + offset_y**2)
        dt = t + 0.001
        F = G * mass1 * mass2 / r**2
        F_x = (-offset_x/r) * F
        F_y = (-offset_y/r) * F
        a_x = F_x / mass1
        a_y = F_y / mass1
        v = 0
        v_x = v + a_x * dt
        v_y = v + a_y * dt
        x = x + (v_x * dt) + (0.5 * a_x * dt**2)
        y = y + (v_y * dt) + (0.5 * a_y * dt**2)
        x1.append(x)
        y1.append(y)
       
        
plt.figure()    
plt.plot(x1,y1)
plt.plot(smallx,smally, marker ='o',linestyle = '-', color = 'b')
plt.plot(bigx,bigy, marker ='o', color = 'r')
plt.axis([-20,20,-20,20])
plt.savefig("project.png")
plt.show()

Now here's my code. I changed t to dt and now my program takes a long time for it to finish and sometimes I gets a memory error. Also did I make the right change with the delta_t equation you're telling me? Can you also please give me an example about the v_x and v_y equation. I really don't get it. :( Sorry! Also, can you give me an example of tuple using my code? I just want to see if it can clean my code because it's a little bit dirty. Thanks for all the help you're giving me! :)
Here's the code my friend made for her group:

import numpy as np
import matplotlib.pyplot as plt 

fig = plt.figure()
fig.patch.set_facecolor('blue')
fig.patch.set_alpha(0.7)

plt.title('Circular Orbit')

G = 6.67*10**(-11) #The Gravitational Constant
M = 1.989*10**(30) #Mass of Sun
m = 1.9*10**(27) #Mass of Jupiter
r = 1.60*10**(22) #radius#for circular orbit 
def j(t):
    return np.cos(t)*h**2/(u*(1+e*np.cos(t)))
def i(t):
    return np.sin(t)*h**2/(u*(1+e*np.cos(t)))

v = 0.091058463
h = r*v
u = G*(M+m)
e = (r*v**2)/(G*M)-1

n = np.arange(0.0, 20.0, 0.001)

plt.figure(1)
plt.plot(j(n), i(n), 'b-', 0, 0, 'y*', 1.58995e+22, 1.11789e+21, 'ro')
plt.show()

I'm not sure if it's correct. It gives you a circular orbit but it doesn't even compute for the velocity and other stuff so I did not copy it and it's way toooo different from my code.

 

[Joffan]

I'm surprised it ever finished. You still need to increment t, like this:
t = t + dt

Otherwise there are a lot of points you haven't addressed from my previous post.

I get the impression you are just throwing bits of code in; you need to take a deep breath, step back and try to understand what is going on.

And yes, I think your friend is not using physics, just drawing circles.

 

[yik-boh]

So I'm right to use this one?

dt = t + 0.001

but I should also put the code t = t + dt?

Also I don't get the idea you're trying to tell about my v_x and v_y. :(

 

[Joffan]

dt = t + 0.001 is nonsense.

dt should be initialized to 0.001, or whatever time step you want to use, at the top. Then in the loop it should be used to increment t, and to calculate the updates to position and velocity (which you already have done).

v_x and v_y should also be initialized at the top, and updated using a_x and a_y respectively. You still have v (on its own), which is not needed.

Given the formulas you use, the position should be updated before the velocity.

 

[Joffan]

As a test set, putting the fixed mass at (0,0) with mass2=1e12, the moving mass at (0,10) with an initial velocity of (2.6,0) and mass1=1000, you should get a fairly circular orbit with a period of about 25.

Increasing the initial velocity to say (3,0) will give you a distinct ellipse. Increasing it to (4,0) will not result in an orbit; the moving mass escapes.

 

[yik-boh]

Super thank you! It's now working. I finally understood the changes you want me to make!

Here's my working code:

import matplotlib.pyplot as plt
from math import *

note = " Note: Use 1000 as mass1 and 1e12 as mass2 to get the proper results. "

print note

i = 0
t = 0
x1 = []
y1 = []
x2 = []
y2 = []
x = 0 #x-coord of the point for the circular orbit
y = 10 #x-coord of the point for the circular orbit
n = 0 #x-coord of the point for the ellipse orbit
m = 10 #y-coord of the point for the ellipse orbit
bx = 0
by = 0
dt = 0.001
v_x = 2.6 #x component of the velocity for the circular orbit
v_y = 0 #y component of the velocity for the circular orbit
w_x = 3 #x component of the velocity for the ellipse orbit
w_y = 0 #y component of the velocity for the ellipse orbit
mass1 = input("input mass 1: ")
mass2 = input("input mass 2: ")
G = 6.67e-11
bigx = x
bigy = y
smallx = bx
smally = by

while i < 10000:
    i = i+1
    while t < 100:
        offset_x = (x - bx) #for the circular orbit
        offset_y = (y - by) #for the circular orbit
        offset_n = (n - bx) #for the ellipse orbit
        offset_m = (m - by) #for the ellipse orbit
        r = sqrt(offset_x**2 + offset_y**2) #for the circular orbit
        s = sqrt(offset_n**2 + offset_m**2) #for the ellipse orbit
        t = t + dt
        F = G * mass1 * mass2 / r**2 #for the circular orbit
        H = G * mass1 * mass2 / s**2 #for the ellipse orbit
        F_x = (-offset_x/r) * F #for the circular orbit
        F_y = (-offset_y/r) * F #for the circular orbit
        H_x = (-offset_n/s) * H #for the ellipse orbit
        H_y = (-offset_m/s) * H #for the ellipse orbit
        a_x = F_x / mass1 #for the circular orbit
        a_y = F_y / mass1 #for the circular orbit
        l_x = H_x / mass1 #for the ellipse orbit
        l_y = H_y / mass1 #for the ellipse orbit
        v_x = v_x + a_x * dt #for the circular orbit
        v_y = v_y + a_y * dt #for the circular orbit
        w_x = w_x + l_x * dt #for the ellipse orbit
        w_y = w_y + l_y * dt #for the ellipse orbit
        x = x + (v_x * dt) + (0.5 * a_x * dt**2) #for the circular orbit
        y = y + (v_y * dt) + (0.5 * a_y * dt**2) #for the circular orbit
        n = n + (w_x * dt) + (0.5 * l_x * dt**2) #for the ellipse orbit
        m = m + (w_y * dt) + (0.5 * l_y * dt**2) #for the ellipse orbit
        x1.append(x) #for the circular orbit
        y1.append(y) #for the circular orbit
        x2.append(n) #for the ellipse orbit
        y2.append(m) #for the ellipse orbit
       
        
plt.figure()    
plt.plot(x1,y1)
plt.plot(x2,y2)
plt.plot(smallx,smally, marker ='o',linestyle = '-', color = 'b')
plt.plot(bigx,bigy, marker ='o', color = 'r')
plt.axis([-20,20,-20,20])
plt.savefig("project.png")
plt.show()

Thanks a lot! :D

 

[Joffan]

Congratulations... even though you never deleted the useless outer loop. :-)

Of course changing the mass of the fixed body will also change the characteristics of the orbits, so properly speaking you don't know in advance of the inputs whether the two orbits you try out are circular, elliptical or non-orbits.