# How to simulate planets orbit curce around the sun ?

How to simulate planets orbit curve around the sun ?

Hi

I'm about to write a java program where I wish to simulate planets orbit curve around the sun, so that I can set the speed of a planet and mass of planet and sun and simulate the curves around the sun (just basic simulations). I only know high-school level physics but I know university level mathematics. I want to build my programs on newtons laws. I know that there are plenty of others programs that can do such simulations, but I wish to know the theory and math behind it so that I can build my program from ground up. I know newton laws of gravity, but I do not understand how to apply them to a simulation, and when I fool around with the formulas I manage to get different results depending on how I interpreted applying the formulas to a simulation. I hope that someone here will be so kind to explain me the theory of doing such a simulation or tell me where I can find books, documents, or websites about it. I have tried searching on the net, but have just got more confused. Thanks.

Last edited:

## Answers and Replies

Hi there,

Keplers laws of planetary motion will explain to you the physics behind the orbit, which I think is probably all you need to get a Johandle on the situation. excuse the cheesy pun.

http://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion

One thing I notice is that you say you can set the speed of the planet and the mass of the planet, but since planets travel in elliptical orbits, the speed of the planet will be a function of the distance from the sun. Setting the speed of the planet will not suffice for a stable model. There is a value called eccentricity which depicts how elliptical an orbit is. an eccentricity of 0 is a circular orbit, but you will come across all of this if you read through the link I posted.

DaveC426913
Gold Member
In terms of simulating:
Decide on your accuracy (some algorthms are harder than others but more accurate).
Decide if you want it 2D or 3D (recommended: 2D)
Decide if the orbiting bodies have significant mass compared to the parent (if no then your job is easer. If yes, you can simulate star clusters).

What you do is make nested iterative loops.

Plot all points in orbits, putting their xy coords in an array
Calc the force vector pulling on them from the parent (prop to mass, inverse prop to sqaure of dist)
From the force, calc the change in velocity.
From that, calc the new position.
Do this for all point, then update the display and repeat.

One thing I notice is that you say you can set the speed of the planet and the mass of the planet, but since planets travel in elliptical orbits, the speed of the planet will be a function of the distance from the sun. Setting the speed of the planet will not suffice for a stable model.

Thanks a lot for your reply. But I thought that the shape of the curve is depending on the total mechanical energy of the system E_mek which is equal to E_kin + E_pot (E_mek = E_kin + E_pot = 1/2*G*(M*m)/r

Therefore I thought that the curve form is dependent on the speed and/or the mass of planet (and/or the mass of the sun)

My plan was to be apple to simulate all sort of curves (comets and so on) including hyperbolic and parabolic curves with the speed and mass as starting conditions. But maybe that is possible with Keplers laws, even though I found it hard to see how to make a simulation out of it!

In terms of simulating:
Plot all points in orbits, putting their xy coords in an array
Calc the force vector pulling on them from the parent (prop to mass, inverse prop to sqaure of dist)
From the force, calc the change in velocity.
From that, calc the new position.
Do this for all point, then update the display and repeat.

... and if you stop at that, you will soon discover that your simulated planets have wandered off, or crashed into the Sun, or did something else equally bizzare. The worst possibility (as regards consequences of the error) is that results may look credible enough to pass cursory examination but be nevertheless wrong.

The reason for the above is accumulated systemic error. Increments (distances between the points on your graph) are by necessity finite. That means that the force, acceleration etc you calculated at the start of next step is wrong for the rest of that step. Even if you are crafty and use a central spot or several of them - that's still not enough. It will only reduce the error but not eliminate it.

To make sure your planets behave, you need to ensure that (at every step of your calculations!) they comply with conservation laws of energy and momentum. Force-adjust the results at every step to make sure your model at least doesn't break the natural laws. Not too much anyway :-)

In terms of simulating:
Decide on your accuracy (some algorthms are harder than others but more accurate).
Thanks a lot for your reply. But does that mean that some algorithms might give different shapes? I would to make it as precise as possible (only 2D for now).

What you do is make nested iterative loops.
Plot all points in orbits, putting their xy coords in an array
Calc the force vector pulling on them from the parent (prop to mass, inverse prop to sqaure of dist)
From the force, calc the change in velocity.
From that, calc the new position.
Do this for all point, then update the display and repeat.
So that should give the same result no matter the step size? If so then I'm not sure why that is, but then I have enough to make the program.

... and if you stop at that, you will soon discover that your simulated planets have wandered off, or crashed into the Sun, or did something else equally bizzare. The worst possibility (as regards consequences of the error) is that results may look credible enough to pass cursory examination but be nevertheless wrong.

The reason for the above is accumulated systemic error. Increments (distances between the points on your graph) are by necessity finite. That means that the force, acceleration etc you calculated at the start of next step is wrong for the rest of that step. Even if you are crafty and use a central spot or several of them - that's still not enough. It will only reduce the error but not eliminate it.

To make sure your planets behave, you need to ensure that (at every step of your calculations!) they comply with conservation laws of energy and momentum. Force-adjust the results at every step to make sure your model at least doesn't break the natural laws. Not too much anyway :-)

Thanks a lot for your reply. That was also what I have been thinking according to that the step size will influence the result. Do You know anywhere (documents, books or something) where it is described how to make the best possible model? (I'm a computer programmer and would like the very best model possible. My plan is to start in 2D (with a single planet) and then move to 3D and then up with a complete solar system simulator)

Last edited:
DaveC426913
Gold Member
... and if you stop at that, you will soon discover that your simulated planets have wandered off, or crashed into the Sun, or did something else equally bizzare.
I've built these several times and they are stable over long periods.

However, what I'm building are not models of solar systems, they're simply generic orbital paths of up to 20 or so stars in a cluster. There's no accuracy to any real system, but they are stable.
[/QUOTE]

DaveC426913
Gold Member
I would to make it as precise as possible (only 2D for now).
Oh. In your OP, you said "basic simulations", which I interpreted as "rough".

Oh. In your OP, you said "basic simulations", which I interpreted as "rough".
Actually I meant with a single planet orbiting around the sun, and when I get that working then with several planets in 2D and then hopefully ending up with a 3D solar system simulator (a solar system space travel program). I just haven't been able to find books and/or documents about how to make such simulations (also I would like to make a simulation where it is possible to add new objects).

Last edited:
I know that there exist ephemerides but I would like to make it as a simulation as I am very interested in astronomy and would like to understand the mechanisms in the solar system.

Last edited:
To make sure your planets behave, you need to ensure that (at every step of your calculations!) they comply with conservation laws of energy and momentum. Force-adjust the results at every step to make sure your model at least doesn't break the natural laws. Not too much anyway :-)

I already knew from experiments that the step size would affect the curve form. Can you (or someone) say more about how to force-adjust at every step (formulas, code fragments or so)?

Also what do you mean with not breaking the natural laws too much anyway? Is it not possible to make a perfect simulation, and if not then why not? And if it is possible to make a perfect simulation where can I find informations about it? I have tried to search for documents and at Amazon for books about it, but I haven't been able to find any documents or books about high quality computer simulations of planets orbiting the sun!

Also my questions is not based on passing any exams, but is the starting point in making a perfect planet simulation program (I will start in 2D with one planet), and therefore I would like to know everything about the physics related to it. I have a degree in computer science (and also I have taken a few courses in physics), so I know how to make graphics etc., but I need to understand the physics and mathematics of planetary computer simulations a lot better.

I want to build my programs on newtons laws.

As DaveC426913 already wrote you have to decide which numerical method you want to use. This may strongly influence the structure of your program. The best choice would be a symplectic multi-step integrator for second order IDEs with variable step width. But there are a lot of simple methods that give sufficient results. I prefer a Runge-Kutta-Nyström method.

Than you should think about the representation of the data. There are several different possibilities such as putting the data into different arrays or encapsulating them into objects. I prefer objects for each body and an array for all bodies. That makes it very easy to remove or add bodies even during a running simulation.

I would not recommend to limit the program to 2D first because it might be very difficult to extend it to 3D later. As you have to use vector algebra in both cases there is actually no reason for such a limitation.

As you want to extend your simulation to a space travel program you have to implement a full simulation of Newtons gravitational law with individual mass for different bodies including the possibility of external forces (e.g. from rocket motors). But you may use a hybrid simulator with Kepler for Sun and planets and Newton for everything else.

Is it not possible to make a perfect simulation, and if not then why not?

It is not possible because there is no global solution for n-body problems with n>2. Thats why such systems can only be approximated numerical and such methods have a limited accuracy.

Last edited:
Hi kjensen. I think http://www.cs.princeton.edu/courses/archive/spring11/cos126/assignments/nbody.html" [Broken] may be what you're looking for. I used it to write a java simulation which I later ported over to VB with added bells and whistles. I learned a lot and had a lot of fun with it.

Last edited by a moderator:

As DaveC426913 already wrote you have to decide which numerical method you want to use. This may strongly influence the structure of your program. The best choice would be a symplectic multi-step integrator for second order IDEs with variable step width. But there are a lot of simple methods that give sufficient results. I prefer a Runge-Kutta-Nyström method.
Thanks a lot for your reply. If I google "symplectic multi-step integrator" then I get a lot of interesting stuff (documents, books and so on). It might be what I am looking for! Do you know any books or documents where it is related to computer programming (any programming language). Thanks.

Hi kjensen. I think http://www.cs.princeton.edu/courses/archive/spring11/cos126/assignments/nbody.html" [Broken] may be what you're looking for. I used it to write a java simulation which I later ported over to VB with added bells and whistles. I learned a lot and had a lot of fun with it.
Thanks a lot for your reply. The link looks really great and even though I will build my program from ground up with a model as precise as possible then it might be worth figuring out the exercise (but it looks rather similar to another programming experiment I have been doing). But I think that the curve form will depend on the step size doesn't it?

Last edited by a moderator:
tony873004
Science Advisor
Gold Member
If you want to learn to walk before you learn to run, I'd start with Euler's method. It's the method DaveC described in post 3. It will give a very intiutive understanding of what you're doing. Leapfrog and RK4 are fine, but more difficult to code. At this point high accuracy is not an issure. You're creating a single planet moving around a stationary star, just to get a feel as to how to do it. So who cares if your simulation isn't accurate enough to predict solar eclipses decades into the future?

Actually Euler's method can be accurate enough to predict solar elipses decades into the future if you take a small enough time step. (delta t). The higher order integrators allow you to achieve good accuracy at larger timesteps, meaning it will take your computer less time to simulate your system decades into the future.

Later, you can upgrade to a higher order integrator. Since an Euler step is performed in the RK4 method, you'll be building upon what you've already done.

If you want to learn to walk before you learn to run ...
Thanks a lot for your reply. Of course you are right, and I will design my program so that I easily can change the 'calculation engines' along the way. But now when I have started a thread about it then I would like to know what the best possible method is, so that I don't have to start another thread later on.

Do you know any books or documents where it is related to computer programming (any programming language).

Unfortunately not. Maybe you will find corresponding sources in the references of your google-results. But you should also be able to write your program without literature. The program code will mainly result from the selected integrator and the representation the physical parameters. For example http://www.drstupid.de/Newton.html" is based on this Runge-Kutta-Nytröm integrator:

$\begin{array}{l} a_1 = a\left( {x_n } \right) \cdot \Delta t \\ a_2 = a\left( {x_n + {\textstyle{{\Delta t} \over 2}}v_n + {\textstyle{{\Delta t} \over 8}}a_1 ,v_n + {\textstyle{1 \over 2}}a_1 } \right) \cdot \Delta t \\ a_3 = a\left( {x_n + {\textstyle{{\Delta t} \over 2}}v_n + {\textstyle{{\Delta t} \over 8}}a_1 ,v_n + {\textstyle{1 \over 2}}a_2 } \right) \cdot \Delta t \\ a_4 = a\left( {x_n + \Delta t \cdot v_n + {\textstyle{{\Delta t} \over 2}}a_3 ,v_n + a_3 } \right) \cdot \Delta t \\ x_{n + 1} = x_n + v_n \cdot \Delta t + {\textstyle{{\Delta t} \over 6}}\left( {a_1 + a_2 + a_3 } \right) \\ v_{n + 1} = v_n + {\textstyle{1 \over 6}}\left( {a_1 + 2a_2 + 2a_3 + a_4 } \right) \\ \end{array}$

[http://theory.gsi.de/~vanhees/faq/gravitation/node62.html" [Broken]]

This algorithm already tells you something about the corresponding program:
As it is a 4-step algorithm intermediate data must be handled. For every step 4 accelerations (for different positions) must be calculated for every body and kept in memory. I solved this problem by an array of 4 accelerations for each body:

Code:
function Body(Mass,sx,sy,sz,vx,vy,vz,Name) {
this.m = Mass ; this.s = new Vector3D(sx,sy,sz) ; this.v = new Vector3D(vx,vy,vz) ;
this.F = new Vector3D(0,0,0) ; this.a = new Array ; for (var k=0;k<4;k++) this.a.push(new Vector3D(0,0,0)) ;
[...]
}

Within each time step the accelerations a - a are calculated using the first four equations of the algorithms. Than the velocities s and v are updated according to the last two equations using the accelerations calculated before. Each step of the algorithm must be calculated for all bodies at once. In doing so the force F is accumulated for each pair of bodies (in order to calculate each force only once).

With other algorithms and other data management you will get different code. If you select a numerical method and your favorite data management you may already have an idea about the structure of your program.

Last edited by a moderator:
If you want to learn to walk before you learn to run, I'd start with Euler's method. It's the method DaveC described in post 3. It will give a very intiutive understanding of what you're doing. Leapfrog and RK4 are fine, but more difficult to code.

You are right. But as a compromise I would suggest http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet". It is almost as simple to understand and to implement as Euler but it conserves the energy.

Last edited by a moderator:
Thank you DrStupid. Your code looks really great. I will give that approach a try!

D H
Staff Emeritus
Science Advisor

Thanks a lot for your reply. If I google "symplectic multi-step integrator" then I get a lot of interesting stuff (documents, books and so on). It might be what I am looking for! Do you know any books or documents where it is related to computer programming (any programming language). Thanks.
Symplectic multi-step methods are a big step for someone just starting out.

I suggest you start with something simple such as Euler's method or symplectic Euler, but design your system so that the most of the details of the integration technique are hidden from the uses of those techniques. This way you can start with something simple and step up to more advanced techniques as you learn.

Besides, accurate ephemeris models typically do not use symplectic methods. JPL uses a variable step-size, variable-order, Adams method to compute their Development Ephemerides (DE series; the latest of which is the DE421); the Russian Institute of Applied Astronomy of Saint-Petersburg uses a similar technique for their Ephemerides of the Planets and the Moon (EPM), as does the French Observatoire de Paris for their INPOP ephemerides. NORAD uses Stoermer-Cowell (aka Gauss-Jackson) to integrate the states of objects that are orbiting the Earth.

There are lots of things to learn about besides integration techniques.

You'll need to learn about how to represent time. Eventually you'll learn about why we have leap seconds. As a starter, I suggest you use TAI (International Atomic Time) as your time standard. If you later want to address relativistic effects you will need to move to a relativistic time scale such as JPL's T_eph or the IERS's Barycentric Dynamic Time. If you want to relate the calculated ephemerides with what you can see from the ground you'll also have to learn about things like Universal Time.

You'll need to learn about how to represent positions in space. There are many choices of reference frames. I suggest the IERS's International Celestial Reference System. If you want to relate the calculated ephemerides with what you can see from the ground you'll also have to learn about things like Terrestrial Reference Frames, and a whole lot about the Earth's rotation.

As a starter, I suggest you use TAI as the dynamic time for your system with coordinates in ICRF, and ignore all of the subtleties. You can add the subtleties as you learn more (and as you need them).

Some key websites:
US Naval Observatory Astronomical Information Center
http://www.usno.navy.mil/USNO/astronomical-applications/astronomical-information-center [Broken]

Standards of Fundamental Astronomy:
http://www.iausofa.org/

International Earth Rotation and Reference Systems Service:
http://www.iers.org

As far as integration techniques go, I suggest that you start with something very simple. The two simplest techniques are basic Euler and symplectic Euler. Basic Euler treats position and velocity as a six-vector and integrates according to
$$\vec u(t+\Delta t) = \vec u(t) + \Delta t\frac{d}{dt}(\vec u(t))$$
This is equivalent to propagating position and velocity according to
\begin{aligned} \vec x(t+\Delta t) &= \vec x(t) + \Delta t\,\vec v(t) \\ \vec v(t+\Delta t) &= \vec v(t) + \Delta t\,\vec a(t) \end{aligned}
Symplectic Euler merely reverses the order of the calculation:
\begin{aligned} \vec v(t+\Delta t) &= \vec v(t) + \Delta t\,\vec a(t) \\ \vec x(t+\Delta t) &= \vec x(t) + \Delta t\,\vec v(t+\Delta t) \end{aligned}

Basic Euler is incredibly inaccurate and unstable. Symplectic Euler is fairly inaccurate but is stable. The only reason to learn about basic Euler is that it is the basis for a lot of other integration techniques. The Adams family of integrators (which includes Stoermer-Cowell) and the Runge-Kutta family of integrators use basic Euler as the underlying mechanism for integration.

A couple of techniques that are a bit better than the Euler techniques are leapfrog, verlet integration, and Beeman's algorithm. A lot of galactic simulations, simulations in chemistry, and robotic simulations use these techniques or similar ones. Note very well: They don't use these techniques because they are good. They use these techniques because anything more advanced would kill the computer while dropping down to Euler integration would kill accuracy.

Beyond these, you'll find two key families, the Adams and Runge-Kutta families of integrators, that predominate in this field. (There are plenty more techniques than these, but these are the dominant ones). Neither Adams nor RK-type integrators are symplectic. The losses are small, particularly if you go to higher order integrators. One problem with symplectic integrators is that it is very, very hard to get beyond second order. It is fairly easy to come up with Adams and Runge-Kutta integrators that are well beyond second order.

Perturbation techniques open the door to a completely different style of integration. The Apollo program used an Encke-Nystrom integrator to accurately propagate the state of the Apollo spacecraft. An orbiting body will follow something that is fairly close to an elliptical path. An Encke-Nystrom integrator takes advantage of this fact, propagating the path as a small deviation from a simple elliptical path. The integrator has to re-anchor the ellipse to on occasion to keep the integration accurate.

Integrating something other than cartesian position and velocity opens the door to yet another style of integration. Variation of parameters techniques integrate Lagrange's planetary equations (or some variant of those planetary equations). In these techniques, it is orbital elements that are integrated rather than position and velocity. The planetary equations provide the partial derivatives of these elements with respect to the perturbing forces. A somewhat related technique is to use the Karhunen-Loève transform. An elliptical path falls out directly from this transform.

Last edited by a moderator:
tony873004
Science Advisor
Gold Member
Thanks a lot for your reply. Of course you are right, and I will design my program so that I easily can change the 'calculation engines' along the way. But now when I have started a thread about it then I would like to know what the best possible method is...
The best possible method might be one that you invent :)
I've done exactly what you've proposed: design code so that I can easily change 'calculation engines'. The code for these integrators is available on my website ( www.gravitysimulator.com ). If you have a Windows computer, you can download it and compare the different routines to see their strengths and weaknesses.

Thank you very much D H. Your answer seems very deep and profound, and it looks as a very good 'recipe' for what I need to make the program. I think that I better make the program in C/C++ as it otherwise will get to slow when I extend it to higher precision! Thanks.

The code for these integrators is available on my website ( www.gravitysimulator.com ). If you have a Windows computer, you can download it and compare the different routines to see their strengths and weaknesses.
Your program looks great and also it runs under Linux with Wine even though i couldn't find the Integrator menu that you have shown on the picture. Also I couldn't find any source code on your website. I personally prefer open source systems and also my program when it is finished one day will be released as open source.