2D orbital Physics model showing strange behaviour

Question

I'm trying to create a 2D environment which simulates basic orbital physics using pygame. However when I run it, it exhibits some weird behaviour (edit: and some weirder behaviour)such as the satellite bouncing off the x and y axes, and the satellite moving in a wave like pattern along the axis toward the sun. Below is the code that calculates the acceleration of a body caused by the central sun. BM is body mass, BP is Body position [x,y], SM is sun mass and SP is sun position [x,y].

def CalcGrav(BM,BP,SM,SP):  
    Dist = SubList(BP,SP)  
    Mass = BM*SM/100  
    for i in range(0,2):  
        if(Dist[i] == 0):  
            Dist[i] += 0.01  
    AV = [-(Mass/Dist[0]),-(Mass/Dist[1])]  
return AV  

The function is then used in the following context:

TraceAcc = CalcGrav(BallMass,TracePosition,SunMass,SunPosition)
TraceVector = SumList(TraceVector, TraceAcc)
TracePosition = SumList(TracePosition,TraceVector)

SumList(A,B) simply adds the items of a list together: (1,3) + (2,6) = (3,9)
The current velocity vector and acceleration vector are stored as lists: [x,y]

Can anyone explain the strange behavior or show me what I've done wrong?
Thanks.

Answer 1

The problem is that you are not constructing a realistic physics model of orbital motion.

What you need to do for a somewhat realistic simulation, is construct a simplified Kepler orbit model. You also will need to work in polar coordinates (radius, angle) and convert to x, y coordinates.

What you will need to do is construct a Kepler orbit equation for an elliptical orbit in code.

The Kepler equation looks like (see the Wikipedia article on orbit modeling for more information)

This describes the radius of the orbit r in terms of the angle of the orbit v (actually Greek nu), where the angle is zero at "periapsis" which means the point of closest approach of orbiting body and central body.

You will need to define your own eccentricity (e) as some value between 0 and 1 - I suggest something close to the Earth's orbital eccentricity, approximately 0.02. That will give you a nearly circular orbit. If you want a more elliptical one, increase the eccentricity, but not more than 1.0.

Similarly a is the semimajor axis, i.e half the longest axis of the ellipse. Either define a value that makes sense for your game coordinates or if you want to compute it using the body masses, you will need to look up some references on Kepler orbits, start with the Wikipedia page above.

Even the Kepler model is vastly simplified from real orbital mechanics, but it will at least produce a semi realistic game simulation.

Answer 2

Turns out it was much simpler than that. I was dividing by individual x,y distances, when what I needed to do was work out the vector towards the sun, calculate the straight line distance between the two objects, then divide both parts of the vector by that distance squared:

def CalcGrav(MM,MP,SM,SP):
    Vect = SubList(SP,MP)
    D = math.sqrt(Vect[0]**2 + Vect[1]**2)**2
    Vect = TimesList(Vect,[1/D,1/D])
    return Vect

Answer 3

Try a smaller step (smaller distance increment). 0.002 rather than 0.01 for instance. I've had a similar problem in my own orbital physics lab and this worked.