Simulating Motion under Gravitation

Even 2 body problems can lead to this

Let's be honest! I am falling in love with the p5.js framework. It makes my algorithms and simulations so beautifully presentable! And that too with so minimum amount of boilerplate.

After building a Mandelbrot set simulator, I thought, why to stop here, let's build something related to Physics, which I loved very much during my +2 days. Then the idea of creating the Gravitation simulator came into my mind.

The Physics Involved

Gravitation comes under (classically) Central forces.

The magnitude of gravitational force between 2 point masses $m_{1}$ and $m_{2}$ separated by a displacement of magnitude $r$ is given by:

$F = G \cdot \frac{m_{1} \cdot m_{2}}{r^{2}}$

where $G$ is the universal gravitational constant, with value $G = 6.674 \cdot 10^{- 11} m^{3} k g^{- 1} s^{- 2}$ .

Combining this with Newton's famous 2nd law, the magnitude of acceleration $a$ of $m_{1}$ is given by:

$a = \frac{F}{m_{1}} = G \cdot \frac{m_{2}}{r^{2}}$

Now, actually in space, every celestial body moves towards every other due to the Gravitation. That'd indeed be an incredible multi-body dynamics problem to solve. But our scope is pretty simple. We'll just consider $m_{1}$ to be movable. Other bodies would be stationary.

Gravitation is always attractive (Shhhh! GTR fans). What that means is, if, by some magic, $m_{2}$ was pinned down at a fixed point, $m_{1}$ will always move TOWARDS $m_{2}$ with the acceleration given above.

Note that the acceleration is independent of $m_{1}$ , so whether we take a rock or a planet, the path traced out by them would be the same.

Approximations and rescalings

The Gravitation is very weak as a force. Using $G$ directly in the simulation is not practical. So for the simulation, I use a different scaling constant, which I pick by hand and constrain the acceleration values within some handpicked bounds, to do away with the problem of division by zero, when $r = 0$ .

Another cheating that I do in the simulation, is that I keep all object's mass as 1. For the simulator, $1$ body of mass $m$ , is the same as, $m$ bodies of mass $1$ located at the same point.

I also work in 2D plane and decompose the accelerations in X and Y axes using the slope of the line joining 2 bodies.

So, if $m_{1}$ (moving mass) is located at $(x_{1}, y_{1})$ and $m_{2}$ (fixed mass) is located at $(x_{2}, y_{2})$ , the components used are (with a scaling factor of 100):

$a_{x} = \frac{100 \cdot (x_{2} - x_{1})}{r^{\frac{3}{2}}}$

$a_{y} = \frac{100 \cdot (y_{2} - y_{1})}{r^{\frac{3}{2}}}$

I also clip $a_{x}$ and $a_{y}$ between some bounds.

The total acceleration of the mass $m_{1}$ is the component-wise sum of accelerations over all the fixed masses. That is to say, at a given time:

${\vec{a}}_{t o t a l} = Σ a_{x} \hat{i} + Σ a_{y} \hat{j}$

I'll drop the sum sign from here on.

The last approximation has to do with how I calculate the velocity and position in the next timestep.

In reality, velocity would be a nice and smooth time integral of the acceleration and the displacement would be a time integral of the velocity.

Here I approximate integral by sum.

I take a small timestep $Δ t$ (in simulation, it has value 1). Then calculate:

${\vec{v}}_{n e x t} = {\vec{v}}_{0} + \vec{a} \cdot Δ t$

${\vec{r}}_{n e x t} = {\vec{r}}_{0} + {\vec{v}}_{n e x t} \cdot Δ t$

Of course, the calculations are done component-wise.

The code, finally

With all the above considerations in mind, let's jump right into the code.

	<!DOCTYPE html>
	<html>
	<head>
	<title>Gravity</title>
	<meta name="charset" content="utf-8">
	<meta name="viewport" content="width=device-width">
	<script src="https://cdnjs.cloudflare.com/ajax/libs/p5.js/0.10.2/p5.js"></script>
	<script src="https://cdnjs.cloudflare.com/ajax/libs/p5.js/0.10.2/addons/p5.sound.min.js"></script>
	</head>

	<body style="display: flex; flex-direction: row;">
	<div style="display:flex; flex-direction:column;">
	<input type="number" max="600" min="0" id="posX" placeholder="posX">
	<input type="number" max="600" min="0" id="posY" placeholder="posY">
	<input type="number" max="1000" min="-1000" id="velX" placeholder="velX">
	<input type="number" max="1000" min="-1000" id="velY" placeholder="velY">
	<button onclick="start()">Start</button>
	<button onclick="reset()">Reset</button>
	</div>
	<script src="sketch.js"></script>
	</body>
	</html>

view raw index.html hosted with ❤ by GitHub

The index.html loads the p5.js library as usual. Notice that I have created fields for initial position and velocity as well as buttons for starting and resetting the simulation.

The main juice of the code is sketch.js:

	let particles = [];
	let position = [100, 100];
	let velocity = [0, 0];
	let started = false;
	let stuck = false;
	let lines = [];
	function start() {
	started = true;
	}
	function reset(){
	started = false;
	stuck = false;
	particles = [];
	lines = [];
	}

	function setup(){
	createCanvas(600, 600);
	}

	const timestep = 1;

	function draw(){
	let mindist = 100000000;
	background(200);
	if (!started){
	if (mouseIsPressed){
	if (mouseX >= 0 && mouseX <= 600)
	particles.push([mouseX, mouseY]);
	console.log(particles)
	}
	position[0] = Number(document.getElementById("posX").value);
	position[1] = Number(document.getElementById("posY").value);
	velocity[0] = Number(document.getElementById("velX").value);
	velocity[1] = Number(document.getElementById("velY").value);
	}else{
	let acc = [0, 0];
	particles.forEach(elt => {
	let dist_sqr = (elt[0] - position[0]) * (elt[0] - position[0]) + (elt[1] - position[1]) * (elt[1] - position[1]);
	if (dist_sqr / (velocity[0] * velocity[0] + velocity[1] * velocity[1]) < 0.1) stuck = true;
	if (dist_sqr < mindist) mindist = dist_sqr;
	acc[0] += constrain(10 * (elt[0] - position[0]) / (dist_sqr * Math.sqrt(dist_sqr)), -10, 10);
	acc[1] += constrain(10 * (elt[1] - position[1]) / (dist_sqr * Math.sqrt(dist_sqr)), -10, 10);
	});

	if (!stuck){
	velocity[0] += acc[0] * timestep;
	velocity[1] += acc[1] * timestep;
	position[0] += velocity[0] * timestep;
	position[1] += velocity[1] * timestep;
	lines.push([...position]);
	}
	}

	particles.forEach(elt => {
	fill(100, 200, 100);
	ellipse(elt[0], elt[1], 10, 10);
	});

	for (let i = 0; i < lines.length - 1; i++){
	stroke(255, 0, 0);
	line(lines[i][0], lines[i][1], lines[i+1][0], lines[i+1][1]);
	}

	stroke(0, 0, 0);
	fill(255, 0, 0);
	rect(position[0], position[1], 5, 5);
	}

view raw sketch.js hosted with ❤ by GitHub

particles array contains the location of the fixed unit masses;
position and velocity correspond to the moving mass.
started and stuck are simulation states.
lines contains the set of points along the path traced by the moving mass.

In the draw function, I check if the mouse is pressed WITHIN THE CANVAS, I add that location to the particles array.

However if started is set to true, I no longer add bodies, rather I proceed on calculating the next position at every frame.

Finally, I draw the fixed particles as green circles and the moving particle as a red body which leaves a red trace.

This logic works fine, until the red body from very far off comes very near a green body. In this process it gains tremendous velocity and shoots out of the screen. To prevent that from happening, the stuck variable is set to true whenever such stuff happens. And once set to true, the particle's position is not updated, it remains stuck.

This sticking behaviour is controlled by the metric in Line 40. The value to which it is set in the code hardly has any effect. Change the threshold to 1 or 10 to see the sticking effect.