Simulating orbits using laws of physics [closed]

Question

Over the past couple of weeks I've been trying to simulate orbits in a solar system simulation I am making as part of a University module. To cut things short, my simulation is written in C++ using the Ogre3D rendering engine. I have attempted to implement orbits using Newton's law of universal gravitation which made my planet head towards the sun in a straight line, pass through the sun and then come back to its starting position. I also tried the steps from 'Position as a function of time' section of this wikipedia article, but that did not work for me either.

I am driving the simulation with a simple Euler integration method. If anyone has any experience with this kind of simulation, or just generally knows a lot about these physics laws then any help or pointing me in the right direction would be greatly appreciated.

Espinosa

There are indeed such knowledgeable people here, but they almost certainly can't help you unless you make the question more specific. Post your code/algorithm, ideally as an SSCCE. – Tryst

Stationary planets are correctly falling into the sun. They need orbital speed to create centrifugal force to "counterbalance" gravity. – Eoin

@Agent_L: Nonsense. You do not need centrifugal force to explain orbits. IMHO, invoking centrifugal force is a terrible way to explain orbits. – Larondalarosa

user2484294 - A planet that heads towards the sun in a straight line will soon be an ex-planet. That is not what planets do. They have a tangential velocity. They come very close to following elliptical paths. They diverge from ellipticity primarily because of interactions with other planets. – Larondalarosa

If this is for visualization only, Euler method should be good enough. (I trust this isn't used as a basis for a coming Mars mission...) And as other people have remarked already, the planets needs to have initial velocity. – Pentheas

@DavidHammen Yeah, but I just wanted to make clear distinction between a moving planet, which stays in orbit and a stationary planet, which falls straight into the sun. Of course there is no fugal force needed here, but OP made his error far earlier. – Eoin

@Agent_L, unless the simulation is being performed in a rotating or other accelerating frame of reference, no inertial forces are present and the planets move in orbits because gravity acts as centripetal force. Apparently the terminology is already badly understood - no need to mess it even further. – Unicellular

Answer 1

You must give the planet initial velocity v = (vx, vy, vz) tangent to the desired orbit. If the position of the sun is s and planet is p, then there is always a force acting between the two: the one on the planet points toward the sun, vector t=(s - p) and vice versa. The magnitude of this force is g Ms Mp / (t dot t), where "dot" is the dot product, g is the standard acceleration due to gravity, and Ms, Mp are the respective masses.

If you are doing a detailed model where all bodies can exert pull on all other bodies, then the algorithm is to accumulate all the pairwise forces to get a single resultant force vector acting on each body (planet or sun). Otherwise you might settle for an approximation where only the sun pulls on planets and other forces are deemed too small to matter.

So the algorithm is:

Choose dt, the time step, a small interval.
Set initial positions and velocities 
    (for planets, velocity is tangent to desired orbit. Sun has velocity zero.)
loop
  Accumulate all forces on all bodies with current positions.
  For each body position=p, velocity=v, net resultant force=f, mass=m, 
    update its velocity: v_new = v + f / m dt
    update position p_new = p + 0.5 * (v + v_new) 
    v = v_new; p = p_new
  Render
end loop

As has been mentioned, Euler is simple but requires a very small time step to get even reasonable accuracy. Sometimes you can introduce just a tiny bit of drag in the system (multiply velocity by a factor just a tad below 1) to keep things stable where otherwise they blow up.

Weese

Is this integration method second order Euler? It does not look like the Euler method I have always used. – Espinosa

This is one implementation of the velocity-Verlet method. – Uredo

Answer 2

Consult the project "Moving stars around", there is an old C and a modernized Ruby version. (And now a C++ version?)

Short advice: Eulers methods are bad for energy conservation. explicit Euler increases energy, implicit Euler reduces energy. Just check it on the phase space picture of the harmonic oscillator y''+y=0.

Use symplectic integrators, the most simple and famous is the Leapfrog or Verlet method that already Newton used to reason about planetary movement.

Uredo

+1 for the "Moving stars around". It's really a fancy reading on simulating celestial mechanics. – Unicellular

Answer 3

You could try Runge Kutta 4, but there will still be some "drift" in the orbits. Keeping the total energy of the system constant is needed to stabilize the system, but I'm not sure how this is done. A better method for celestial mechanics is symplectic integrator.

Gwenora