gravity

I've just finished writing all my US college essays!

This is a massive accomplishment. I should be relieved. I am.

What to do now? Study for my prelims that are coming up? Nah. It's time to write a blog post.

No, but, in all seriousness. I've been wanting to write here for a while, but I never got the time to. I've had to write a billion essays for my US apps, but they made me realise something important: I'm not that terrible at writing.

Exposition aside, I want to write about something really cool that I'd explored in part some months ago: gravity simulations!

Wait, what? Gravity simulations? From you? Aren't you the pure maths guy who thinks physics is too practical to be enjoyable?

That's true. I don't know how I ended up here either, but I somehow did. I think it happened at midnight. I was learning about magnetism for my school physics exams, and I started thinking:

Why don't magnetic monopoles exist? While that question alone probably deserves its own blog post, I'd asked this question to myself because in our syllabus, we learn electrostatics before magnetism. So I was already familiar with basic formulas like:

$$ F = \frac{1}{4\pi\epsilon_0} \frac{q_1 q_2}{r^2} $$

I noticed a stark contrast in difficulty while doing magnetism from when I was doing electrostatics: everything was just so much harder here: and I didn't understand why. Conceptually, they're the same, right? Just attraction and repulsion, just like before!

I didn't feel like studying much. This is probably when mass came to mind. I thought: wait a minute, screw magnetism, aren't charges very similar to bodies with mass? I mean, we have Newton's law of gravitation:

$$ F = G \frac{m_1 m_2}{r^2} $$

that looks virtually identical to charges. That is interesting.

Of course, there are obvious differences, like how charges are signed but mass isn't, or how electrostatic force can be repulsive but gravitational force can't. I won't delve much into this, because what I really want to write about is related more to gravity: specifically that formula.

I think I was just bored. I remember recently learning about HTML Canvas and decided to use my skills for good. I opened VSCode and quickly coded up a small program that let me spawn objects. It looked something like this (click to add dots):

That's cool. Now, it's time to spice things up. In order to add gravity, I'd need three extra attributes: mass, vx, and vy. So far I have color, x, y, and radius. Color and radius are purely cosmetic. On the other hand, x and y are needed for position, mass is needed to determine how much inertia an object has, and vx and vy are needed to store the velocity of the body, which we'll need for actually simulating motion.

Alright. Let's add those.

document.addEventListener("click", (e) => {
  b.push({mass: 5e30, vx: 0, vy: 0, color: '#ff0000', x: e.clientX, y: e.clientY, radius: 10});
});

At this point, you may have questions. Why is the mass 5e30? Why is this object so massive?

Well... good question. The truth is, I tried many masses before this. I eventually realised that gravity would only work for really big ones. Blame the gravitational constant, not me.

Either way, I'd decided that I wanted my simulation to be 'as accurate as possible' so, naturally, I felt obligated to use astronomical masses.

So that tracks.

What's next? Oh right, I need to add the actual... gravity itself. This should be simple. Just apply Newton's formula on all $O(n^2)$ pairs, choosing a small enough $\Delta t$, right?

Not so fast. Before const G = 6.6743e-11; would make any sense, I'd have to determine how I intended to convert between pixels and kilometers. You see, the radius attribute I included before was purely visual and in pixels: I didn't really care about the physical signficance of it. On the other hand, all the calculations I'll need to perform for gravity very much care about units: so I'd probably have to as well.

How do I decide a scale?

I know that I definitely want the earth and the sun to fit in one frame. How would I have fun otherwise? Okay, so I need to be able to represent 1.5e8.

I decided to go with a diagonal distance of 400 pixels corresponding to 1.5e8 kilometers. More formally:

$$ 400\sqrt{2} \text{ pixels} = 1.5 \times 10^8 \text{ km} $$ $$ \implies 1 \text{ pixel} \approx 265,165 \text{ km} $$

With this done, I can now code the nested for loop:

for (let i = 0; i < b.length; ++i) {
  for (let j = i + 1; j < b.length; ++j) {
    let r = 265165 * Math.sqrt((b[i].x - b[j].x) * (b[i].x - b[j].x) + (b[i].y - b[j].y) * (b[i].y - b[j].y));
    let F = G * (b[i].mass / r) * (b[j].mass / r);
    let F_vec = [265165 * (b[i].x - b[j].x), 265165 * (b[i].y - b[j].y)];
    F_vec[0] /= r, F_vec[1] /= r;
    F_vec[0] *= F, F_vec[1] *= F;
    let a_ix = F_vec[0] / b[i].mass, a_iy = F_vec[1] / b[i].mass, a_jx = F_vec[0] / b[j].mass, a_jy = F_vec[1] / b[j].mass;
    b[i].vx -= 1 / hz * a_ix / 265165, b[j].vx += 1 / hz * a_jx / 265165;
    b[i].vy -= 1 / hz * a_iy / 265165, b[j].vy += 1 / hz * a_jy / 265165;
  }
}

Okay. Okay. What just happened? What even is this code?!

Don't worry. It wasn't that easy coding it up either, and definitely didn't happen all at once. Let's go line by line.

let r = 265165 * Math.sqrt((b[i].x - b[j].x) * (b[i].x - b[j].x) + (b[i].y - b[j].y) * (b[i].y - b[j].y));

This makes sense. r is the distance between the two bodies. The Math.sqrt part is just Euclidean distance, and multiplying by 265165 converts pixels to kilometers.

let F = G * (b[i].mass / r) * (b[j].mass / r);

Alright, nice. This is just $F = G \frac{m_1 m_2}{r^2}$! So far, no issues whatsoever.

let F_vec = [265165 * (b[i].x - b[j].x), 265165 * (b[i].y - b[j].y)];

Hm, what? F_vec? Why do we need that? Oh. We're in two dimensions. Alright. I respect vectors, so this is cool, still very principled stuff. But wait, how is force equal to just $\Delta s$ in km?

Oh okay, it's mutated in the next lines:

F_vec[0] /= r, F_vec[1] /= r;
F_vec[0] *= F, F_vec[1] *= F;

It's obvious now: we initialise F_vec as a unit vector pointing from one object to the other, and then simply scale it up by F: the magnitude of force we calculated earlier.

let a_ix = F_vec[0] / b[i].mass, a_iy = F_vec[1] / b[i].mass, a_jx = F_vec[0] / b[j].mass, a_jy = F_vec[1] / b[j].mass;

This is also very simple: we're just inverting $F = ma$. i and j refer to the two objects, and x and y are the two dimensions.

b[i].vx -= 1 / hz * a_ix / 265165, b[j].vx += 1 / hz * a_jx / 265165;
b[i].vy -= 1 / hz * a_iy / 265165, b[j].vy += 1 / hz * a_jy / 265165;

While the opposite signs might be confusing at first, we're only calculating $F_{ji}$, and by Newton's third law, i will experience an equal and opposite force in the other direction. We divide by 265165 to convert back to pixels, and the 1/hz part ensures that we... dampen the force by a bit.

While I could've probably adjusted mass or distance or something else about the objects I chose, I found the dt option to be the easiest at the moment, at the cost of being less 'principled'.

Let's not forget to actually update the distances based on the velocity. Adding these two lines before drawing the object suffices.

b[i].x += b[i].vx;
b[i].y += b[i].vy;

After all that hard work, we finally get a somewhat decent simulation!

... Wow.

While I knew what was coming, using it for the first time felt miraculous. This was not due to gravity or anything, but due to the fact that I've tried to create smooth animations programmatically before, and they've never worked. However, all of a sudden, I have an animation of a circle moving towards another circle, and it looks like it was generated by a Bezier curve. How can such beautiful visuals emerge from such simple rules?

I first tried this, obviously.

You can recreate this: try placing two objects.

That was cool. But let's make this more interesting.

This one is a bit trickier. You can't just place three objects, you need to get one moving towards the left decently fast and then place one to the right of it and you'll see... no spoilers.

I asked and the simulation delivered.

Seeing that was sufficient to get me curious enough to Google this phenomenon. That is when I came across the Wikipedia article for the two body problem.

I read this line:

By contrast, the three-body problem (and, more generally, the n-body problem for n > 2) cannot be solved generally, except in special cases.

By contrast. This means the two-body problem... can be solved generally, in all cases. Okay, what the hell... I need to try.

I then rushed to grab a pen and some paper. I should remind you: it was still midnight.

Let's do a quick runthrough, because my lucidity didn't last for long.

Given $\vec{F}$, acceleration ($\vec{a}$) is known. We want position $\vec{r}$.

$$ \begin{aligned} \vec{v} &= \int \vec{a} \, dt \\ \vec{r} &= \iint \vec{a} \, dt \, dt \end{aligned} $$

The force on body A is:

$$ \vec{F}_A = G \frac{m_a m_b}{|\vec{r}_b - \vec{r}_a|^2} \frac{\vec{r}_b - \vec{r}_a}{|\vec{r}_b - \vec{r}_a|} $$

So the position of A at time $t$ is:

$$ \vec{r}_a(t) = \frac{G m_a m_b}{m_a} \iint \frac{\vec{r}_b(t) - \vec{r}_a(t)}{|\vec{r}_b(t) - \vec{r}_a(t)|^3} \,dt\,dt $$ $$ \vec{r}_b(t) = \frac{G m_a m_b}{m_b} \iint \frac{\vec{r}_a(t) - \vec{r}_b(t)}{|\vec{r}_a(t) - \vec{r}_b(t)|^3} \,dt\,dt $$

Or, in more standard form, we can differentiate both equations twice and get:

$$ \frac{d^2 \vec{r}_a}{dt^2} = \frac{G m_a m_b}{m_a} \frac{\vec{r}_b - \vec{r}_a}{|\vec{r}_b - \vec{r}_a|^3} $$ $$ \frac{d^2 \vec{r}_b}{dt^2} = \frac{G m_a m_b}{m_b} \frac{\vec{r}_a - \vec{r}_b}{|\vec{r}_a - \vec{r}_b|^3} $$

I tried, but I wasn't able to solve this. If anyone knows how to, please let me know: I'd appreciate it.

Any way, that's not the last thing I did. Next, I realised: I have a gravity simulator (*) anyway, I have to try some orbits.

After a quick Google search for the right initial conditions for a stable three body configuration, I made this.

(Click to spawn extra mass. Click canvas then use Arrow Keys to move bodies)

I also tweaked hz (time is now a lot slower) and the masses a bit (to be more precise, I made the three initial bodies really heavy compared to the spawned ones), so you'll notice changed behaviour on object spawning. But we now have a stable three body configuration!

I've had my fun for now. If I were to continue this, it'd almost certainly be about the two differnetial equations from before, because I really want to know how I'd solve them.

Some disclaimers:

If you were expecting this blog to head in a concrete direction or with a satisfying conclusion, sorry for not providing that. The truth is, I don't think I've myself 'concluded' this line of thought either, so, in that way, now we're both suffering.

The update loop where I calculate forces and then apply them is, I later found out, called Euler integration. It is also one of the worst ways to numerically integrate in physics simulations because it doesn't conserve energy.

This is also why, if you keep the 'stable' three body orbit open for long enough, you may discover that it isn't so 'stable' after all. The reality, again, is that, in order to actually get an orbit that's stable forever,

(*) Ignoring the orbits issue, the gravity simulator I've coded here is still inaccurate due to the same issues just discussed. But it's probably accurate enough for 5-10 seconds and it looks cool regardless.

Also, if you want to play with the last simulation more, here's a full-screen link.

Thanks for reading!