Create an Orbital Simulator

Redbelly98

Maybe the best thing right now is to use some basic algorithm, that you can easily understand, to get something working first. Then worry about making it more accurate or complex as a later project.

bkelly

I found a simulator here: http :// kornelix . squarespace . com / galaxy /
(it seems I don't have URL priviledge, but google on "stellar simulation")
It runs under Linux so I fired up that computer and installed it. It runs well, but not much in the way of features. I haven't reviewed the code closely, but it seems to just use the F=MA type of approach. If I can get version 1.0 running, then I can look into better accuracy for close approaches.

There are two primary segments of this projects
Calculating the positions
Displaying the field
In this forum I want to discuss the calculations.

What type of coordinate system would be best, polar, Cartesian, or something else I don't know about? It seems to me that Cartesian would be better. If for no other reason, I have never messed with polar coordinates. But if polar would be much better, I would like to know.

If Cartesian, I suspect that it would be easier to deal with if the origin is in one of the corners rather than in the middle. My first thought would be the lower left corner. Thoughts about selecting the range: Astronomical Units, millions of miles, something else?

BTW: This site seems to respond much faster under Linux with Firefox than under XP with IE-7.

Phrak

bkelley's the one you want to talk to about code. I wrote mine a few moons ago before they has invented rocks. Did you solve or notice the walk-away problem where circular orbits are really outward spirals? You can check if it's happening with a running total-energy value.

I hope I didn't PO lonedragon excessively...

Redbelly98

bkelly,

In Cartesian coordinates it's easier to convert from the (x,y) calculated positions to displayed pixel locations, as only simple scaling factors are involved. Might also be easier to put the origin at pixel (0,0) in the display in the upper left, but if you don't it's still pretty straightforward to rescale from (x,y) to pixel location.

I remember an algorithm a former physics teacher gave me. Better than Euler, not as accurate as Runge-Kutta, but pretty straightforward in terms of understanding it.

In 1-D it works like this ...

Do a Taylor series expansion to find [tex] x_{n+1} [/tex] given [tex] x_n [/tex] , where [tex]x_n = x(n \cdot \Delta t)
[/tex] :

[tex]
x_{n+1} = x_n + v_n \cdot \Delta t + \frac{1}{2} a_n \cdot \Delta t^2 + \frac{1}{6}x'''_n \cdot \Delta t^3 +[/tex] ...

Similarly,

[tex]
x_{n-1} = x_n - v_n \cdot \Delta t + \frac{1}{2} a_n \cdot \Delta t^2 - \frac{1}{6}x'''_n \cdot \Delta t^3 + [/tex] ...

where [tex] v_n, a_n, x'''_n [/tex] are the velocity, acceleration, and 3rd derivitive of x, respectively, at time [tex] t= n \cdot \Delta t [/tex]

Add these two equations:

[tex]
x_{n+1} + x_{n-1} = 2 \cdot x_n + a_n \cdot \Delta t^2 +[/tex] ...

or

[tex]
x_{n+1} = 2 \cdot x_n - x_{n-1} + a_n \cdot \Delta t^2 + [/tex] ...

Note the cancellation of the odd-order terms, so that the error is actually of order [tex]\Delta t^4[/tex], not [tex]\Delta t^3[/tex]. All you need to calculate [tex]x_{n+1}[/tex] are the previous 2 values of x, as well as the acceleration at time [tex]n \cdot \Delta t[/tex], which you calculate from the Forces and masses in the system.

In 2-D you'd need to calculate the x- and y-components of the accelerations, then apply the algorithm to find [tex]x_{n+1}[/tex] and [tex]y_{n+1}[/tex] .

There might be a common name for this algorithm, if anybody knows it feel free to post the info.

Redbelly98

I didn't check circular orbits, that would have been a good test. What I did notice was that if two objects got very close, they could suddenly go flying off the screen at incredibly high speed.

p.s. I just looked over my code, and found that I used the algorithm I described in message #21. Somehow I thought I later changed this to Runge-Kutta but apparently not.
I can't figure out, from a quick glance at the code, just what I'm doing to the force to make it "behave" at close range. I would have to find some 3-year-old handwritten notes that I have somewhere to figure that out again.

Shukie

Anyone feel like sharing one of their simulators? Would be interesting.

Phrak

The guys at NASA, NSA, and whoever, should have devoted quit a few man-years at this. So they may have found a different method than interating t, and updating V_i and X_i. There's the Hamiltonian, if you've heard of it. It's an alternative, equivalent to Newtonian mechanics. What I'm trying to say, is that it may be possible to iterate values in some other mathematical space, rather than velocity and position. Angular momentum and Kinetic energy perhaps; I have no idea really. Just offinging possibilities.

D H

This is a direct consequence of using the Euler method. The Euler method displays this problem even with two bodies and even without any eccentricity.

The method you described is the basic verlet method. As you noted, it is better than Euler. An even better approach and still very simple to understand is Heun's method.

The cited simulator updates the velocity and then updates the position using the updated velocity. This is not a very good integrator. A simple improvement is to update the position using the average of the velocity before and velocity after the velocity update. This simple improvement is Heun's method, a second-order Runge-Kutta integrator. Heun's method is considerably better than the midpoint method, another second-order Runge-Kutta integrator.

*All* of the Runge-Kutta techniques as well as the verlet method "just use the F=Ma" type of approach. Some just do it a lot better than others.

That is always a good test. A simulator that can't do a very good job reproducing a circular orbit isn't worth very much at all. Highly elliptical orbits and hyperbolic orbits are also good tests because the close approaches will stress even good integrators.

D H

NSA, as far as I know, doesn't address this issue. OTOH, NASA certainly does. Hamiltonian-based techniques are called symplectic integrators; google this phrase for more info. The advantage is that these techniques inherently do a better job at conserving energy than do non-symplectic integrators (e.g., the Euler method and the standard RK4 method). The basic verlet method is a symplectic integrator, for example, and it does a much better job than does the Euler method even though both evaluate state derivatives once for each propagation interval. I believe JPL uses symplectic integration techniques for developing ephemerides.

There are many, many other orbital propagation techniques. One approach is to use Keplerian elements (e.g., [itex]a,e,i,\Omega,\omega,M[/itex]) rather than position and velocity. In the case of two bodies, all Keplerian elements but the mean anomaly are constant and the mean anomaly is a linear function of time. Propagation is easy and exact. Making sense out of the elements (e.g., plot positions on a computer screen) takes a little work, but not too much. Now consider a stellar system which comprises a single star and several planets. The planets' orbits can still be described in terms of Keplerian elements, but these now exhibit time-varying behavior. The variations are small, so perturbation techniques can be put into play. NASA used such an approach during the Apollo era, and still uses it in some arenas today.

Phrak

I hadn't considered perturbations from elliptical orbits, but then again, I haven't considered orbital simulations in a long while.

An interesting possibility to avoid close-approach divergence for any pair of bodies, might be to switch to elliptical orbits + perturbations (for the pair only) based upon some nearness condition, then switch back. It may even be possible to survive a direct hit where the mass centers pass through one another.

Phrak

In the above I failed to mention the possibility of close encounders of hyperbolic and parabolic trajectories where these might also be degenerate (head on collisions).

lzkelley

I'm also trying to make an n-body orbital simulator, its going quite well - i've recently upgraded from euler to fourth order runga-kutta.
I've been using Cartesian coordinates for all of my work so far, and for the actual calculations i plan to continue - but i need to be able to use keplerian orbital elements instead of state vectors for initial conditions. I was wondering if anyone out there had experience or direction on how to convert between orbital elements and state vectors besides drowning myself in geometry. Thanks!

csprof2000

Regardless of what everyone is telling you, your method will work fine. Just be prepared to use very small time steps if you are using Euler's method.

Besides inaccuracies that can occur due to the approximate nature of the solution, the near-collision "explosion" can be solved by simply defining the force piecewise. Calculate the distance between each pair of bodies (this must be done anyway) and see if it's less than some minimum value (this can be as simple or as complicated as you like). Then decide whether or not you want to "do" the force there, and act accordingly. Otherwise, you could use a similar method to use conservation of momentum to have them collide (elastically or inelastically) like billiard balls. Or have them lump together and run until all masses have lumped together.

There's a lot of potential for fun and good programming in what you're talking about.Good graphics can be done in Java, C, C++, C#, Visual Basic, etc. For 3-D, you may even look into DirectX or OpenGL.

The programming won't be hard, but there are some things you can do to make it easy on yourself. Look for libraries with functions / objects you can use. I recommend taking the opportunity to explore OO approaches, if you're not already familiar with these. You can make objects for vectors and massy particles and this will greatly ease the notation. Have fun with it.

Ulysees

May I ask, for a given accuracy of 10^-5 say, and for a given orbit lasting 1 year, which integration method is the fastest?

csprof2000

I would suspect that you can't do much better than the fourth order adaptive stepsize Runge Kutta algorithm.

Ulysees

Does that mean repeating with a smaller time interval until all 3 steps of 4th order Runge-Kutta agree within an error margin?

csprof2000

I'm not sure precisely what happens, but I'm pretty sure that a timestep is taken with one RK method and then another is taken with one which is of one higher order, and if the difference is within the error bars, it goes with it. Otherwise, it tells you what the timestep should be to get your accuracy. It has the benefit of remaining as accurate as you tell it (within reason) and not being *too* when the forces are negligibly small (for instance, why use a timestep of h = 1.0e-10 both when the force is large and changing quickly AND when the force is next to nothing and staying put?)