N-body simulation test with sun, earth and moon

[Heimdall]

Hi !

first, excuse my english, I'm just a poor french student... and we are so bad in languages here... anyway..

I'm trying to make a little program which could calculate the motion (and some other little things) of n bodies linked with gravitationnal interaction.
(my final aim will be to simulate planets orbits arount binary stellar systems)..

In fact I have already wrote the code, but I don't understand some things... that I hope you'll be able to explain to me ;)


I use the 4th order runge-kutta method to solve the equations of the motion, with a fixed step for the moment (I'll write a method with an adaptable step later).

I have tested my code with the three bodies system : Earth moon and sun. using the JPL orbital parameters for my starting conditions.

I suppose that at time t=0, the orbits are elliptic so I can calculate positions and velocities with the semi-major axis 'a', the excentricity 'e' and the inclinaision angle 'i'.

At time t=0 I put the planets at their apogee with the formula Rmax = a*(1+e).

so for example the cartesian coordonnates for Earth at t=0 are :

x_t = a*(1+e)
y_t=0
z_t=0

for the moon, the inclinaison i is 5.16 deg, so :
x_m = x_t + a_m*(1+e_m)*cos(i)
y_m = 0
z_m = a_m*(1+e_m)*sin(i)

(a_m and e_m are respectively semi major axis and excentricity of the moon).

for the velocities, they're purely according to Oy with the module :

v = 2*pi*a*sqrt(1-e²)/(T*(1+e))

(with T the period)

of course, I have added the Earth velocity for the moon..


but my result, seems to be incorrect, since I find a minimal distance earth-sun of only 149.5 millions of km, whereas it should be around 147...

You can see my results and the graphics here : http://nicolas.aunai.free.fr/cours/magistere/3c/tl2a2s/tl2a.htm

(please do not pay attention to the x and y labels, I made some mistakes with my plotting program)

thanks if you can help me to solve this mystery :)

 

[pervect]

The very first thing I would suggest is cutting your timestep in half, and seeing if the results change. I dont' know if you tried this or not - I looked at your webpage, but I don't speak French ...

 

[Heimdall]

thank you for your answer. In fact my problem was in the starting conditions, it seems that the velocity I gave to the Earth was a little bit too high...

now the results are better but I still don't understands some things.

for exemple, why does the perigee of the moon is oscillating with an approximatly 400 days period (I checked this number on the web, it seems to be correct since I found 412 days exactly).

Then I don't understand why my Earth's Perihelion is 400000km too big when the moon is orbiting in direct direction (I find around 147.5 millions of km) whereas it seems to be perfectly correct when the moon is orbiting in oposite direction (I find here around the good value 147.1 millions of km)

you can see the plots for a 2 years simulation in direct direction http://nicolas.aunai.free.fr/cours/magistere/3c/2adirect.png

and the same plots for the oposite motion : http://nicolas.aunai.free.fr/cours/magistere/3c/2aretro.png

 

[pervect]

I'm guesssing a little bit, but I think the period of the oscillation is the sum of the lunar orbital period plus the Earth-sun orbital period.

 

[Heimdall]

I'm guesssing a little bit, but I think the period of the oscillation is the sum of the lunar orbital period plus the Earth-sun orbital period.

hum, I don't know why I have said that I had a period around 400 days... In fact the period of approximatly 200 days.

I've made a plot representing both the Earth-sun distance and the Earth moon distance, in a 10 years simulation, to see if there was an evident link between the two distances... apparently it is not as evident as I could expect... you can see it http://nicolas.aunai.free.fr/cours/magistere/3c/perigeeoscill.png

hum hum.. ?

 

[pervect]

OK, it looks like you have about 7 oscillations in 4 years, if I'm interpreting your diagram correctly. Which pretty much shoots down my idea.

At this point I don't have a solution, just some ideas that don't work. My first guess was that if we used the equivalent one dimensional problem, the moon would have a characteristic "resonant frequency"

I'm not sure if you're familair with the equivalent one-dimensonal problem, it's just a way of formulating the equation for the radial motion of an object

$$r'' = -GM/r^2 + r \dot{\theta}^2 = -GM/r^2 + (L/m)^2/r^3$$

Here L is the angular momentum of the orbiting object, we've made use of the constancy of L to remove the rate of change of theta from the equation and to make it an equation only in 'r'.

You can actually envision this equation as the 1 dimensional motion of an object in a 1d potential well, (the potential is the integral of the right hand side of the above equation) i.e. r'' = -dV(r)/dr

However, this quickly lead to the result that for a nearly circular orbit, the period of oscillation of 'r' would be one orbital period, i.e. one month for the moon.

Next, I looked at the graphs, and it reminded me of the output of an electrical circuit, a "mixer". I began to think that since a linear analysis dindn't seem to be solving the problem, it might be due to nonlinear effects. A linear circuit that is fed sine waves will produce only sine waves at the same frequency. A non-linear circuit that is fed sine waves will produce sine waves at the original frequencies and also generate output at the sum and difference frequencies.

However, there doesn't seem to be any particular simple relation to the output frequency and the two input frequencies.

I'll think about this a little more, but right now I don't have any great ideas. I do think that the right general approach is to model 'r' of the moon as a spring mass-system, possibly with a non-linear spring, using the 1 dimensional equivalent problem. The influence of the sun should then be a periodic pertubation of this, much as if we were applying an external force to this spring mass-system, or "shaking it" up and down.

 

[tony873004]

I think one problem you are running into is that it seems you are computing your x,y,z positions and velocities from the orbital elements. But since the Moon has mass it can not be considered a test particle. Its presence is continuously changing Earth's instantaneous orbital elements. This will introduce small errors. You can use JPL's Horizons page to get vectors directly for your program without having to start from orbital elements. These are the state vectors of Earth, Moon, and Sun for January 1, 2005 00:00:00 GMT. Assuming that the Sun is not fixed in your simulation, 0,0,0,0,0,0 will only be its position and velocity at T=0.

x,y,z (kilometers)
x-dot, y-dot, z-dot (kilometers / second)

Sun:
0,0,0
0,0,0

Earth (relative to Sun):
-27067884.6257174
144587718.527863
-1969.59490492195

-29.7609559402974
-5.58170729426966
7.82453662445448E-04

Moon (relative to Sun):
-27440663.8165208
144729193.335639
24808.236184314

-30.074128500212
-6.5095675566316
-0.057508367552749

Also, make sure you have good values for Earth, Moon and Sun mass. JPL Horizons gives:
Sun: 1.9891 x 10^30 kg
Earth: 5.9736 x 10^24
Moon: 734.9 x 10^20


As far as your time step, I would recommend using 1 second per step. Then in an effort to keep the formulas to a minimum inside the loop part of your program, don't adjust your time step there. Adjust it externally. For example, if you want a time step of x seconds per step, multiply the velocities of Sun, Earth, and Moon by x, and multiply their masses by x^2. Using this method, you avoid doing timestep math in your loop. You can keep your original values is separate variables if you wish, so they remain unaltered.