{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Black Hole rendering with SageMath\n",
"\n",
"### Florentin Jaffredo"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"This notebook is a step-by-step implementation of a basic rendering engine in curved spacetime. The objective is to obtain a somewhat realistic image of the accretion disk around a black hole.\n",
"\n",
"The technique consists in launching lightlike geodesics toward the past from a single point (the virtual camera), using the [geodesic integrator](http://doc.sagemath.org/html/en/reference/manifolds/sage/manifolds/differentiable/integrated_curve.html) of SageMath. To reduce computation time, the spacetime is assumed be spherical symmetric; this reduces the number of required geodesics to produce an image of $n_x\\times n_y$ pixels from about $O\\left(n_x n_y\\right)$ to $O\\left(\\sqrt{n_x^2+n_y^2}\\right)$.\n",
"\n",
"This work relies heavily on the [SageManifolds Project](https://sagemanifolds.obspm.fr/). Advanced SageMath notions will also be used throughout this notebook, like Cython compilation and multithreading.\n",
"\n",
"This notebook requires a version of SageMath at least equal to 8.5:\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'SageMath version 8.7, Release Date: 2019-03-23'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"version()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"The code is separated into 9 parts.\n",
"\n",
"* Declaring the spacetime\n",
"* Launching a geodesic\n",
"* Launching a lot of geodesics!\n",
"* Figuring out where it intersects with the accretion disk\n",
"* Adding thickness to the disk\n",
"* Using black-body radiation and converting spectra to RGB\n",
"* First relativistic effect: Doppler effect\n",
"* Second relativistic effect: aberration (forward focalisation)\n",
"* Conclusion\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Configuration\n",
"\n",
"This notebook can be quite ressource hungry to run. For that reason different configurations options are provided. It is recommended to start with the lowest one to check that everything works properly. You can of course adapt the number of CPUs to your needs."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**First configuration:** will run in less than a minute on a 4-core laptop. Produces tiny images with no details (no secondary image)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# n_cpu = 4 # 4 Go Ram minimum\n",
"# n_geod = 100\n",
"# nx, ny = 180, 90"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Second configuration:** will run in about 5 minutes on a workstation, produces a reasonably sized image:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"n_cpu = 8 # 8 Go Ram minimum\n",
"n_geod = 1000\n",
"nx, ny = 720, 360"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Third configuration:** will run in 30 minutes on the Google Cloud Compute Engine. Produces a 4K image showing tiny details on the secondary disk images."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# n_cpu = 36 # 144 Go Ram minimum\n",
"# n_geod = 30000\n",
"# nx, ny = 4000, 2000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Additional preliminaries: display objects with $ \\LaTeX $ where possible, and silence deprecation warnings that arise from a few third-party packages:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"%display latex\n",
"import warnings\n",
"warnings.simplefilter('ignore', DeprecationWarning)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Declaring the spacetime\n",
"\n",
"Let's start slow by declaring the spacetime we'll use for rendering: it is the ***Schwarzschild spacetime***.\n",
"\n",
"It is important to use a coordinate system that is regular at the horizon. Here we use the *Eddington-Finkelstein coordinates*.\n",
"\n",
"Let $m$ be the mass of the black hole (that we'll take equal to 2 later). \n",
"\n",
"We also add a restriction to ensure that nothing touches the central singularity, and we set the metric $g$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"M = Manifold(4, 'M', structure='Lorentzian')"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"t: (-oo, +oo); r: (1, +oo); th: (0, pi); ph: (-oo, +oo)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C. = M.chart(r't r:(1,+oo) th:(0,pi):\\theta ph:\\phi')\n",
"C.coord_range()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"m = var('m')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[ 2*m/r - 1 2*m/r 0 0]\n",
"[ 2*m/r 2*m/r + 1 0 0]\n",
"[ 0 0 r^2 0]\n",
"[ 0 0 0 r^2*sin(th)^2]"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = M.metric()\n",
"g[0,0] = -(1-2*m/r)\n",
"g[0,1] = 2*m/r\n",
"g[1,1] = 1+2*m/r\n",
"g[2,2] = r^2\n",
"g[3,3] = (r*sin(th))^2\n",
"g[:]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"g = (2*m/r - 1) dt*dt + 2*m/r dt*dr + 2*m/r dr*dt + (2*m/r + 1) dr*dr + r^2 dth*dth + r^2*sin(th)^2 dph*dph"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also define a 3-dimensional Euclidean space $E$ to plot some results, using a map $\\phi: M \\rightarrow E$:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"M --> E^3\n",
" (t, r, th, ph) |--> (x, y, z) = (r*cos(ph)*sin(th), r*sin(ph)*sin(th), r*cos(th))"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E. = EuclideanSpace()\n",
"phi = M.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])\n",
"phi.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Launching a geodesic\n",
"\n",
"[Geodesic integration](http://doc.sagemath.org/html/en/reference/manifolds/sage/manifolds/differentiable/integrated_curve.html) was first implemented in SageMath in 2017 and perfected in 2018 to support fast integration and event handling (used to detect the singularity in our case).\n",
"\n",
"To introduce the method, let's plot an orbit around a black hole.\n",
"\n",
"To do that, we need to find a starting point $p$ as well as an inital velocity vector $v$. It can be quite troublesome to find a suitable one, but here is a free one:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"p = M((0, 14.98, pi/2, 0))\n",
"Tp = M.tangent_space(p)\n",
"v = Tp((2, 0, 0.005, 0.05))\n",
"v = v / sqrt(-g.at(p)(v, v))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$v$ is defined as a member of the tangent space at $p$. The last line is used to normalize $v$ as a unit timelike vector."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next is the definition of the geodesic. We need to pass a symbolic variable for the proper time (which will not be used). The starting point is deduced from the velocity vector (as the point where the velocity vector is defined)."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"tau = var('tau')\n",
"curve = M.integrated_geodesic(g, (tau, 0, 3000), v)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The integration should be very fast. Don't forget to give some numerical value to $m$ here."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"sol = curve.solve(step = 1, method=\"ode_int\", parameters_values={m: 2})\n",
"# sol = curve.solve(step = 1, parameters_values={m: 2})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting the solution requires an interpolation. This is automatically done in the next line."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"interp = curve.interpolate()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following cell plots the result using the mapping we provided previously. We also add a grey sphere at $r_s = 2m = 4$ (the event horizon) to give a scale."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"P = curve.plot_integrated(mapping=phi, color=\"red\", thickness=2, plot_points=3000)\n",
"P = P + sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"P.show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see that it look nothing like an ellipse, as we are used to in classical celestial mechanics. At this step, you can try adding an angular momentum to the black hole--in other words going from Schwarzschild to Kerr--by setting a non-zero angular momentum in the definition of the manifold ($J=1$ works fine). When this is the case, the orbits are not even included in a plane. Don't forget de revert back your changes before proceeding to the next part."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Launching a lot of geodesics!\n",
"\n",
"Of course one geodesic is not enough for us, we'll need at least a few hundred of them.\n",
"\n",
"Because we don't need to compute the equation again each time, we simply copy the previous declaration of the geodesic while changing the initial point and velocity.\n",
"\n",
"It will be useful here to introduce the Python module `multiprocessing` and progress bars as widgets:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"import multiprocessing\n",
"from ipywidgets import FloatProgress\n",
"from IPython.display import display"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It wouldn't be a great idea to set \"1 job = 1 geodesic integration\". Indeed, that would mean copying the geodesic declaration a few hundred times, which would be quite slow. What is done instead is seperating geodesics into batches using the following function:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def chunks(l, n):\n",
" \"\"\"Yield successive n-sized chunks from l.\"\"\"\n",
" for i in range(0, len(l), n):\n",
" yield l[i:i + n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The number of batches per CPU in not very important. If set to 1, some CPUs may run faster than other ones and stay idle at the end. If too high, too much time will be spent copying the curve setting. I found 3 to be a good value."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"n_batches_per_cpu = 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also redefine the previous geodesic to our new needs: fewer steps and the ability to check for chart boundaries when integrating. The $v$ in this case will not be used; it will always be overwritten before starting any integration."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"curve = M.integrated_geodesic(g, (tau, 0, 200), v, across_charts=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When using `multiprocessing`, functions can only accept a single argument. To overcome this limitation, each argument will be a tuple (curve, start index, number of curves to integrate)."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"args = []\n",
"start_index = 0\n",
"\n",
"for chunk in chunks(range(n_geod), n_geod//(n_batches_per_cpu*n_cpu)):\n",
" args += [(loads(curve.dumps()), start_index, len(chunk))]\n",
" start_index += len(chunk)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next line prints the list of arguments. We can check that each of the 100 geodesics are correctly set. Our little trick allowed us to only define 13 geodesics (about 3 per core, as we wanted; note, the exact result here will depend on what you used for `n_cpu` at the beginning)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(Integrated geodesic in the 4-dimensional Lorentzian manifold M, 984, 16)\n",
"25\n"
]
}
],
"source": [
"print(args[-1])\n",
"print(len(args))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now comes a question: which vector can be used as the starting 4-velocity?\n",
"\n",
"We need a past-oriented lightlike vector pointing toward the center but with a linearly increasing angle. The 3 space components are already imposed. The time component must then be chosen so that the total vector is lightlike.\n",
"\n",
"Let $p$ be the initial point and $v$ the initial 4-velociy, with an unknown time coordinate $dt$ ($y$ depends on the angle, it is a known quantity)."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"dt, y, r0 = var('dt, y, r0')"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"p = M((0, r0, pi/2, 0))\n",
"Tp = M.tangent_space(p)\n",
"v = Tp((dt, -1, 0, y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The norm of $v$ is currently given by:"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"r0^2*y^2 + dt^2*(2*m - r0)/r0 - 4*dt*m/r0 + (2*m + r0)/r0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.at(p)(v, v)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need to find $dt$ so that this expression is equal to 0 (lightlike condition). this is easy:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[dt == -(sqrt(-(2*m*r0 - r0^2)*y^2 + 1)*r0 - 2*m)/(2*m - r0), dt == (sqrt(-(2*m*r0 - r0^2)*y^2 + 1)*r0 + 2*m)/(2*m - r0)]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol = g.at(p)(v, v).solve(dt)\n",
"sol"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As expected, there are two solutions: one past-oriented and one future-oriented. In fact, in our case it does not matter, given that the Schwartzschild spacetime is static."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next cell defines the function that will be called by `multiprocessing`. It starts by unpacking the arguments, setting an empty dictionary as the result, and defining the starting position.\n",
"\n",
"The initial velocity is then overwritten using the formula above, the integration is performed, and the result is added to the dictionary."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"def calc_some_geodesics(args):\n",
" \"\"\"\n",
" Compute nb geodesics starting at index n0\n",
" \"\"\"\n",
" curve, n0, nb = args\n",
" res = {}\n",
" r = 100\n",
" posi = [0, r, pi/2, 0]\n",
" p = M(posi)\n",
" Tp = M.tangent_space(p)\n",
" for i in range(n0, n0+nb):\n",
" # starting vector\n",
" dy = i*0.006/n_geod\n",
" v = Tp([sol[0].rhs()(r0=r, y=dy, m=2).n(), -1, 0, dy])\n",
" # overwrite the starting vector\n",
" curve._initial_tangent_vector = v\n",
" # integration with m=2\n",
" curve.solve_across_charts(step=0.2, parameters_values={m:2})\n",
" # copy and clear solution\n",
" res[i] = (p.coord(), curve._solutions.copy())\n",
" curve._solutions.clear()\n",
" return res"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`geo` will keep the numerical solutions. I like to see `pool` as a hole in which I can throw some jobs. `multiprocessing` will then magically do them for me using every resource available on the computer."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "e691872fcfbb487b836f396f5cc38a19",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"FloatProgress(value=0.0, max=1000.0)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"geo = {}\n",
"pool = multiprocessing.Pool(n_cpu)\n",
"\n",
"# progress bar display\n",
"f = FloatProgress(min=0, max=n_geod)\n",
"display(f)\n",
"\n",
"for i, some_res in enumerate(pool.imap_unordered(calc_some_geodesics, args)): # do and wait\n",
" # progress bar update\n",
" f.value += len(some_res)\n",
" # update result\n",
" geo.update(some_res)\n",
"\n",
"# clean exit\n",
"pool.close()\n",
"pool.join()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If, for any reason, you don't want to use parallel computing, you can replace the previous cell with this one:"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"# geo = calc_some_geodesics((c, 0, n_geod))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now try to visualize those geodesics.\n",
"Next cell will plot 20 of them."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# add the sphere\n",
"P = sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"\n",
"#Â cycle through the solutions\n",
"for i in range(0, n_geod, 5*n_geod/100): \n",
" #Â set solution\n",
" curve._solutions = geo[i][1]\n",
" #Â do interpolation\n",
" interp = curve.interpolate()\n",
" # plot the curve\n",
" P += curve.plot_integrated(mapping=phi, color=[\"red\"], thickness=2, plot_points=150, \n",
" label_axes=False, across_charts=True)\n",
"\n",
"#Â show the result \n",
"P.show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that some fall inside the black hole toward the singularity. That's not an issue because the integration is automaticaly stopped when the geodesic leaves the chart domain defined in part 1."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Intersection with the accretion disk\n",
"\n",
"Time to transform those simulated light-rays into an image. To do this, we first need to compute the intersection between each geodesic and the accretion disk.\n",
"\n",
"For this example, the disk spans from $r=8$ to $r=50$, and is tilted by an angle $\\alpha = - \\frac{\\pi}{20}$."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"disk_min = 12\n",
"disk_max = 50\n",
"alpha = -pi/20"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's plot the disk on top of the last figure.\n",
"\n",
"(We cheat a little bit here and use a flattened torus.)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"D = sage.plot.plot3d.shapes.Torus((disk_min+disk_max)/2,\n",
" (disk_min-disk_max)/2).scale(1,1,0.01).rotateY(-pi/20)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"(P+D).show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The same but tilted on the X-axis by an angle $\\beta=\\frac{\\pi}{3}$. As explained earlier, the final image will be obtained by computing for each pixel : \n",
"\n",
"* Which geodesic best describes the light-ray\n",
"* Which angle $\\beta$ at which the disk should be tilted\n",
"* The intersection between the disk and that geodesic\n"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"