Solar system simulator JAVA source code

[sefiroths]

Hi, I'm developing a very simple solar system simulator for educational purpose.
I had taken from wikipedia mass, velocity of planet and applied universal gravitation formula and applied to all celestial corps. It works well, but when I add the moon... the moon is like a planet, rotating around the sun.
I don't know if it is the right place to ask... sorry for that
I have made it in java, I hope that is clear what I have made
thanks

 

[BvU]

Hello sefiro,

the moon is like a planet, rotating around the sun

I think I know what you mean. Does the program ?

i.e.: How does the program know Luna's trajectory is mainly affected by Terra ?

 

[nasu]

The moon does revolve around the sun. If you look at the trajectory form a reference attached to the Sun, the trajectory looks like a large circle made from smaller arcs. You can get an idea here, for example:
https://phet.colorado.edu/sims/html/gravity-and-orbits/latest/gravity-and-orbits_en.html

 

[DrStupid]

i.e.: How does the program know Luna's trajectory is mainly affected by Terra ?

The calculation of the accelerations in lines 69-87 includes all interactions between all bodies according to Newton's law of gravitation.

But I don't understand why the step-width is increased with each step in line 145 (s.timestep+=0.005;). And I'm not sure if the integrator is correct.

 

[sefiroths]

@DrStupid
Yes this is an error and seem that resolve the problem!
I have watched the code many times without see that "+"!
many thanks

 

[DrStupid]

Yes this is an error and seem that resolve the problem!

Have you also checked the integrator? The position step of the Euler method looks like this:

$x_{k + 1} = x_k + v_k \cdot \Delta t + {\textstyle{1 \over 2}}a_k \cdot \Delta t^2$

but your code results in

$x_{k + 1} = x_k + v_k \cdot \Delta t + a_k \cdot \Delta t^2$

And even the correct Euler method is not suitable for this problem. I would suggest the Leapfrog algorithm because it conservs energy but is still quite simple. It constists in the Euler step for the position and

$v_{k + 1} = v_k + {\textstyle{1 \over 2}}\left( {a_k + a_{k + 1} } \right) \cdot \Delta t$

for the velocity.

 

[sefiroths]

I don't know if I can make something like I've written. I wanted to use this example in a high school and they don't know what is an integration, derivate or Euler method... so I thought something like this:
1) resultant a at istant k is:

a_k = ∑a_i

given by attraction of every planet
2) now that I have the a_k I can calculate the new velocity of the planet:

v_k+1=v_k+a_kΔt

3) now that I have the velocity, for a Δt small I can consider the velocity of the planet constant(?) so I use the Rectilinear motion formula for the position:

x_k+1=x_k+v_k+1Δt

that leads to what you have written:

x_k+1=x_k+v_k⋅Δt+a_k⋅Δt²

Do you think I can say something like that or I have made some bad mistake?
thanks

 

[BvU]

Do you think I can say something like that

You want to stick to the first derivative:

xk+1=xk+vk⋅Δt +ak⋅Δt2

 

[sefiroths]

I have searched something and found this:
https://kof.zcu.cz/st/dis/schwarzmeier/gravitational_simulation.html
see (4.11), seem what I have implemented... is it correct?

 

[BvU]

That's a lot to read ...
4.11 is OK, but then he embarks on the midpoint method ...

 

[BvU]

I'm not fluent in Java. Could you elaborate on

It works well, but when I add the moon ...

-- which makes me suspicious of numerical problems: how long does it 'work well' ? What are the time steps ? Do $\vec a$(Luna) and $\vec v$(Luna) come out as expected ?

[edit] found out this 5 ms time step.
With speeds in terms of 10³ km/s, does this mean you follow trajectories in steps of (e.g. earth) 150 km ?
Positions are in units of 10⁷ km ?

 

[DrStupid]

I have searched something and found this:
https://kof.zcu.cz/st/dis/schwarzmeier/gravitational_simulation.html
see (4.11), seem what I have implemented... is it correct?

Yes, that seems to be the Euler method for a system of 1st order ODEs. The error of this methods results in a loss of energy. What I referred to as Euler method (I'm not completely sure if this is correct) results in an increasing energy. The truth is somewhere in the middle. That's why your souce introduces the midpoint method. It is much better, but needs two accelerations per step. The Leafrog methos needs only one acceleration per step and is therefore as fast as Euler but much more accurate.

In order to implement Leapfrog into your code you need additional variables to keep the acceleration. The best place would be the Planet object:

    public Planet( /* insert your parameters here */ ){
  
       // insert your code here
  
       this.accX = 0;
       this.accY = 0;
       this.old_accX = 0;
       this.old_accY = 0;
   }

In order to change the order for the calculation of positions, velocities and acceleration it is a good idea to write separate procedures:

    void move(){
       for(Planet p1 : lista){
  
           // update position according to dr = v·dt + a·dt²/2
      
           p1.x += timestep * (p1.velX + timestep * p1.accX / 2);
           p1.y += timestep * (p1.velY + timestep * p1.accY / 2);
      
           // keep accelerations for velocity calculation
      
           p1.old_accX = p1.accX;
           p1.old_accY = p1.accY;
       }
   }

   void getAccelerations(){
       double dx, dy, dz, D, A;
       for(Planet p1 : lista){
           p1.accX = 0;
           p1.accY = 0;
           for(Planet p2 : lista){
               if (p1 != p2){
                   dx = p2.x - p1.x;
                   dy = p2.y - p1.y;
                   D = Math.sqrt(Math.pow(dx, 2) + Math.pow(dy, 2));
                   A = G * p2.mass / Math.pow(D, 2);
                   p1.accX += dx * A / D;
                   p1.accY += dy * A / D;
               }      
           }
       }
   }

   void accelerate(){
       for(Planet p1 : lista){

           // update velocity according to dv = (a + old_a)·dt/2
      
           p1.velX += (p1.accX + p1.old_accX) * timestep / 2;
           p1.velY += (p1.accY + p1.old_accY) * timestep / 2;
       }
   }

Now the update of the positions looks like this:

    void updatePositions(){
       move();
       getAccelerations();
       accelerate();
   }

As the accelerations are required at the begin of each step, there must be a corresponding initialisation:

    public SolarSystem(){

       // insert your code here

       getAccelerations();
   }

I didn't check if you have direct access to the lista array with an index in Java. If this is possible than you can easily double the speed of the getAccelerations procedure by use of Newton's third law.

PS: Here is a comparison of the methods:

red: $x_{k + 1} = x_k + v_k \cdot \Delta t + a_k \cdot \Delta t^2$ (your code)
blue: $x_{k + 1} = x_k + v_k \cdot \Delta t + {\textstyle{1 \over 2}}a_k \cdot \Delta t^2$
green: Leapfrog

 

[sefiroths]

@BvU
I have left the program running all night, also changed timestep with 10, 1, 0.5... seems that speed corresponds (scaled by my program
//1 sdu = 10^10 m // scaled distance unit
//1 smu = 10^29 kg //scaled mass unit
//1 stu = 10^5 s
)
In next version I'll use coords, vels, ecc... without any scaling factor
However if sonething went wrong I think I had to see a planet falling in the sun or go away from the sun...
@DrStupid
This sound very interesting, I'll try to implement all equations to see the comparison.
How can I make my program go in this spyral way?
What I can see by my program is:

Thanks

 

[BvU]

However if sonething went wrong I think I had to see a planet falling in the sun or go away from the sun...

Fully agree. So, probably no numerical problems.

What I can see by my program is:

And that looks excellent also. Somewhat puzzling that in post #1 you remark

but when I add the moon...

which gave the impression you think there's something wrong ?

NIce check: how many moons in a year ?

 

[sefiroths]

In the first post I had a problem that the moon had an orbit like the earth, on post #5 i wrote that I had solved the problem, then DrStupid point my attention on a potentially wrong calculation that should lead a falling planet over the sun.
Now I'm searching how to enphatize this error seeing the planet falling down, but after hour of iterations the planet stay on his orbit...
For the check I'll change the program to calculate it

 

[DrStupid]

How can I make my program go in this spyral way?

Increase the step width to 10000 seconds.

 

[sefiroths]

tried and this corrupt even other simulations with other algorithm. in next version I'll implement all methods of calculation to make it comparable
Thanks a lot, I'll provide all source code after all, and perhaps I'll make a javascript version so we can play with parameters online

 

[DrStupid]

perhaps I'll make a javascript version so we can play with parameters online

Maybe you are interested in my javascript version.

 

[sefiroths]

Thanks a lot, I'll check out