1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Python Help in Plotting

  1. Mar 21, 2012 #1
    1. The problem statement, all variables and given/known data

    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.


    2. Relevant equations

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

    3. The attempt at a solution

    Code (Text):
    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()


     
     
  2. jcsd
  3. Mar 21, 2012 #2
    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" :wink: , as opposed to the orbit I think you want.
     
  4. Mar 21, 2012 #3
    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?
     
  5. Mar 21, 2012 #4
    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.
     
  6. Mar 22, 2012 #5
    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?
     
  7. Mar 22, 2012 #6
    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:

    Code (Text):

            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.
     
  8. Mar 22, 2012 #7
    Code (Text):



    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.
     
  9. Mar 22, 2012 #8
    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.
     
    Last edited: Mar 22, 2012
  10. Mar 22, 2012 #9
    Code (Text):


    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:

    Code (Text):


    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.
     
  11. Mar 22, 2012 #10
    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.
     
  12. Mar 23, 2012 #11
    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. :(
     
  13. Mar 23, 2012 #12
    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.
     
  14. Mar 23, 2012 #13
    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.
     
  15. Mar 24, 2012 #14
    Super thank you!! It's now working. I finally understood the changes you want me to make!

    Here's my working code:

    Code (Text):



    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
     
  16. Mar 24, 2012 #15
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Python Help in Plotting
Loading...