(blog/bhs) Added an initial explanation of math and physics behind it

11 files changed, 289 insertions, 0 deletions

diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/Schwarzschild-metric-in-Schwarzschild-coordinates-wikipedia.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/Schwarzschild-metric-in-Schwarzschild-coordinates-wikipedia.png
new file mode 100644
index 0000000..3964615
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/Schwarzschild-metric-in-Schwarzschild-coordinates-wikipedia.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/Schwarzschild-radius-wikipedia.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/Schwarzschild-radius-wikipedia.png
new file mode 100644
index 0000000..fccf54f
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/Schwarzschild-radius-wikipedia.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/algorithm.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/algorithm.png
new file mode 100644
index 0000000..40d4c86
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/algorithm.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/antonelli-formula-expanded.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/antonelli-formula-expanded.png
new file mode 100644
index 0000000..48eb7d4
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/antonelli-formula-expanded.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/isco-wikipedia.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/isco-wikipedia.png
new file mode 100644
index 0000000..2b01135
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/isco-wikipedia.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/normalize-vector-wikipedia.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/normalize-vector-wikipedia.png
new file mode 100644
index 0000000..517e142
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/normalize-vector-wikipedia.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/photon-sphere-wikipedia.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/photon-sphere-wikipedia.png
new file mode 100644
index 0000000..840572c
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/photon-sphere-wikipedia.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/riccardo-antonelli-schwarzschild-black-hole-force-field-formula.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/riccardo-antonelli-schwarzschild-black-hole-force-field-formula.png
new file mode 100644
index 0000000..846b0d3
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/riccardo-antonelli-schwarzschild-black-hole-force-field-formula.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/specific-angular-momentum-wikipedia.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/specific-angular-momentum-wikipedia.png
new file mode 100644
index 0000000..61b3383
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/specific-angular-momentum-wikipedia.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/img/variables-expanded.png b/src/blog/difficulties-of-making-a-black-hole-shader/img/variables-expanded.png
new file mode 100644
index 0000000..484e2f1
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/img/variables-expanded.png
diff --git a/src/blog/difficulties-of-making-a-black-hole-shader/index.md b/src/blog/difficulties-of-making-a-black-hole-shader/index.md
new file mode 100644
index 0000000..7c78c05
--- /dev/null
+++ b/src/blog/difficulties-of-making-a-black-hole-shader/index.md
@@ -0,0 +1,289 @@
+# Difficulties of making a black hole shader
+
+## Contents
+
+This year I took a university elective course for WebGL programming.
+Back in August I had to create and submit some sort of a project, so I chose to create an interractive and animated black hole demo.
+The primary idea was to show how the gravitational pull of black holes bend light and how the result looks from different points of view.
+
+During development I found a lot of hurdles, primarily with the actual bending of light.
+Information online is kind of hard to understand and find, so I'll try my best to lay out everything I know.
+
+I'll assume you don't understand a lot of physics, but are familiar with vectors and operations with them.
+Also this is in no way an exploration of best practices and good WebGL development.
+
+## Context
+
+This section serves as a quick refresher.
+
+### Brief introduction to computer graphics
+
+Imagine a photo of a bottle: you can see the bottle, but you cannot rotate it, or change it's size, or modify it in any way, without having to download the image, fire up a comlicated image editing program and spend hours upon hours doing your modifications.
+And it becomes even harder if you want to add or remove objects, say a table on which the bottle sits.
+
+The main idea of computer graphics is to define your objects with some few but "abstract" values: the bottle is X tall, the base is Y wide, it has Z color, etc. and then a computer program has to take those values and "generate" the actual image.
+Or to be more precise, for every pixel on the screen where you should be able to see the bottle, the computer calculates what the precise color of that pixel is.
+
+This allows us to dynamically change objects: if I wanna make the bottle shorter, I make X a smaller number, and the computer recalculates everything using the new height.
+If I want to add or remove objects, I just give or remove the specific parameters that define those objects and tell the computer to recalculate the image again.
+
+As you might imagine, it's extremely difficult to calculate everything needed to show an image, and it takes an absurd amount of system resources and time.
+That's why people decided to create a separate processing unit, just for rendering graphics: the [GPU](TODO).
+
+A GPU is optimized on the hardware level for parallel computation (calculate multiple stuff at the same time), fast math with 3D points and matricies, and so on.
+The tradeoff is that a GPU doesn't implement most of the normal stuff a [CPU](TODO) does.
+
+For this reason, programs written for GPUs are called "shaders".
+Most often shaders are written in some "higher" level programming language, which is then translated via a graphics API or library to something the actual GPU can understand.
+For comparison, in normal softare development you write a program (not a shader), which is then translated to something the CPU can understand via a compiler or interpreter.
+
+### Brief introduction to WebGL
+
+One of the most popular such libraries is [OpenGL](TODO) with it's associated shader programming languages [GLSL](TODO) and [GLSL ES](TODO).
+In the 2000s people decided it would be a good idea to add a graphics API/library to web browsers.
+Skipping over a lot of history, [WebGL](TODO) was born and widely adopted.
+
+It serves as a "layer above" OpenGL: you write JavaScript code with which the browser tells OpenGL what to do.
+You would get a DOM object in which the graphics will be rendered, and a "context" object to which you'll pass down your shaders, shader parameters and change graphics settings.
+
+## Theory behind black holes
+
+As we know from physics, objects with mass create a gravitational field which "pulls" inwards everything around them.
+However, [Einstein's general theory of relativity](https://en.wikipedia.org/wiki/General_relativity) showed that this isn't really a force which pulls, like with (eletro)magnetism, but more of an effect: mass changes the curvature of [spacetime](TODO), which by extension changes the geodesic's "curavture".
+
+### Explanation of geodesics and gravitational lensing
+
+Let's back up a bit, first think of a flat surface, like a table or a (flat) sheet of paper.
+What is the shortest distance between two points?
+Obviously, that's a straight line.
+
+However, what if our surface isn't flat?
+Imagine two points on a sphere, what is the shortest distance between them?
+You might be inclined to say that's also a straight line, but it really isn't: it's an arc of a circle!
+
+![](https://upload.wikimedia.org/wikipedia/commons/thumb/c/cb/Illustration_of_great-circle_distance.svg/480px-Illustration_of_great-circle_distance.svg.png)  
+[*Image by CheCheDaWaff on Wikipedia*](https://commons.wikimedia.org/w/index.php?curid=48187293)
+
+This circle is called a [great circle](https://en.wikipedia.org/wiki/Great_circle).
+Looking at all possible great circles, we find they have an interesting property: their center (and radius) matches with the center (and radius) of the sphere.
+
+Another big difference between this arc and a "line" is that lines can be parallel, but arcs cannot.
+Again, the shortest distance between to points is a line *segment*, extend that segment into infinity and it becomes a line.
+On a flat surface, we know that if you take two points and create the line passing through them, then there exists at least two other distinct points, whose line, crated from intersecting them, does not ever intersect the first line.
+However, as we said, on a sphere all great circles have the same center and radius, so it doesn't matter which two pairs of points we take, their great circles must always intersect[^parallels-of-latitude]!
+
+In non-flat surfaces, the segment which has the shortest distance between two points is called a geodesic (in flat surfaces, a line is practically a geodesic, on spheres, great circles are geodesics).
+We can now also define what it means to "go straight forwards": in flat surfaces we would say that means following a straight line, while on non-flat it means following a geodesic.
+
+Going back to physics, what the General theory of relativity tells us is that gravity bends spacetime, the actual space, or 3D surface, in which the universe exists.
+In effect, the geodesics in our space are also curved.
+Therefore, if you are always travelling forward, perfectly following a geodesic, if you start passing an object of great mass, it's gravity will bend that geodesic, so even though you are always going forwards, you'll actually "turn" or follow a curve.
+
+This is the exact effect a black hole produces: it's mass is so great and it bends space(time) **so much**, light can take a really sharp curve.
+We call this [graviational lensing](TODO):
+
+![](https://upload.wikimedia.org/wikipedia/commons/0/02/Gravitational_lens-full.jpg)  
+*(notice how the gray arrows show the path light travels from a galaxy that is behind (obstructed by) the black hole)*
+
+I spent a good amount of time explaning that gravity doesn't curve object's path but spacetime itself, however, this distinction really matters for physics, in our simulation it requires significantly less computation and programmer effort if we assume the former.
+
+Then were all of these explanations a waste of time?
+Not really, since physics works with the spacetime curvature, we'll have to do conversions/modifications to the formulae, so they behaves like a force that pulls on the object in motion.
+
+### Finding a formula
+
+Now, how do we calculate the exact curvature a black hole produces?
+General theory of relativity provides us with the [Einstein field equations](https://en.wikipedia.org/wiki/Einstein_field_equations), a complicated set of equations which allows us to calculate how spacetime curves around an object with large mass.
+[Schwartzschild](TODO) solved these equations, assuming that the body with mass doesn't spin[^schwm-assumptions].
+This solution is called the [Schwarzschild metric](https://en.wikipedia.org/wiki/Schwarzschild_metric).
+A black hole for which this (these) assumption is true is called a "Shwartzschild black hole".
+
+This is the simplest solution to those equation, and a Schwartzschild black hole is the simplest type of black hole, so we're gonna use that.
+Now, by definition, here is the Schwartzschild metric in Schwartzschild coordinates:
+
+![](./img/Schwarzschild-metric-in-Schwarzschild-coordinates-wikipedia.png)
+
+**YOUCH!**
+
+You might notice there are derivatives in this equation, as you might imagine, it's unbearably slow to calculate it, especially if we want our demo to run at at least 30 frames per second!
+There is also another problem: Schwartzschild coordinates, our GPU works with normal cartesian coordinates, so we'll also have to do transform the coordinates.
+We'll also need to transfrom the curved spacestime coordinates to cartesian, as if spacetime wasn't curved around the black hole.
+
+I personally have zero clue how to deal with this monster.
+Thankfully [Riccardo Antonelli](http://rantonels.github.io/starless/) was able to do all of our coordinate conversions and simplify the living crap out of that monster (using some approximations), turning it in a pretty small equation:
+
+![](./img/riccardo-antonelli-schwarzschild-black-hole-force-field-formula.png)
+
+Question now is: how do we use it?
+This turned out to be a much harder problem than it should have, I actually had to refer to [Dmitry Brant](https://dmitrybrant.com/2018/12/11/ray-tracing-black-holes)'s usage of [the formula for his own project](https://github.com/dbrant/BlackHoleRaytracer) as well as [Ricardo Antonelli's implementation](https://github.com/rantonels/starless/blob/7a0c3c4efdc0666f53dcfa48d310fb2164be2d1a/tracer.py).
+
+### Applying the formula to "bend" light's path
+
+First, how are we going to represent light?
+The simplest way is to think of light as a small particle (or spaceship) travelling through space.
+We'll move our particle by a tiny amount, and before moving it, we'll calculate how much the black hole changes it's direction and move the particle in that new direction.
+
+Referring at the formula again:
+
+![](./img/riccardo-antonelli-schwarzschild-black-hole-force-field-formula.png)
+
+The <code>F<sup>&rarr;</sup></code> returns a force vector (`F` for "force" and <code><sup>&rarr;</sup></code> for "vector"): any force acting on an object has a strength (magnitude) and a direction.
+Imagine you are pushing a box: you are applying some force to it in a "forward" direction.
+If you pull the box, you'll apply some force in a "backward" direction.
+
+This force vector specifies how much and in what direction exactly the black hole is pulling on the light "particle".
+
+Although the function takes in only one argument, `r`, that doesn't mean that `h` and <code>r<sup>\^</sup></code> are some universal constants: they also need to be calculted and "passed" to the function.
+
+`r` is the distance from the center of the black hole (the singularity) to the **current** position of the light "particle".
+
+`h` is the magnitude of the ["specific angular momentum"](https://en.wikipedia.org/wiki/Specific_angular_momentum) of the light "particle".
+
+"Normal" linear momentum specifies how much "force" at what speed it takes to stop (or change) the speed and direction of an object in linear motion (that's why it's in kg\*m/s).
+Imagine a big heavy train going at a fast speed.
+How do you forcefully stop it?
+
+You can put a *really* heavy barrier on the tracks, in front of the train.
+Since the train has a lot of momentum, this barrier will have to be much much heavier (with much more mass) than the train itself (when the train would be stationary).
+
+Alternatively, you can "launch" something against the train, meaning you can take a second train, and make it go in the opposite direction of the first train, so they hit head-on.
+Since the second train would also be moving, that would add to it's momentum, and the second train wouldn't need to be much heavier than the first train.
+In fact, if the second train travels at the same speed as the first, then the mass of both trains would suffice to be equal, and if it travels twice as fast, it would suffice if the second train has half the mass of the first.
+
+Angular momentum specifies how much "force" at what speed it takes to stop (or change) the speed and directino of a rotating object.
+Imagine the train scenario but on a circular track.
+
+The specific angular momentum is the angular momentum of an object, divided by it's mass.
+It kind of tells us which property is more prevalent in the momentum, the mass or the speed: for the same angular momentum, a high *specific* angular momentum means our object is small and moves fast, while for a small *specific* angular momentum out object is large and moves slowly.
+
+In our case a light "particle" is being pulled by a stationary object of large mass (a Shwarzschild black hole), so rather than finding the angular momentum[^angular-momentum] and then dividing by the mass[^light-has-no-mass], we'll use another formula:
+
+![](./img/specific-angular-momentum-wikipedia.png)
+
+where **in this specific formula only** `r` is the vector from the center of the big mass (the singularity of the black hole) to the light "particle" and `v` is the velocity (direction and amount of speed) the light is travelling with.
+
+The cross product returns a vector, so we'll need to take the magnitutde of the result and plug that into the formula.
+It's also important that **we calculate it only once, in the beginning**.
+So, we shouldn't be calculating the magnitude of the specific angular momentum every time the light "particle" is moved, rather we calculate it once, after it gets close enough to be affected by the black hole.
+
+<code>r<sup>&#94;</sup></code> is the vector from the center of the black hole to the light particle as a [unit](https://en.wikipedia.org/wiki/Unit_vector) vector (hence the `^` on top).
+
+This means that we normalize the vector, meaning divide it by it's magnitude (length):
+
+![](./img/normalize-vector-wikipedia.png)
+
+Finally, putting all needed variables together in one formula, we get:
+
+<!--
+\begin{align*}
+& \overrightarrow{F}(r_l,h) = -\frac{3}{2} \; h^2 \; \frac{\hat{r_l}}{\|r_l\|^5} \quad \text{, where}\\
+& h = \|r_l \cross v_l\|\\
+& \hat{r_l} = \frac{r_l}{\|r_l\|}\\
+& \|x\| = \text{magnitude of }x\\
+& r_l = \text{vector from black hole center to light "particle"}\\
+& v_l = \text{velocity of light "particle"}
+\end{align*}
+-->
+![](./img/variables-expanded.png)
+
+where velocity of light "particle" means a vector, with the direction in which the light "particle" will travel, and a magnitude equal to the distance (speed) the light "particle" will travel in.
+As you may notice, this can be simplified, if we substitude <code>r<sup>&#94;</sup></code>:
+
+<!--
+\begin{align*}
+& \overrightarrow{F}(r_l,h) = -\frac{3}{2} \; h^2 \; \frac{r_l}{\|r_l\|^6} \quad \text{, where}\\
+& h = \|r_l \cross v_l\| \text{ (calculated right before being affected by black hole)}\\
+& \|x\| = \text{magnitude of }x\\
+& r_l = \text{vector from black hole center to light "particle"}\\
+& v_l = \text{velocity of light "particle"}
+\end{align*}
+-->
+![](./img/antonelli-formula-expanded.png)
+
+Finally, we know what to give the formula, and what it will return.
+Current question is: how do we use it's output?
+For that we'll need to get familiar with black holes.
+
+### Explaning black holes
+
+Physics regards black holes as an infinitessimal (infinitely small) point with some huge amount of mass.
+From the mass we can derive three radii, which will be the backbone of our black hole: the event horizon, the photon sphere and [ISCO](https://en.wikipedia.org/wiki/Innermost_stable_circular_orbit).
+
+The event horizon, also called the Schwarzschild radius, is the distance from the center of the black hole (singularity), which once passed, cannot be escaped.
+In other words, once something gets too close, and passes this distance, then it cannot ever escape, even if it's light.
+The formula for it is:
+
+![](./img/Schwarzschild-radius-wikipedia.png)
+
+where `G` is the [gravitational constant](https://en.wikipedia.org/wiki/Gravitational_constant) and `M` is the mass of the black hole (in kg).
+
+The photon sphere is the distance from the singularity, in which photons start orbiting the black hole (and can either escape or pass the event horizon).
+Formula:
+
+![](./img/photon-sphere-wikipedia.png)
+
+where `G` and `M` are the same as before, `c` is the speed of light in vacuum and <code>r<sub>s</sub></code> is the afroementioned Schwarzschild radius.
+
+Finally, the innermost stable circular orbit (ISCO), officially it's the distance, before which you can orbit the black hole, but once passed, you'll get sucked in and either get shot out or forever lost past the event horizon.
+For our purposes, it's fine to approximate it as the maximum radius after which spacetime isn't bent[^isco] and objects aren't affected.
+Formula:
+
+![](./img/isco-wikipedia.png)
+
+Putting them all together, we should expect light to travel like this:
+![](https://upload.wikimedia.org/wikipedia/commons/2/27/Black_Hole_Shadow.gif)
+*Animation by Hugo Spinelli on [Wikipedia](https://commons.wikimedia.org/w/index.php?curid=138168679)*
+
+### Putting everything together
+
+Finally, now we have everything we need to divise an algorithm for our light.
+For light we'll store two vectors: one for it's position and one for it's velocity.
+For a black hole we'll store a vector and a quantity: the position of the black hole's center and it's mass.
+For convinience, we could just store the Schwarzschild radius, since all of our radii depend only on it.
+
+Once the distance between the light and the center is less than or equal to the ISCO, we can begin augmenting light's velocity.
+On small steps, we move the light forward a bit with the velocity and then update the velocity with the black hole's pull.
+
+If light enters the event horizon, then we know it cannot escape, and stop moving it.
+If it passes the ISCO then it has left the black hole and will never approach it again.
+
+If light enters the photon sphere, then we'll have to "track" how long it's been there, if it decides to orbit in a stable manner, then it will stay there indefinetly.
+The simplest way to "track" this is just to limit how many times we can move the light forward.
+
+An algorithm doing this would look a bit like this (in pseudocode, inspired by Python):
+
+<!--
+\begin{align*}
+& p_{cbh} = \text{position of center of black hole}\\
+& r_s = \text{Schwarzchild distance from } p_{cbn} \\
+& r_{ISCO} = 3 r_s \\
+\\
+& p_l = \text{position of light "particle"}\\
+& v_l = \text{velocity of light "particle"}\\
+\\
+& r_l = \| p_l-p_{cbh} \| \quad \text{(distance between black hole and light)}\\
+& \text{If } r_l \leq r_{ISCO} : \\
+& \quad h = \| r_l \cross v_l \|\\
+& \quad \text{Do maximum 100 times}:\\
+& \qquad \text{If } r_l < r_s : \\
+& \qquad \quad \text{return "Light never escapes"} \\
+& \qquad \text{If } r_l > r_{ISCO} : \\
+& \qquad \quad \text{return "Light escaped black hole"} \\
+& \qquad p_l = p_l + v_l \quad \text{(move light)}\\
+& \qquad r_l = \| p_l-p_{cbh} \| \quad \text{(update distance to $p_{cbh}$)}\\
+& \qquad v_l = v_l + \overrightarrow{F}(r_l,h) \quad \text{(update light velocity)} \\
+& \quad \text{return "Light stuck in photon sphere, never escapes"}
+\end{align*}
+-->
+
+![](./img/algorithm.png)
+
+as a reminder, here is the definiton of <code>F<sup>&rarr;</sup></code>:
+
+![](./img/antonelli-formula-expanded.png)
+
+[^parallels-of-latitude]: You might think about the latituginal lines on a globe, they are indeed parallel (equidistant), but if you take any two points on the same such line, the distance between them won't be the smallest.
+[^schwm-assumptions]: He also assumed the body didn't have an electric charge and the [cosmological constant](TODO) is zero. These are way outside of my understanding and our purposes, so I skiped them for simplicity.
+[^isco]: It's not really true that gravity stopps affecting spacetime after the ISCO, even though things can orbit Earth that doesn't mean that Earth doesn't bend spacetime. However, for the purposes of a simple shader, this assumption produces perfectly useable results.
+[^angular-momentum]: Angular momentum is actually [defined using the cross product](https://en.wikipedia.org/wiki/Specific_angular_momentum#Definition)
+[^light-has-no-mass]: Actually, light has no mass, so this version of the formula wouldn't work. You might think that makes no sense, since we cannot divide by zero, so the specific angular momentum shouldn't exist, but then with the cross product we find a value for it. You would be correct: Antonelli does some funky approximations, combining Newtonian physics with General theory of relativity, to make the formula that simple. Someone in physics will probably vomit looking at it, but for our purposes it's more than perfect.