I'm getting strange results when calculating the physics between two objects.

My goal is to build a galaxy simulator that can simulate how galaxies consisting of 100s of stars might look when interacting. It's a fun little side project and I'm not even sure how it'll turn out (which is part of the fun). I have a really limited background in math and physics so this is a pretty big challenge for me atm. I'm starting with one star and a static object to attract it to nail down the basic physics.

The problem I'm facing is that, as star one approaches the static object it's is accelerating as expected and when it passes the object I expect the velocity to slowly bleed to zero and the star to fall back in. However, after passing the static object the star is behaving oddly. Depending on my delta time variable It position sometimes shoots off crazy fast in one direction or another. I'm too new to all of this and there are too many variables for me to follow where exactly the issue lies.

As an FYI this is a little web app written in javascript and not an actual video game. However, the concepts are similar to calculating video game physics and most web developer forums aren't helpful. Here is a link to the code in codepen. NOTE: If the “try me” button doesn't appear try shrinking your browser down in size or scrolling down to find the button.

Finally, the rendering engine (canvas) doesn't appear to be the problem as I can see the issue in the console as being related to the calculations and not he rendering.

let frameCount = 0;
let stars;
let dt;
let gravConst = 6.674 * Math.pow(10, -11);
// let solarMass = Math.pow(10, 30);
// let lightYear = 9.5 * Math.pow(10, 15);
let solarMass = 1;
let lightYear = 1;
let yearCounter;

stars = [
 { x: 250, y: 75, m: 1, v: 0 },
 { x: 275, y: 75, m: 1, v: 0 }
];
dt = 50000;
yearCounter = 0;

const render = () => {
let canvas = document.getElementById("canvasRef");
const context = canvas.getContext("2d");
let animationFrameId;
frameCount++;
animationFrameId = window.requestAnimationFrame(render);
draw(context, frameCount);
};

const draw = (ctx, frameCount) => {
// console.log(frameCount)
ctx.clearRect(0, 0, ctx.canvas.width, ctx.canvas.height);
ctx.fillStyle = "#ffffff";
calculateGravAccel(stars[0], stars);
drawStar(stars, ctx);
};

let calculateGravAccel = (star, stars) => {
// Force on current star by all other star sources (newtons)
let m1 = stars[0].m * solarMass;
let m2 = stars[1].m * solarMass;
let distance = (stars[0].x - stars[1].x) * lightYear;
let xDistance = stars[0].x - stars[1].x;
let f = gravConst * (m1 * m2) / Math.pow(distance, 2);
console.log(f);
// Convert velocity to account for the fact that the force is now coming from behind
// New velocity of current star
star.v = distance < 0 ? star.v + (f / m1) * dt : star.v - (f / m1) * dt;
// Position
star.x = star.x + star.v * dt;
// Ellapsed Time
yearCounter += dt;
};

const drawStar = (stars, ctx) => {
console.log(`x:${stars[0].x} y: ${stars[0].y} v: ${stars[0].v}`);
ctx.beginPath();
ctx.arc(stars[0].x, stars[0].y, 5, 0, 2 * Math.PI);
ctx.arc(stars[1].x, stars[1].y, 5, 0, 2 * Math.PI);
ctx.fill();
};

Do you know C++?

@taby , Not really. Like I said I'm a web developer and this is sharpening my javascript skills. However, the problem itself is very game dev like.

OK, well it's not super clear where you're going wrong.

The code below should be easy enough to decipher.

Here is my acceleration function:

custom_math::vector_3 acceleration(const custom_math::vector_3 &pos, const custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
	custom_math::vector_3 grav_dir = sun_pos - pos;

	float distance = grav_dir.length();

	grav_dir.normalize();
	custom_math::vector_3 accel = grav_dir * (G*sun_mass / pow(distance, 2.0));

	return accel;
}

Some integrators are:

// https://en.wikipedia.org/wiki/Symplectic_integrator
// Also known as Verlet integration
void proceed_symplectic2(custom_math::vector_3 &pos, custom_math::vector_3 &vel)
{
	static const double c[2] = { 0, 1 };
	static const double d[2] = { 0.5, 0.5 };

// 	pos += vel * c[0] * dt; // first element c[0] is always 0
	vel += acceleration(pos, vel, ang_vel) * d[0] * dt;

	pos += vel * c[1] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[1] * dt;
}

// https://www.gamedev.net/forums/topic/701376-weird-circular-orbit-problem/?do=findComment&comment=5402054
// https://en.wikipedia.org/wiki/Symplectic_integrator
void proceed_symplectic4(custom_math::vector_3 &pos, custom_math::vector_3 &vel)
{
	static double const cr2 = pow(2.0, 1.0 / 3.0);

	static const double c[4] =
	{
		1.0 / (2.0*(2.0 - cr2)),
		(1.0 - cr2) / (2.0*(2.0 - cr2)),
		(1.0 - cr2) / (2.0*(2.0 - cr2)),
		1.0 / (2.0*(2.0 - cr2))
	};

	static const double d[4] =
	{
		1.0 / (2.0 - cr2),
		-cr2 / (2.0 - cr2),
		1.0 / (2.0 - cr2),
		0.0
	};

	pos += vel * c[0] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[0] * dt;

	pos += vel * c[1] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[1] * dt;

	pos += vel * c[2] * dt;
	vel += acceleration(pos, vel, ang_vel) * d[2] * dt;

	pos += vel * c[3] * dt;
//	vel += acceleration(pos, vel, ang_vel) * d[3] * dt; // last element d[3] is always 0
}

void proceed_Euler(custom_math::vector_3 &pos, custom_math::vector_3 &vel)
{
	vel += acceleration(pos, vel, ang_vel) * dt;
	pos += vel * dt;
}

inline void proceed_RK2(custom_math::vector_3 &pos, custom_math::vector_3 &vel)
{
	custom_math::vector_3 k1_velocity = vel;
	custom_math::vector_3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
	custom_math::vector_3 k2_velocity = vel + k1_acceleration * dt*0.5;
	custom_math::vector_3 k2_acceleration = acceleration(pos + k1_velocity * dt*0.5, k2_velocity, ang_vel);

	vel += k2_acceleration * dt;
	pos += k2_velocity * dt;
}

void proceed_RK4(custom_math::vector_3 &pos, custom_math::vector_3 &vel)
{
	static const double one_sixth = 1.0 / 6.0;

	custom_math::vector_3 k1_velocity = vel;
	custom_math::vector_3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
	custom_math::vector_3 k2_velocity = vel + k1_acceleration * dt*0.5;
	custom_math::vector_3 k2_acceleration = acceleration(pos + k1_velocity * dt*0.5, k2_velocity, ang_vel);
	custom_math::vector_3 k3_velocity = vel + k2_acceleration * dt*0.5;
	custom_math::vector_3 k3_acceleration = acceleration(pos + k2_velocity * dt*0.5, k3_velocity, ang_vel);
	custom_math::vector_3 k4_velocity = vel + k3_acceleration * dt;
	custom_math::vector_3 k4_acceleration = acceleration(pos + k3_velocity * dt, k4_velocity, ang_vel);

	vel += (k1_acceleration + (k2_acceleration + k3_acceleration)*2.0 + k4_acceleration)*one_sixth*dt;
	pos += (k1_velocity + (k2_velocity + k3_velocity)*2.0 + k4_velocity)*one_sixth*dt;
}

You're probably familiar with the proceed_Euler function, which is what you're kind of doing now. I would recommend that you ditch the Euler integrator, and use the symplectic or Runge-Kutta integrators intead – if you don't, then you'll find that the orbits don't form properly.

@taby Awesome response. I really appreciate the effort!

@taby Also, how did you get your codeblocks to highlight function names? edit: nevermind it took a second for mine to do that as well.

EDIT: Posted code has wrong math - ignore!

I did such planets simulation for some fun, but i had no static bodies.
So there is no gravity constant, instead each planet affects all others.
I've left integration to the physics engine (RK4), so given code only calculates the force for each body:

	void StepPlanets ()
	{
		for (int i=0; i<planets.size(); i++)
		{
			Body *bI = planets[i];
			float massI, massJ;
			sVec3 comI, comJ;
			comI = BodyGetCenterOfMass (bI);
			massI = BodyGetMass (bI);

			for (int j=0; j<planets.size(); j++) if (j != i)
			{
				Body *bJ = planets[j];
				comJ = BodyGetCenterOfMass (bJ);
				massJ = BodyGetMass (bJ);

				sVec3 d = comJ - comI;
				sVec3 f = d * (massI + massJ) / d.SquaredLength();

				((BodyData*)BodyGetUserData(bI))->force += f;
				((BodyData*)BodyGetUserData(bJ))->force -= f;
			}
		}
	}

@JoeJ The gravitations constant, as far as I can tell, is necessary in determining the force applied to each body by another body. So in your sVec3 f = d * (massI + massJ) / d.SquaredLength(); shouldn't the first d actually be swapped out with the gravitational constant?

It's worth noting that there are cases where G = 1, like when using geometrized units.

You could also limit the speed and acceleration. For a photon, using the reference (e.g. Euler) integrator:

void proceed_Euler(custom_math::vector_3 &pos, custom_math::vector_3 &vel)
{
	custom_math::vector_3 accel = acceleration(pos, vel);

	vel += accel * dt;

	if (vel.length() > speed_of_light)
	{
		vel.normalize();
		vel *= speed_of_light;
	}

	pos += vel * dt;
}

To limit the acceleration, perform the following directly after the gravitational acceleration is calculated:

if (total_accel.length() > speed_of_light*2)
{
	total_accel.normalize();
	total_accel *= speed_of_light*2;
}

You could also use a Lennard-Jones potential, where the interaction is attractive at large distances, and repulsive close up. Have fun!!

The problem I'm facing is that, as star one approaches the static object it's is accelerating as expected and when it passes the object I expect the velocity to slowly bleed to zero and the star to fall back in. However, after passing the static object the star is behaving oddly. Depending on my delta time variable It position sometimes shoots off crazy fast in one direction or another. I'm too new to all of this and there are too many variables for me to follow where exactly the issue lies.

Getting very high forces when bodies are close is normal. When the distance between two bodies is zero, the inverse square function pushes the values towards infinite ( 1 / 0 == lots ).

The bodies won't slow down when they get close just from the gravity. You need to do something when they get close …. you can justify this however you like ( ie. drag from gasses when they are near, or a collision ) – but you will need to do something to get them to slow down, and fall inwards.

When I was doing this kind of thing, I had circles colliding, so the bodies never got close enough to go nuts. Also when I applied drag, they would eventually slow down enough to fall into each other:

drag_force = normalize( velocity ) * pow( length( velocity ), 2.0 ) * drag_constant * -1.0;

From you description, there is nothing wrong with what you are doing, it is just how the system behaves!

It is a good idea to clamp the lower limit of distance, so that the number never gets huge:

distance = max( length(p1 - p2), 0.01 );

force = normalize( p1 - p2 ) * ( 1.0 / pow( distance, 2.0 ) );

Add some other conditions to make it do what you want before it ever gets to this, but it is a good safety. In reality something else would happen in these scenarios – like the body would break apart if it ever had huge forces.

Also have you seen this? (227) Illustris Simulation: Most detailed simulation of our Universe - YouTube