# B Solar system simulator JAVA source code

1. Mar 28, 2017

### 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

#### Attached Files:

• ###### src.zip
File size:
2.6 KB
Views:
87
2. Mar 28, 2017

### BvU

Hello sefiro,

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 ?

3. Mar 28, 2017

### nasu

4. Mar 28, 2017

### DrStupid

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.

5. Mar 28, 2017

### 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

Last edited: Mar 28, 2017
6. Mar 28, 2017

### DrStupid

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$

$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.

7. Mar 30, 2017

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

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

8. Mar 30, 2017

### BvU

You want to stick to the first derivative:

9. Mar 30, 2017

### sefiroths

10. Mar 30, 2017

### BvU

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

11. Mar 30, 2017

### BvU

I'm not fluent in Java. Could you elaborate on
-- 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 ?

 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: Mar 30, 2017
12. Mar 30, 2017

### DrStupid

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 (Text):

public Planet( /* insert your parameters 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 (Text):

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 (Text):

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

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

Code (Text):

public SolarSystem(){

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

Last edited: Mar 30, 2017
13. Mar 31, 2017

### 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

14. Mar 31, 2017

### BvU

Fully agree. So, probably no numerical problems.
And that looks excellent also. Somewhat puzzling that in post #1 you remark
which gave the impression you think there's something wrong ?

NIce check: how many moons in a year ?

15. Mar 31, 2017

### 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

16. Mar 31, 2017

### DrStupid

Increase the step width to 10000 seconds.

17. Apr 1, 2017

### 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

18. Apr 1, 2017

### DrStupid

Maybe you are interested in my javascript version.

19. Apr 4, 2017

### sefiroths

Thanks a lot, I'll check out