# Simple particles simulation (gravity)

1. Jul 21, 2013

### qwe

i believe i'm using the correct equations, but the orbit won't work. here is the relevant code for a 2 particle simple orbit simulation. i'm not sure what i'm doing wrong

Code (Text):

...

timeStep = 1.0
gravityConstant = 6.67384 * 10.0 ^ -11.0

...

// (initialize particles)
// p[n].m: mass, p[n].x: initial x position, p[n].vz: initial velocity vector along z axis

p[0].m = 1.0
p[1].m = 1.0

p[0].t = 0
p[0].x = 10.0

p[1].t = 0
p[1].x = -10.0

// initial velocity = sqrt(m2 ^2 * G / m1 + m2 * distance)

p[0].vz = sqrt( ( p[1].m * p[1].m * gravityConstant ) / ( p[0].m + p[1].m * 10.0 ) )
p[1].vz = sqrt( ( p[0].m * p[0].m * gravityConstant ) / ( p[1].m + p[0].m * 10.0 ) ) * -1

...

// (each loop, for each particle... p[n] is current particle, p[m] is other particle)

// calculate distance
dx# = p[m].x - p[n].x
dy# = p[m].y - p[n].y
dz# = p[m].z - p[n].z
distance# = sqrt( dx# * dx# + dy# * dy# + dz# * dz# )
// calculate force
force# = (gravityConstant * p[n].m * p[m].m) / (distance# * distance#)
// calculate acceleration, include timestep
acceleration# = (force# / p[n].m) * timeStep
// calculate unit vector and calculate magnitude of vector (acceleration)
p[n].vx = p[n].vx + dx# * (acceleration# / distance# )
p[n].vy = p[n].vy + dy# * (acceleration# / distance# )
p[n].vz = p[n].vz + dz# * (acceleration# / distance# )

Last edited: Jul 21, 2013
2. Jul 21, 2013

### Delta Kilo

Without going through the code in details:
1. I don't see where you update particle position. I assume you just omitted it for brevity.
2. Gravitation force between 2 1kg masses at 10 meters is not exactly very strong.
3. With a time step of 1s your are going to have a lot of very small steps, running into precision and error accumulation issues
4. Euler method sux. Bite the bullet and code up 4th order Runge-Kutta, then you can increase your step size to a sensible value.

3. Jul 22, 2013

### qwe

am i doing the timestep correctly? the bigger the timestep, the more acceleration applied. that doesn't seem right, the orbit shape changes based on the timestep, but it shouldn't

4. Jul 23, 2013

### Khashishi

Numerical integration suffers problems with accuracy and stability. The simple method you are using is called Euler method and isn't good. As Delta said, you should try 4th order Runge-Kutta.