Modeling stellar orbits

Recommended Posts

I'm toying with simulating stellar bodies in my spare time right now, and I wanted to be able to do orbits with actual gravity such as in newton's equation mentioned here:

http://en.wikipedia.org/wiki/Newton's_law_of_universal_gravitation

along with F=ma of course. It all appears to work correctly, except that I can't achieve an actual orbit. I've gotten plenty of close calls, but in the end the minor stars slingshot around the major stars and end up flinging off into space. I logged some of the accelerations on the bodies that are being slingshot around and found that what appears to be happening (obviously) is that they are picking up acceleration on their way around the star, and can't be slowed down in a reasonable time. The values for their accelerations usually become steady around 1/10^-7 units per second (I'm just using pixels as it is in 2D) and VERRRRRY SLOWWWLY start to turn around and head back toward the major star(s). I have a feeling that this would eventually perform a complete orbit, but it would probably take me years of waiting to see it happen. I'm wondering if maybe I'm doing something wrong or leaving a factor out that contributes to the gravitational pull so that stars can orbit in a more observable manner. Here is the code I'm using for the two major functions dealing with this every frame:

//updates a body's internal statevoid Body::update(float delta) {	m_velocity.x += (m_force.x * delta) / (double)m_mass;	m_velocity.y += (m_force.y * delta) / (double)m_mass;	m_visual.Move(m_velocity.x * delta, m_velocity.y * delta);	m_location.x += m_velocity.x * delta;	m_location.y += m_velocity.y * delta;}//calculate and apply the force due to gravity void Body::applyPull(Body& otherBody) {	double forceOfPull = (GRAVITY_CONSTANT)*(double)(otherBody.getMass()*m_mass);	auto dist = this->distanceFrom(otherBody);	float forceAngle = std::atan2(dist.y, dist.x);	forceOfPull /= (dist.x * dist.x + dist.y * dist.y);	double forceX = forceOfPull * std::cos(forceAngle);	double forceY = forceOfPull * std::sin(forceAngle);	//this can be sped up if needed, we dont actually need to get the angle as we have the ratio	otherBody.addImpulse(forceX, forceY);}

EDIT: Okay, I started zeroing my m_force vector at the end of each update() call so that the objects didn't GAIN force each frame and now I'm getting orbits finally which makes me very happy. I'm still very interested in what people may have to recommend me changing / adding to this model so far.

So, anyone have any ideas what I can do to improve this? I know there's a LOT more that can be factored in, and many equations that deal with this sort of thing, but I wanted to start simple so if there's a quick recommendation that will help me solve my issue I would appreciate it.

[Edited by - ChugginWindex on November 20, 2010 2:59:51 PM]

Share on other sites
h4tt3n    1974
Hi there,

First off, I'd recommend you to get rid of all trig functions like atan2, sine, and cosine, and instead use vectors. This minimizes the calculation cost to just one cpu-heavy function call - sqrt() - per body pair per cycle.

Besides that, I'd recommend you to make it deal with more than two bodies. This allows for much more interesting simulations :-)

When you've handled that, implementing a better integration algorithm is the thing, but let's leave that for now.

If you'd like some very simple code samples to get inspired by, I don't mind digging something out for you. Have fun!

Cheers,
Mike

Share on other sites

I'd love any samples you could give me. I'm currently at uni for a Software Development degree but my physics is less than perfect so any inspiration would be helpful. Specifically I'm intrigued by what you recommended I should do during body interactions, such as dealing with more than two at a time. Currently (poorly) I am just iterating a list of bodies twice, and having each one apply pull on each other one, and I'm actually doing the same for collision as well. I'm well aware of how terrible this is but I was afraid if I didn't get SOMETHING running I was going to lose motivation like I have in the past. Now that I'm more motivated I'd like to revise this a lot, but I'm stuck on how I can make it more efficient. I figure some sort of BSP tree is in order, but I'm stuck on how I partition correctly right now because I don't technically define any boundaries to my world. Also, the whole iterating-every-body-twice approach seemed really inefficient when I was just dealing with collision as many objects would be far enough apart that it was useless to check them, but now that I'm having every body apply their force of gravity to every other object, is it really that costly to just check for their collisions while I'm at it?

Other than that, I'm interested in how you would recommend I use vectors for my acceleration calculations. In my code that I posted there's a comment I wrote to myself about how I didn't think I really needed to do the trig functions you were talking about, because I already have the tangent of my angle (y / x). If this is all you were talking about, then yes I planned to revise that soon I just didn't want to think a lot when I wrote it initially because it was very late. If there's something else you are thinking then I would be interested to hear it.

Thanks again!

Share on other sites
JoeCooper    350
Are you thinking of making a game with this?

Share on other sites
I haven't decided yet. I've never made a game before so I'm not sure if I'm in a position to think up a bunch of rules and such to make it into a REAL game. But I AM interested in the concept a lot. I just don't really know...

To be honest, this is all because I played with one online that was written in flash or something similar I believe. I got frustrated with all of its shortcomings and lag and felt that I'd take a shot at making one in native code that would be more responsive. Not to mention to have something I could add features to I felt the original was missing.

If it turns into a game that'd be cool, but I don't have a lot of ideas at the moment for how to make it have lasting appeal beyond its interest as a simulation, and it's not really a priority at the moment.

Share on other sites
Numsgil    501
The issue is that when you numerically integrate the motion (ie: when you turn acceleration into velocity into changes in position), there's some error involved. The tiny bit of error either adds or removes a bit of energy from the system, and the result is that the planets either spiral into their stars or fling out to the far reaches of space.

You can't really 100% solve this issue, but you can make it stable enough by toying with an adaptive step size. Basically, step your simulation forward by two half steps, and compare it to stepping forward with one whole step. If the difference in velocity/position is too large, divide your step size by two and try again. If the error is low, try doubling the step size and try again. The goal is to find the "goldilocks" zone of step size that is as large as possible while still providing stable behavior.

Obviously there's some added computational cost, since you're basically stepping the simulation at least 3 times per actual simulation step (potentially much more), so it may or may not be appropriate for what you have in mind.

Share on other sites
JoeCooper    350
You might go have a look at what "Orbiter Sim" does about this. Range Kutta something something. The fellow put a lot of work into making it reasonably accurate from hi-speed, tight LEO orbits to the rest of the Solar System.

Quote:
 I'm not sure if I'm in a position to think up a bunch of rules and such to make it into a REAL game. But I AM interested in the concept a lot. I just don't really know...

I've put a lot of thought into building a game like this, but I don't have time at the moment to start or join a project.

But if you're interested, I can dump my ideas here and if you decide to go ahead with a game, you can read through to see if there's anything useful.

Share on other sites
I don't see how posting your ideas would hurt ;)

However like I said, at this point I'm really interested in what improvements to the core functionality I might want to look at. The physics behind it and such. I'm also interested in any other projects like this that people here might either have going or have in a back room collecting dust. I'm not saying I want people's code, but I'd like to see what others have done.

Share on other sites
JoeCooper    350

I spent some time drafting a design but there were a bunch of projects and that one got de-prioritized. I've been experimenting with other things. Nevertheless, here's what I had come up with. Won't hurt, but it's a bit of typing so I thought I'd ask.

*Orbital maneuvers are counter-intuitive. I've found that games like KSpaceDuel can be rather hard to control. A common solution is to remove gravity. Instead, the solution I had was to remove rotation control; your vessel points "forward" while goin' round and round, and you can either accelerate or de-accellerate. Most orbital maneuvers can be done this way. So you have a speed up and speed down button. Having it draw the vessel's predicted path ahead for half an orbit could give instant feedback when you press either button, so that the controls are easy to master.

*Some games go the extra mile and remove combat maneuvering altogether. This is because the simplicity renders combat maneuvering superfluous. IRL, the need to perform combat maneuvering stems from all the complexities...

e.g. the AIM-9L air to air missile can only lock on to an enemy from his rear, so it can see the engines.

There are monochromat IR seeker heads that can't differentiate between engines and flares, and more expensive dichromats which can, and require more advanced electronic countermeasures to confuse, or even lasers to blind them.

AIM-7 needs the target to be continually painted by the jet's radar...

Aim-9X is a heat seeker that can get a lock head-on, so if you're fighting someone who has a much smaller budget and has only rear-IR, than the best place to be is _in front of him_, where your missile has a fundamental physics advantage.

When we move to space, some of those advantages are gone, but some issues are added, like delta-v to reach the target, kinetic energy on arrival (anti-sats like SM-3 use kinetic warheads) and so on. A bullet will have more kinetic energy on impact if fired down the gravity well than up. Might not be enough to matter IRL, but of course you can take all the creative liberties you want.

SO. When you've clearly define the weapon's abilities and constraints, you've defined good positions and bad positions.

Then you can dog-fight for a missile lock.

I think that minimally qualifies as a game. There's a lot further to go to make it interesting though. It'll come up on my list sooner or later and I'll flesh it out...

I had a few visions for the facade. One was to do it in the Saturn System between a nuclear-electric vessel and a lox-methane. The other was a duel ("for some reason") in LEO with either a 2020 or 1960s facade. It's _mostly_ the same tech, but the 1960s version would have a guy next to you with a notebook and pencil as your flight computer.

Of course, IRL, space fighters are a silly idea, but with some careful cherry-picking, the complexities of RL can allow for an interesting game.

[Edited by - JoeCooper on November 17, 2010 8:46:58 PM]

Share on other sites
Wow, that's quite a lot of ideas!

Thanks for that, whether you go through with it or not, it's useful to have some ideas on the back burner since, like I mentioned, I've never done a game before in my life...

My code is crap since I wrote it without a care in the world for efficiency or readability, but while I'm refactoring it all I figure I might as well upload my little orbital app to see if anyone's got some ideas after seeing it. I'd be interested in ideas to either add features to the simulator itself. Also, if anyone who's possibly more familiar with the physics and nature of planetary orbits than I am manages to look at this, I'd like to know if it looks believable. I know it's not going to be 100% accurate, but if the orbits look like I'm doing something wrong in my physics I'd like to hear it!

Otherwise, enjoy and I look forward to any early comments before the project gets any larger!

Share on other sites
JoeCooper    350
I fetched it thinking it was source code but its a program. My Windows computer is boned, I'll look at it after I'm done reinstalling everything.

Oh, also, while you're playing with it you might throw in atmospheres. Aerodynamic drag is a big deal. Space stations are kept at low altitudes, in the thin upper atmosphere, where the wind drag brings down debris and destroys meteoroids. Much safer. But if they don't supply it with fuel continually, it'll sink and the headwind will destroy it within only half a year.

Share on other sites
Hmm, atmospheres might be interesting once I get to a more detailed visual representation of these stellar bodies. You'll see what I mean if you open the program haha. Also at the moment my bodies are basically just stars, I'm working to set up a little polymorphic system so that bodies can be either gasseous or solid, and then doing stars and planets like that to give more diversity. That way I could handle collisions differently even if I wanted to. When I get that far, atmospheres will probably be a must. It's a very cool idea.

Share on other sites
Emergent    982
Quote:
 Original post by NumsgilThe issue is that when you numerically integrate the motion [...]You can't really 100% solve this issue, [...]

Two words: variational integrator. The intro chapters to Matt West's Thesis are a good starting point.

The downside is you get an implicit scheme, but that's really not so bad.

Share on other sites
h4tt3n    1974
Quote:
Original post by Emergent
Quote:
 Original post by NumsgilThe issue is that when you numerically integrate the motion [...]You can't really 100% solve this issue, [...]

Two words: variational integrator. The intro chapters to Matt West's Thesis are a good starting point.

The downside is you get an implicit scheme, but that's really not so bad.

Well, maybe that's a little complicated for a first simple experiment. Anyway, your system won't loose or gain energy if you use a symplectic integrator. A really good simple one is velocity verlet. I've implemented some higher order symplectic integrators too, and they did a really good job.

ps. Haven't forgotten about those code samples - still looking for them :-)

Cheers,
Mike

Share on other sites
Alright, I've got velocity verlet integration in place instead of Euler's integration. I'm noticing some strange behavior however. My major stars are being thrown around a lot more by their minor stars than they used to be. Maybe this is accurate? I don't really know. I just see them moving around with more speed now...

I saw that there's a similar integration technique for position that does not use velocity at all. I tried implementing and switching between updating position by velocity, and by this new verlet integration and saw pretty much the same thing. Should my major stars be moving more with velocity verlet integration than they were with Euler's?

In case I'm just messing up, here is my new code for updating everything:

//update velocity using velocity verlet integration based on the accelerations of both this time and last timestep	auto lastAccel = m_lastAppliedForces;	lastAccel.x /= m_mass;	lastAccel.y /= m_mass;	m_velocity.x += 0.5 * (lastAccel.x + m_appliedForces.x / m_mass) * timeDelta;	m_velocity.y += 0.5 * (lastAccel.y + m_appliedForces.y / m_mass) * timeDelta;	m_position.x += m_velocity.x * timeDelta;	m_position.y += m_velocity.y * timeDelta;	//since this is called at the end of all interactions on the body this frame, clear out the forces	//after updating them as the last timestep's forces so they can be recalculated each frame	m_lastAppliedForces.x = m_appliedForces.x;	m_lastAppliedForces.y = m_appliedForces.y;	m_appliedForces.x = 0;	m_appliedForces.y = 0;

In that page, they talk about a technique that looks awfully similar to my original Euler method, only they update velocity by half of the delta time, then update the position, and then update the velocity again by the other half delta. Is this better than the verlet technique I am currently using? They claim the accuracy and error is better.

[Edited by - ChugginWindex on November 18, 2010 11:19:11 AM]

Share on other sites
Emergent    982
Quote:
 Original post by h4tt3nps. Haven't forgotten about those code samples - still looking for them :-)

From me? Of a variational integrator? If so, then I had forgotten, but I wouldn't mind digging one up...

Share on other sites
h4tt3n    1974
Quote:
Original post by Emergent
Quote:
 Original post by h4tt3nps. Haven't forgotten about those code samples - still looking for them :-)

From me? Of a variational integrator? If so, then I had forgotten, but I wouldn't mind digging one up...

Oh, that was for ChugginWindex - I promised him some code samples. But hey, throw in whatever you got as well. I'd like to see how that variational integrator works, haven't played around with those yet. :-)

cheers,
Mike

Share on other sites
h4tt3n    1974
Ok, here's something...

http://www.jernmager.dk/stuff/gravity_code_examples.zip

They're very simple and commented line by line. Sorry for the basic dialect - I really should redo these in c++ once I get the time :-/

Cheers,
Mike

Share on other sites
Oh hey! I stumbled across those myself not too long ago. The methods you used seemed very similar to what I had at the time.

Share on other sites
JoeCooper    350
I got around to running this. It's a really fun toy! I wish I could center on an entity; they drift around.

Share on other sites
Yeah, I wanna add the ability to lock on to a star at some point. I'm working on trying to get a rewind ability working so that I can aid myself in figuring out why every once in a while my physics go haywire.

Share on other sites
h4tt3n    1974
Try using the velocity verlet instead, here in pseudocode:

// half velocity update ( old accel)Velocity += 0.5 * Accel * delta// acceleration updateAccel = Force / Mass// half velocity update (new accel)Velocity += 0.5 * Accel * delta// position updatePosition += Velocity * delta + 0.5 * Accel * delta * delta

If it still goes haywire, the error is somewhere else in your code.

Cheers,
Mike

Share on other sites
Thanks for that h4tt3n, your 2-step velocity update is actually the same as what I had in my last code sample, but your position update is more accurate than mine, and it helped significantly actually. Orbits seem to be much more stable!

EDIT: Actually, while it seems more "stable" in the sense that things look more realistic, I'm noticing that all my stars eventually fall into the middle star now...perhaps I need to tweak my gravitational constant?

[Edited by - ChugginWindex on November 20, 2010 10:28:23 AM]

Share on other sites
h4tt3n    1974
Quote:
 Original post by ChugginWindexActually, while it seems more "stable" in the sense that things look more realistic, I'm noticing that all my stars eventually fall into the middle star now...perhaps I need to tweak my gravitational constant?

Hmmm, that's strange. Your G constant probabl hasn't anything to do with that. It sounds like your system is loosing energy. Could you in pseudocode describe exactly how you loop through your stars in order to obtain a gravitational pull for each unique pair of stars. It should look something like this:

for a = 1 to number_of_stars - 1 for b = a to number_of_stars  // calculate pull bewt. stars a & b next bnext a

(look up "hand shake problem")

cheers,
Mike

Share on other sites
Here's pseudocode for the update function of my BodySystem class, which represents the "universe":

for(bodyIter = allBodies.begin(); bodyIter != allBodies.end(); bodyIter++) {   bodyIter->update(deltaTime); //this calculates the acceleration, velocity etc. EXACTLY as you posted before h4tt3n, at the moment}//now I calculate physicsfor(bodyIter = allBodies.begin(); bodyIter != allBodies.end(); bodyIter++) {   for(bodyIter2 = allBodies.begin(); bodyIter2 != allBodies.end(); bodyIter2++) {      if(bodyIter == bodyIter2) continue; //skip itself      bodyIter->applyGravitationalPull(bodyIter2);      //check for collisions, not important to note here I think   }}

For good measure, here's my applyGravitationalPull() method inside the Body class

void Body::applyGravitationalPull(Body& otherBody) {	double forceG = (GRAVITY_CONSTANT) * otherBody.getMass() * m_mass;	vector2f dist = GetDistanceBetweenBodies(*this, otherBody);	forceG /= (pow(dist.x, 2) + pow(dist.y, 2)); //finish gravity calculation	        //normalize the distance vector	double absDistance = sqrt(pow(dist.x, 2) + pow(dist.y, 2));	dist.x /= absDistance;	dist.y /= absDistance;		//generate the final force vector	otherBody.applyImpulse(dist.x * forceG, dist.y * forceG);}