Plotting an orbit in Python, but when x<0, trajectory suddenly goes linear

Question

I'm plotting in a simple 2D plane the path taken by a body in motion past another gravitationally attractive body.

Every loop, the time is incremented by a second and the body's new position is calculated and printed out. I then paste the results into a spreadsheet and graph it.

Things look ok until the body's x component becomes negative - and then the trajectory goes linear, and scoots off to the top left.

This is all for my own solo entertainment, I'm no student. So after scratching my head for a bit I've finally lumped for asking someone for help. I've probably missed something obvious. I suspect my trigonometry is lacking something.

I'm using Python 2.7.10

import sys
import os
import math

mu = 4.0*(10**14)
massAst = 1
earthRadius = 6371000.
alt = 100000.
r = earthRadius+ alt
rTheta = 270.
rAngtoX = math.radians(rTheta)
tInc = 1 ## increment time by 1 seconds - one update of pos&velocity per second of travel
calcTime = 1100 ## simulation runtime (86400 seconds = 1 day) 10 mins
t = 1 ## integral of time t to use in the calcs in the loop.
printbell = 120 ## print results now
printclock = 0
hourCount = 0

## Initialise velocity vectors for Asteroid:
uAstX = 1500.
uAstY = 0.
vAstX = 0.
vAstY = 0.
## Displacement
dAstX = r*math.cos(rAngtoX)
dAstY = r*math.sin(rAngtoX)

for i in range(0, calcTime):
    acc = -1*(mu/r**2)

    accX = acc*math.cos(rAngtoX)
    accY = acc*math.sin(rAngtoX)

    vAstX = uAstX + accX*t ## new value for velocity in X direction
    vAstY = uAstY + accY*t ## and in Y

    deltaDAstX = uAstX*t + 0.5*accX*(t**2) ## change in position over this time interval
    deltaDAstY = uAstY*t + 0.5*accY*(t**2)

    dAstX = dAstX + deltaDAstX
    dAstY = dAstY + deltaDAstY

    uAstX = vAstX
    uAstY = vAstY
    ## Now calculate new angle and range
    ## tan(theta) = dAstY/dAstX, so:
    rAngtoX = math.atan(dAstY/dAstX) ##+(2*3.141592654)
    ##print 'theta:', math.degrees(rAngtoX)
    r = dAstY/math.sin(rAngtoX)
    ## if i == print
    print dAstX, '  ', dAstY

To format your code correctly, copy-paste it into the question box, then select everything you pasted (including the parts that look okay in the preview) and hit Ctrl-K or the button with the braces. This is very important for Python code, as we need to see the indentation exactly as it is in the code you actually ran. — user2357112 supports Monica, Oct 03 '18 at 19:29

score 1 · Answer 1 · answered Oct 03 '18 at 20:32

As dAstX approaches zero, dAstY/dAstX will approach a division-by-zero... Which will cause all sorts of problems (roundoff issues at the very least).

I'd recommend keeping the x/y components of distance/velocity/acceleration separate. The distance between the objects is important, of course, but that can be calculated using r=sqrt(dAstX**2 + dAstY**2).

score 0 · Answer 2 · answered Oct 03 '18 at 20:34

I'm not familiar with the math, but I took your code, modified it to run on my machine, plottet the data using the seaborn library and I came up with this:

This is the code I used:

import math
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt


def calc(calcTime):
    mu = 4.0*(10**14)
    earthRadius = 6371000.0
    alt = 100000.0
    r = earthRadius+ alt
    rTheta = 270.0
    rAngtoX = math.radians(rTheta)
    t = 1               # integral of time t to use in the calcs in the loop.

    # Initialise velocity vectors for Asteroid:
    uAstX = 1500.0
    uAstY = 0.0
    # Displacement
    dAstX = r*math.cos(rAngtoX)
    dAstY = r*math.sin(rAngtoX)

    point_list = []
    for i in range(0, calcTime):
        acc = -1*(mu/r**2)

        accX = acc*math.cos(rAngtoX)
        accY = acc*math.sin(rAngtoX)

        vAstX = uAstX + accX*t                  # new value for velocity in X direction
        vAstY = uAstY + accY*t                  # and in Y

        deltaDAstX = uAstX*t + 0.5*accX*(t**2)      # change in pos over time interval
        deltaDAstY = uAstY*t + 0.5*accY*(t**2)

        dAstX = dAstX + deltaDAstX
        dAstY = dAstY + deltaDAstY

        uAstX = vAstX
        uAstY = vAstY
        # Now calculate new angle and range
        # tan(theta) = dAstY/dAstX, so:
        rAngtoX = math.atan(dAstY/dAstX)            #+(2*3.141592654)
        # print 'theta:', math.degrees(rAngtoX)
        r = dAstY/math.sin(rAngtoX)
        # if i == print
        if i % 5 == 0:
            print('{:05d} | {:15.2f} | {:15.2f}'.format(i, dAstX, dAstY))
            point_list.append((i, dAstX, dAstY))

    df = pd.DataFrame(data=point_list, columns=['i', 'x', 'y'])
    return df


if __name__ == '__main__':
    df = calc(950)
    sns.scatterplot(data=df, x='x', y='y')
    plt.show()

My analysis: the spaces between dots gets bigger each time (note: I only plottet every 5th point to make the image more readabl,e in my opinion). From a physics perspective, that would indicate that the object is gaining speed.

Isn't it possible that your calculations are correct and that the object is leaving the orbit because it gained enough speed (aka energy) to leave the gravitational field of the center mass (aka the earth) ?

As I said, I'm not familiar with the specific math, but to me it makes sense that the object could break free from the orbit with the speed it is gaining in the half turn.

Plotting an orbit in Python, but when x<0, trajectory suddenly goes linear

2 Answers2