Java Solar system simulator JAVA source code

AI Thread Summary
The discussion revolves around developing a solar system simulator in Java for educational purposes. The primary issue addressed is the incorrect behavior of the moon, which is mistakenly treated as a planet orbiting the sun instead of revolving around Earth. Participants suggest that the problem lies in the numerical integration method used, specifically the Euler method, which can lead to energy loss and inaccuracies. Recommendations include switching to the Leapfrog algorithm for better energy conservation and accuracy, as well as refining the code to properly calculate the gravitational interactions between celestial bodies. The developer plans to implement these suggestions and share the source code, potentially creating a JavaScript version for interactive use.
sefiroths
Messages
9
Reaction score
0
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
 

Attachments

Technology news on Phys.org
Hello sefiro, :welcome:

sefiroths said:
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 ?
 
BvU said:
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.
 
  • Like
Likes BvU
@DrStupid
Yes this is an error and seem that resolve the problem!
I have watched the code many times without see that "+"!
many thanks
 
Last edited:
sefiroths said:
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.
 
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:

ak = ∑ai

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

vk+1=vk+akΔ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:

xk+1=xk+vk+1Δt

that leads to what you have written:

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

Do you think I can say something like that or I have made some bad mistake?
thanks
 
sefiroths said:
Do you think I can say something like that
You want to stick to the first derivative:
sefiroths said:
xk+1=xk+vk⋅Δt +ak⋅Δt2
 
  • #10
That's a lot to read ...
4.11 is OK, but then he embarks on the midpoint method ...
 
  • #11
I'm not fluent in Java. Could you elaborate on
sefiroths said:
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 103 km/s, does this mean you follow trajectories in steps of (e.g. earth) 150 km ?
Positions are in units of 107 km ?
 
Last edited:
  • #12
sefiroths said:
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:

Code:
    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:

Code:
    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:

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

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

Code:
    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

WSMWm6A.gif
 
Last edited:
  • Like
Likes BvU
  • #13
@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:
Schermata 2017-03-31 alle 10.37.11.png

Thanks
 
  • #14
sefiroths said:
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.
sefiroths said:
What I can see by my program is:
And that looks excellent also. Somewhat puzzling that in post #1 you remark
sefiroths said:
but when I add the moon...
which gave the impression you think there's something wrong ?

NIce check: how many moons in a year ?
 
  • #15
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
 
  • #16
sefiroths said:
How can I make my program go in this spyral way?

Increase the step width to 10000 seconds.
 
  • #17
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
 
  • #18
sefiroths said:
perhaps I'll make a javascript version so we can play with parameters online

Maybe you are interested in my javascript version.
 
  • Like
Likes BvU
  • #19
Thanks a lot, I'll check out
 

Similar threads

Back
Top