Modeling stellar orbits

Mike Mayo · 2010-12-12T12:37:34

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_gravitationalong 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]

Math and Physics Programming

Started by ChugginWindex November 15, 2010 06:15 PM

94 comments, last by taby 13 years, 4 months ago

ChugginWindex

187

Author

November 19, 2010 08:56 PM

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.

h4tt3n

1,974

November 20, 2010 12:51 AM

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

ChugginWindex

187

Author

November 20, 2010 10:28 AM

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]

h4tt3n

1,974

November 20, 2010 11:06 AM

Quote:Original post by ChugginWindex
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?

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

ChugginWindex

187

Author

November 20, 2010 12:43 PM

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);}

h4tt3n

1,974

November 20, 2010 01:20 PM

All looks ok at first glance. What exactly does the ApplyImpulse function do? Remember we're applying *forces* here, not impulse. If you've been doing so, then that would explain the sim failure.

Cheers,
Mike

[Edited by - h4tt3n on November 20, 2010 3:20:16 PM]

JoeCooper

350

November 20, 2010 01:21 PM

Quote:I'm noticing that all my stars eventually fall into the middle star now...perhaps I need to tweak my gravitational constant?

Does this happen when you only have two stars - a center one and a satellite?

Systems are often not stable. The Solar System has had a very long time to settle out, and there's still evidence of a messier early history (e.g. formation of Moon) and the planets still effect each other (Jupiter effects our orbit and climate in 100,000 year cycles).

The satellites of the Jovian system are locked in resonance with each other, so actually achieving stable orbits with many bodies is probably tricky.

ChugginWindex

187

Author

November 20, 2010 01:36 PM

ApplyImpulse is actually just bad language on my part. It's been about a year since I took Physics and I didn't feel like looking up terminology when I originally wrote that function (impulse is the change in acceleration isn't it?).

Either way, the ApplyImpulse() function actually just adds the forces that it is passed to the total X and Y forces on the object that frame. I should probably just rename it.

It's come quite a ways (technically) since when I uploaded an EXE, so here's a newer version if someone wants to take a look at the behavior that I'm talking about: Download

My apologies to all of you non-windows users, I'm Win7 only at the moment but if someone's really interested I've got no problem putting up a link to download the source as well, especially if it helps me solve this problem!

Also JoeCooper, it happens with as many stars as I put in there, whether its 2 or 20. Run the exe and you'll see for yourself.

Thanks for all the help from you guys too, I appreciate it.

JoeCooper

350

November 20, 2010 01:44 PM

I have two here - spawned mass of about 130 and 5, respectively - and it just keeps goin round and round in what appears to be a stable orbit. There's even a coresponding wobble on the parent star, which is neat. Should it take a very long time?

But I do have the previous version.

Maybe some test rig could show useful info. Autospawn two stars in what should be a stable orbit, have a monitor of the satellite's energy.

ChugginWindex

187

Author

November 20, 2010 01:47 PM

That's actually some stuff I'm currently working on. I need more information visible to me as the sim runs. But yeah you should try it with the new download I put up to see what I mean. The version you have is actually just Euler's integration method, whereas the new one uses velocity verlet integration.

EDIT: ...the wobble. The wobble should be just that, right? A wobble? In the newest version the stars don't just wobble, they end up shifting ever onward slowly in the same general direction as the "satellite" star's initial velocity when it was spawned. If I wasn't losing energy, the star would be staying in relatively the same place, correct? But just with a little wobble? Perhaps that's another sign that something is wrong...

Modeling stellar orbits

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Modeling stellar orbits

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines