2D Orbits making flower petal shapes

[Wheffle]

So I decided on a whim to see if I could program a simple 2D model of a satellite orbiting a large body. It's very simple code that calculates the force of gravity every step on the ship and applies the force to its current velocity vector. I altered the gravity constant so that calculations wouldn't involve giant numbers and I'm also ignoring the mass of the satellite, I figured "G*m1*m2" is a constant anyway. For a relatively circular orbit, it came out looking nice:

There are clearly some irregularities in the orbit. So I pushed the "ship" into an elliptical orbit, and I got this pattern:

So what's causing the orbit to "flower petal" like that? Does that happen in nature? Is it a result of modeling orbits on a two dimensions vs three? Are the shortcuts from my calculations causing problems? I suppose the body itself should technically be wobbling a bit. How does KSP get around this kind of issue?

[Bill Phil]

I think you don't have all of the Keplerian elements properly defined.

[pincushionman]

What exactly are you adding to the velocity vector? What equations are you integrating, and what kind of integrator did you write?

[Camacha]

It looks like it is precessing, which happens in real life too. Check out the precession of Mercury:

[Bill Phil]

It looks like it is precessing, which happens in real life too. Check out the precession of Mercury:

http://physics.ucr.edu/~wudka/Physics7/Notes_www/img249.gif

That's what I thought, too.

I think that the reference frame isn't accounting for the precession.

[K^2]

It'd be nice to see your code to actually compute the forces and update the position of the satellite. I have a few guesses on what might be causing this, but I can't be sure on which one it is without seeing the code.

[Findthepin1]

Is it simulating the satellite's effect on the planet?

[YNM]

Probably something wrong in the acceleration calculation ? I don't know for sure, but it could also be due to computing steps... KSP slightly suffers from that you know.

[Robotengineer]

It looks like it is precessing, which happens in real life too. Check out the precession of Mercury:

http://physics.ucr.edu/~wudka/Physics7/Notes_www/img249.gif

IIRC, precession wasn't something that can be produced in newtonian physics, the precession of Mercury is one of the reasons Relativity was needed.

[Alias72]

So what's causing the orbit to "flower petal" like that? Does that happen in nature? Is it a result of modeling orbits on a two dimensions vs three? Are the shortcuts from my calculations causing problems? I suppose the body itself should technically be wobbling a bit. How does KSP get around this kind of issue?

1. Can't tell without code

2. Orbits process in nature but Newtonian physics wont simulate that

3. 2 vs 3 dimension does not matter as all orbits occur on an orbital plane (2d surface)

4. maybe. Must see code

5. yes there should be planetary wobble but for immensely disproportionate masses it should be negligible.

6. KSP uses conic sections during timewarp (the orbit is "locked" in advance) also It probably is more careful with its numbers.

without seeing the code the only thing I can suggest is to look for errors in numerical simulation. orbital mechanics is prone to cascading error.

[problemecium]

Make sure you're dividing by r² rather than r, and that your timestep isn't too big.. Otherwise, I'm just going to go with the precession hypothesis.

[Jovus]

IIRC, precession wasn't something that can be produced in newtonian physics, the precession of Mercury is one of the reasons Relativity was needed.

You do not recall correctly. Precession is a perfectly Newtonian phenomenon, and frankly usually to be expected in all but the simplest situations. (Two bodies attracted to one another by an inverse square law is one of those simplest situations.)

As to whether or not the OP's 'problem' is merely that the code correctly predicts a precession he didn't expect, or if the program's calculations are incorrect...well, I can't answer that without seeing the code.

[Wheffle]

So it was quite silly not to put any kind of code up. I apologize for that. I'm actually using Game Maker: Studio, but here's the pseudo-code:

//script for getting force magnitude of gravity for a ship/planet pair
public double get_gravity_force(planet, ship)
{
    var G = 0.001; //constant I chose

    var x1 = planet.x;
    var y1 = planet.y;
    var x2 = ship.x;
    var y2 = ship.y;

    var r2 = (point_distance(x1, y1, x2, y2))^2;

    if (r2 == 0) return 0;

    //ignoring ship mass (G*m1*m2), this might be causing issues
    return (G*planet.mass/r2);
}

//script that runs every step for the ship
//ship has local variables "x", "y", "v" and "direction"
{
   //get gravity force magnitude and direction
   var fg = get_gravity_force(planet, this);
   var fg_dir = point_direction(x, y, planet.x, planet.y);

   //find x offset caused by current velocity and gravity
   var xx = v*(cos(direction)) + fg*cos(fg_dir));
   var yy = v*(sin(direction)) + fg*sin(fg_dir));

   //change direction and velocity for next step
   direction = point_direction(x, y, x+xx, y+yy);
   v = point_distance(x, y, x+xx, y+yy);

   //update position
   x = x + xx;
   y = y + yy;
}

"point_direction" and "point_distance" are provided scripts from the IDE, and the IDE actually provides scripts not used in this pseudo-code for the trigonometry, which I'm using.

[K^2]

This is terribly inefficient, because point_direction computes atan2, and then you run and compute sin and cos from that. You should work out how to avoid this. Hint: What are the definitions of sin and cos and can you compute them from your data?

Otherwise, everything looks good. Precession you are seeing is almost certainly the result of time increment being too large and your use of Euler Method for integration. Try reducing time step by re-writing your position update as xx = vx*dt + 0.5*ax*dt*dt. Then try different values for dt starting with 1 and going down. By the way, your missing factor of 0.5 is also wrong, but it just works as arbitrary increase of planet's gravity, which won't affect the types of orbit you get.

To get even better results, you can try using better integration method. Nothing works great for gravity, but almost anything will be better than Euler Method. Velocity Verlet is always my first recommendation, because it's really easy to implement and it gives great improvement.

[Yourself]

    var r2 = (point_distance(x1, y1, x2, y2))^2;

In GML, ^ is the bitwise xor operator, not exponentiation. This essentially results in your force being inversely proportional to the distance (for distances significantly larger than 2), which results in the characteristic petal pattern of your orbits. I used Game Maker in the early 2000s and people were making this mistake back then, too. You could use the power method instead or just compute the distance first and then multiply it by itself. You could also manually compute distance by just summing squares because it avoids an intermediate square root, but worrying about optimization before the code works is pretty useless.

Make it work, make it right, make it fast, in that order.

[K^2]

In GML, ^ is the bitwise xor operator, not exponentiation.

And it doesn't throw an exception on being used with a real variable, rather than integer? That's just silly. Good catch, though.

[Wheffle]

Good eye with the bitwise operator issue. I've been using GML for a while, super brain fart there. Thanks for all the helpful replies!

Edit:

So yeah, with that bitwise operator mishap, looks like gravity was linear based on radius. After fixing that and tweaking my gravity constant, everything seems to be modeling correctly! The problem at its heart seemed to be the lack of inverse square law being applied.

Edited November 13, 2015 by Wheffle

[problemecium]

HAH! CALLED IT!

In other news, I'm glad you got this working ^^

[Yourself]

And it doesn't throw an exception on being used with a real variable, rather than integer? That's just silly. Good catch, though.

Oh, that's not the end of the silliness with that language. In this case it does that because GML only has one numeric type and that's a real (specifically a double). As I recall there are only two types in the whole language: real and string. In fact, the data structures within the language aren't really data types, you're provided with a numeric ID (which is how all resources work in Game Maker) which is passed as an argument to various functions that manipulate the data structure. The IDs aren't unique and they're created sequentially, so the first instance of any resource you create has ID 1 (or 0, I can't remember). People can (and have) gotten into all sorts of trouble because of this.

Good eye with the bitwise operator issue. I've been using GML for a while, super brain fart there. Thanks for all the helpful replies!

Edit:

So yeah, with that bitwise operator mishap, looks like gravity was linear based on radius. After fixing that and tweaking my gravity constant, everything seems to be modeling correctly! The problem at its heart seemed to be the lack of inverse square law being applied.

You're not the first person I've helped with this problem. The first time I saw this exact issue was probably 14 years ago.

[cantab]

I did think it looked like a non-inverse-square pattern, having seen them in a flash gravity sim I played around it. But that oddity of syntax is devious.

[Yourself]

I did think it looked like a non-inverse-square pattern, having seen them in a flash gravity sim I played around it. But that oddity of syntax is devious.

Well, it's not the most devious. Syntactically GML is based largely on the C family of languages* (since we're talking syntax this means C, C++, Java, C#, etc.) and the ^ operator is bitwise XOR in all of those. In fact, very few (popular) languages use ^ as an exponentiation operator. More often than not it's bitwise XOR. I'm not entirely sure where the popularity of using ^ for exponentiation derives from exactly. I'd wager it comes from LaTeX since it uses the ^ operator is used to typeset superscripts. However, I don't know enough of its history to say whether that was inspired by something else.

*GML also supports an alternative syntax based on Delphi. This leads to some serious craziness since it understands Delphi operators for assignment (:=) and equality (=) as well as the C-style operators for assignment (=) and equality (==). These two syntax styles can be mixed interchangeably and often leads people to suggest that using = in an if statement would result in a bug:

if( a = b )

Since in most C-like languages, = is assignment. The thing about GML is that the meaning of = is context dependent. In that case, since it's parsing an expression, it interprets = is the Delphi-like equality operator and ends up working like its supposed to. But this can lead to other craziness:

a = b = 1

In C-like languages this would be multiple assignment; it would assign the value 1 to the variable b and then assign the new value of variable b (which is 1) to a. In GML that's not what this does. Unlike C an assignment in GML is a statement, not an expression (much like it is in Python), so assignments can only appear as a statement. The grammatical structure of an assignment is a variable name followed by an assignment operator (which can be either := or =) followed by any expression. Since the = operator acts as equality inside expressions, this code actually compares the value of b with 1 and stores the result into a.

[Wheffle]

It's one of those quirks that I've known about for years but continue to brain fart on. I might blame it on having to use scripting languages like Matlab for classes a lot right now, or it just might be my brain being stubborn.

I have a love-hate relationship with GML. My first programming was done with GML when I was young, rather sneakily teaching me some basics without me really knowing it. There's value in that. It does drive me nuts very often though. They have improved over the years though, as they now treat arrays more intelligently, and they actually have an 'undefined' data type, pointers (used with a very small set of functions related to graphics) and enums. Up an onwards! lol.

Anyway, I actually got a small 2D gravity simulation working that computes the orbit all at once KSP-style and allows you to perform prograde and retrograde velocity changes. I'm using very hacky methods for producing the orbital path, so it breaks in extreme orbits (steep falls and such), but it's a start! I was surprised actually how involved the math is for orbital mechanics, which made me take some shortcuts (I guess, rocket science isn't easy?). I might go back and do it properly if I feel up to it.

For anyone interested:

My Janky 2D Orbital Simulator (windows)