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.
The technique consists in launching lightlike geodesics toward the past from a single point (the virtual camera), using the geodesic integrator 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)$.
This work relies heavily on the SageManifolds Project. Advanced SageMath notions will alo be used throughout this notebook, like Cython compilation and multithreading.
This notebook requires a version of SageMath at least equal to 8.5:
version()
'SageMath version 8.7, Release Date: 2019-03-23'
The code is separated into 9 parts.
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.
First configuration: will run in less than a minute on a 4-core laptop. Produces tiny images with no details (no secondary image).
# n_cpu = 4 # 4 Go Ram minimum
# n_geod = 100
# nx, ny = 180, 90
Second configuration: will run in about 5 minutes on a workstation, produces a reasonably sized image:
n_cpu = 8 # 8 Go Ram minimum
n_geod = 1000
nx, ny = 720, 360
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.
# n_cpu = 36 # 144 Go Ram minimum
# n_geod = 30000
# nx, ny = 4000, 2000
%display latex
Let's start slow by declaring the spacetime we'll use for rendering: it is the *Schwarzschild spacetime*.
It is important to use a coordinate system that is regular at the horizon. Here we use the Eddington-Finkelstein coordinates.
Let $m$ be the masse of the black hole (that we'll take equal to 2 later).
We also add a restriction to ensure that nothing touches the central singularity, and we set the metric $g$.
M = Manifold(4, 'M', structure='Lorentzian')
C.<t, r, th, ph> = M.chart(r't r:(1,+oo) th:(0,pi):\theta ph:\phi')
C.coord_range()
m = var('m')
g = M.metric()
g[0,0] = -(1-2*m/r)
g[0,1] = 2*m/r
g[1,1] = 1+2*m/r
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g[:]
g.display()
We also define a 3-dimensional Euclidean space $E$ to plot some results, using a map $\phi: M \rightarrow E$:
E.<x, y, z> = EuclideanSpace()
phi = M.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
phi.display()
Geodesic integration 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).
To introduce the method, let's plot an orbit around a black hole.
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:
p = M((0, 14.98, pi/2, 0))
Tp = M.tangent_space(p)
v = Tp((2, 0, 0.005, 0.05))
v = v / sqrt(-g.at(p)(v, v))
$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.
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).
tau = var('tau')
curve = M.integrated_geodesic(g, (tau, 0, 3000), v)
The integration should be very fast. Don't forget to give some numerical value to $m$ here.
sol = curve.solve(step = 1, method="ode_int", parameters_values={m: 2})
# sol = curve.solve(step = 1, parameters_values={m: 2})
Plotting the solution requires an interpolation. This is automatically done in the next line.
interp = curve.interpolate()
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.
P = curve.plot_integrated(mapping=phi, color="red", thickness=2, plot_points=3000)
P = P + sage.plot.plot3d.shapes.Sphere(4, color='grey')
P.show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)
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 an 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.
Of course one geodesic is not enough for us, we'll need at least a few hundred of them.
Because we don't need to compute the equation again each time, what is done is simply copying the previous declaration of the geodesic while changing the initial point and velocity.
It will be usefull here to introduce the Python module multiprocessing
and progress bars as widgets:
import multiprocessing
from ipywidgets import FloatProgress
from IPython.display import display
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 in batches using the following function:
def chunks(l, n):
"""Yield successive n-sized chunks from l."""
for i in range(0, len(l), n):
yield l[i:i + n]
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 spend copying the curve setting. I found 3 to be a good value.
n_batches_per_cpu = 3
We also redefine the previous geodesic to our new needs: fewer steps and the ability to check for chart boundaries when integrating. The $v$ used here will not be used, it will always be overwritten before starting any integration.
curve = M.integrated_geodesic(g, (tau, 0, 200), v, across_charts=True)
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 curve to integrate).
args = []
start_index = 0
for chunk in chunks(range(n_geod), n_geod//(n_batches_per_cpu*n_cpu)):
args += [(loads(curve.dumps()), start_index, len(chunk))]
start_index += len(chunk)
The next line prints the list of arguments. We can check that each of the 100 geodesic is correctly set. Our little trick allowed us to only define 13 geodesics (about 3 per core, as we wanted)
print(args[-1])
print(len(args))
(Integrated geodesic in the 4-dimensional Lorentzian manifold M, 984, 16) 25
Now comes a question: which vector can be used as the starting 4-velocity?
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.
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).
dt, y, r0 = var('dt, y, r0')
p = M((0, r0, pi/2, 0))
Tp = M.tangent_space(p)
v = Tp((dt, -1, 0, y))
The norm of $v$ is currently given by:
g.at(p)(v, v)
We need to find $dt$ so that this expression is equal to 0 (lightlike condition). this is easy:
sol = g.at(p)(v, v).solve(dt)
sol
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.
Next cell defines the function that will be called by multiprocessing
. It starts by unpacking the arguments, setting an empty dictionnary as the result, and defining the starting position.
The initial velocity is then overwritten using the formula above, the integration is performed, and the result is added to the dictionnary.
def calc_some_geodesics(args):
"""
Compute nb geodesics starting at index n0
"""
global sol, M
curve, n0, nb = args
res = {}
r = 100
posi = [0, r, pi/2, 0]
p = M(posi)
Tp = M.tangent_space(p)
for i in range(n0, n0+nb):
# starting vector
dy = i*0.006/n_geod
v = Tp([sol[0].rhs()(r0=r, y=dy, m=2).n(), -1, 0, dy])
# overwrite the starting vector
curve._initial_tangent_vector = v
# integration with m=2
curve.solve_across_charts(step = 0.2, parameters_values={m:2})
# copy and clear solution
res[i] = (p.coord(), curve._solutions.copy())
curve._solutions.clear()
return res
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 ressource available on the computer.
geo = {}
pool = multiprocessing.Pool(n_cpu)
# progress bar display
f = FloatProgress(min=0, max=n_geod)
display(f)
for i, some_res in enumerate(pool.imap_unordered(calc_some_geodesics, args)): # do and wait
# progress bar update
f.value += len(some_res)
# update result
geo.update(some_res)
# clean exit
pool.close()
pool.join()
FloatProgress(value=0.0, max=1000.0)
If, for any reason, you don't want to use parallel computing, you can replace the previous cell with this one:
# geo = calc_some_geodesics((c, 0, n_geod))
We can now try to visualize those geodesics. Next cell will plot 20 of them.
# add the sphere
P = sage.plot.plot3d.shapes.Sphere(4, color='grey')
# cycle through the solutions
for i in range(0, n_geod, 5*n_geod/100):
# set solution
curve._solutions = geo[i][1]
# do interpolation
interp = curve.interpolate()
# plot the curve
P += curve.plot_integrated(mapping=phi, color=["red"], thickness=2, plot_points=150,
label_axes=False, across_charts=True)
# show the result
P.show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)
We can see that some falls 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.
Time to transform those simulated light-rays into an image. To do this, we need to first compute the intersection between each geodesic and the accretion disk.
For this exemple, the disk span from $r=8$ to $r=50$, and is tilted by an angle $\alpha = - \frac{\pi}{20}$.
disk_min = 12
disk_max = 50
alpha = -pi/20
Let's plot the disk on top of the last figure.
(We cheat a little bit here and use a flatten torus).
D = sage.plot.plot3d.shapes.Torus((disk_min+disk_max)/2,
(disk_min-disk_max)/2).scale(1,1,0.01).rotateY(-pi/20)
(P+D).show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)
The same but tilted on the X-axis by an angle $\beta=\frac{\pi}{3}$. As explained earlyer, the final image will be obtained by computing for each pixel :
(P+D.rotateX(pi/3)).show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)
The geodesics are formatted in a strange way because of the solver used. The next line makes it easier to use. geo
is now a list of list of coordinates (and not a dictionnary of strange things).
geo = [list(geo[i][1].values())[0][0][1].tolist() for i in range(len(geo))]
To detect the intersection between the disk and a geodesic, the only solution is to parse the list of successive coordinates. This is done in the following function.
For each point of the curve, two rotations are performed (manually for speed purposes) before checking the point coordinates.
def intersection(curve, alpha, beta):
"""
Return True if the curve intersect the disk comprised between dmin and dmax
tilted of angles alpha and beta
"""
n = len(curve)
r, theta, phi = curve[0][2:5]
x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)
z = z*cos(alpha)+x*sin(alpha)
for i in range(1, n):
# done in 3 lines for speed consideration
r = curve[i][2]
theta = curve[i][3]
phi = curve[i][4]
# conversion to cartesian:
x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
# rotation around the X-axis:
y2, z2 = y2*cos(beta)-z2*sin(beta), z2*cos(beta)+y2*sin(beta)
# rotation around the Y-axis:
x2, z2 = x2*cos(alpha)-z2*sin(alpha), z2*cos(alpha)+x2*sin(alpha)
if z!=z2: # needed to prevent a division by zero next line
t = z/(z-z2) # if 0<=t<1 then the curve intersect the disk between the points i and i-1
if t>=0 and t<1 and curve[i][2]>disk_min and curve[i][2]<disk_max: #check coordinates
return True
x, y, z = x2, y2, z2
return False
The only problem with this function is the speed at which it operates. Let's try it on an example:
%time intersection(geo[28], alpha=pi/5, beta=0)
CPU times: user 1.6 s, sys: 36.7 ms, total: 1.64 s Wall time: 1.58 s
On my computer, it takes more than one second. You can imagine that if this must be done for each pixel, we would have to wait for days for even a single image.
By the way, you can argue that everything this function does could be done more efficiently, using for example numpy
matrix operations. That's true, but in the next parts we'll do more inside the loop than just return True
.
The solution: compile this function!
Sage let us use Cython quite freely, so let's do that.
%%cython
from libc.math cimport cos, sin
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef intersection(list curve, float alpha, float beta, float dmin, float dmax):
"""
Return True if the curve intersect the disk comprised between dmin and dmax
tilted of angles alpha and beta
"""
cdef float x, y, z
cdef float x2, y2, z2
cdef float r, theta, phi
cdef int n, i
cdef float t
cdef float sinalpha, cosalpha
cdef float sinbeta, cosbeta
sinalpha = sin(alpha)
cosalpha = cos(alpha)
sinbeta = sin(beta)
cosbeta = cos(beta)
n = len(curve)
r, theta, phi = curve[0][2:5]
x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)
z = z*cos(alpha)+x*sin(alpha)
for i in range(1, n):
r = curve[i][2]
theta = curve[i][3]
phi = curve[i][4]
x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
y2, z2 = y2*cosbeta-z2*sinbeta, z2*cosbeta+y2*sinbeta
x2, z2 = x2*cosalpha-z2*sinalpha, z2*cosalpha+x2*sinalpha
if z!=z2:
t = z/(z-z2)
if t>=0 and t<1 and curve[i][2]>dmin and curve[i][2]<dmax:
return True
x, y, z = x2, y2, z2
return False
%time intersection(geo[28], pi/5, 0, disk_min, disk_max)
CPU times: user 291 µs, sys: 0 ns, total: 291 µs Wall time: 292 µs
Boom!
x10000 speedup in my case!
Time to automatize the process and make an image. We'll once again use multiprocessing
. One job this time will correspond to one row of pixel. This time we'll use numpy
.
The result will be an array of RGB colors.
import numpy as np
import matplotlib.pyplot as plt
data = np.zeros( (ny, nx, 3), dtype=float )
def render_row(x):
"""
Render a single row of the image
"""
res = np.zeros((ny,3)) # result row in RGB format
for y in range(ny): # for each pixel in the row
beta = atan2(y-ny/2,x) # beta angle
r = sqrt(x**2+(y-ny/2)**2) # pixel distance to the center of the image
ind_geo = int(r/400*n_geod*720/nx) # index of the geodesic to use. values are obtained by trial and error.
if ind_geo<n_geod: # don't bother if the index is too big.
if intersection(geo[ind_geo], float(pi/20), float(beta),
float(disk_min), float(disk_max)):
res[ny-1-y,:] = [255, 0, 0] # red if intersects; black otherwise
return x, res # also return the argument, which is not remembered otherwise
This time, it is very important to use the argument maxtasksperchild
to kill the unused processes. This prevent a memory leak.
We enclose it inside a function because we'll be calling it multiple times.
def render():
global data
# display progress bar
f = FloatProgress(min=0, max=nx)
display(f)
pool = multiprocessing.Pool(n_cpu, maxtasksperchild=int(10))
for i, rows in enumerate(pool.imap_unordered(render_row,
range(-nx/2, nx/2))): # do and wait
# update progress bar
f.value+=1
# copy the result in the right place
x, res = rows
data[:, x+nx/2, :] = res
# clean exit
pool.close()
pool.join()
render()
FloatProgress(value=0.0, max=720.0)
Using SciPy, we can now convert it to an actual image:
import scipy.misc as smp
img1 = smp.toimage(data)
img1
/home/eric/sage/8.7/local/lib/python2.7/site-packages/sage/repl/ipython_kernel/__main__.py:3: DeprecationWarning: `toimage` is deprecated! `toimage` is deprecated in SciPy 1.0.0, and will be removed in 1.2.0. Use Pillow's ``Image.fromarray`` directly instead. IPKernelApp.launch_instance(kernel_class=SageKernel) /home/eric/sage/8.7/local/lib/python2.7/site-packages/scipy/misc/pilutil.py:381: DeprecationWarning: `bytescale` is deprecated! `bytescale` is deprecated in SciPy 1.0.0, and will be removed in 1.2.0. bytedata = bytescale(data, high=high, low=low, cmin=cmin, cmax=cmax)
If you have seen the movie Interstellar, you will know that this looks a bit like a (monochrome) black hole, a very simple one for sure, but we'll improve it step by step int the following parts.
First we will be adding transpareny and thickness to the disk. For that, let's assume the disk is optically thin. This means light can only be added and never obstructed. What's more, the intensity will be proportional to the total length of emmisive medium traversed.
This length depends on 2 factors : the profile of the disk and the angle $\theta$ between the disk and the light ray:
$\qquad d = \frac{f(r)}{sin(\theta)}$
The computation of the angle $\theta$ is not trival. The fastest way to obtain it is the perform a change of frame, which will locally (at a single point) give us a Minkowsky metric (orthonormal frame), in which angles are easy to compute.
To remind you, here is our metric:
g[:]
Or for a point of the disk:
r0, phi = var('r_0, phi')
p = M((0, r0, pi/2, phi))
g.at(p)[:]
We are not that far off. Some rescaling and mixing the first two lines should do the trick.
# default frame
fr = C.frame()
# create an automorphism field
aut = M.automorphism_field()
# some symbolic variables
a, b, c = var('a, b, c')
# let's try with the simplest matrix possible
aut.add_comp()[:] = [[a, 0, 0, 0], [b, c, 0, 0], [0, 0, 1/r0, 0],
[0, 0, 0, 1/r0]] # only b is off-diagonal
fr2 = fr.new_frame(aut, 'f2')
In this new frame, the metric at $p$ looks like this:
g.at(p)[fr2.at(p), :]
This lets us find the value of a,b and c:
c = sqrt(r/(r+2*m))
a = sqrt(((r+2*m)/(r)))
b = -2*a*m/(2*m+r)
Let's check:
aut2 = M.automorphism_field() # new automorphism field
aut2.add_comp()[:] = [[a, 0, 0, 0], [b, c, 0, 0], [0, 0, 1/r, 0], [0, 0, 0, 1/(r*sin(th))]]
fr3 = fr.new_frame(aut2, 'f3')
g[fr3, :]
Success!
We now have an orthonormal frame everywhere (this is only possible because this frame doesn't correspond to a system of coordinates. Otherwise the spacetime would be flat).
Don't forget that vectors are contravariant tensors, so the components of the 4-velocity will be transformed by this matrix:
aut2.inverse()[:]
The angle $\theta$ can now be computed using:
$$\theta = \arctan\left(\frac{\overrightarrow{v}_{\bot}}{\overrightarrow{v}_{\parallel}}\right)$$
Because this formula is using 3-vectors, it was indeed important to have been in an orthonormal frame.
We also need to perform the rotation of angle $\alpha$ and $\beta$ for the 4-velocity. In the next cell is a Cython version optimized for speed.
%%cython
from libc.math cimport cos, sin, sqrt
cpdef tuple spherical_to_xyz(float dr, float dtheta, float dphi, float r,
float theta, float phi, float alpha, float beta):
"""
Convert spherical coordinates to cartesian and apply the
two rotations at the same time.
"""
cdef float dx, dy, dz
dx = (-cos(beta)*cos(theta)*sin(alpha)-(sin(alpha)*sin(beta)*sin(phi)-
cos(alpha)*cos(phi))*sin(theta))*dr+\
(r*cos(beta)*sin(alpha)*sin(theta)-(sin(alpha)*sin(beta)*sin(phi)-\
cos(alpha)*cos(phi))*r*cos(theta))*dtheta+\
(-(cos(phi)*sin(alpha)*sin(beta)+cos(alpha)*sin(phi))*r*sin(theta))*dphi
dy = (cos(beta)*sin(phi)*sin(theta)-sin(beta)*cos(theta))*dr+\
(r*cos(theta)*cos(beta)*sin(phi)+r*sin(beta)*sin(theta))*dtheta+\
(r*cos(phi)*cos(beta)*sin(theta))*dphi
dz = (cos(alpha)*cos(beta)*cos(theta)+(cos(alpha)*sin(beta)*sin(phi)+\
cos(phi)*sin(alpha))*sin(theta))*dr+\
(-r*cos(alpha)*cos(beta)*sin(theta)+(cos(alpha)*sin(beta)*sin(phi)+\
cos(phi)*sin(alpha))*r*cos(theta))*dtheta+\
((cos(alpha)*cos(phi)*sin(beta)-sin(alpha)*sin(phi))*r*sin(theta))*dphi
return (dx, dy, dz)
cpdef tuple xyz_to_spherical(float dx, float dy, float dz, float x,
float y, float z):
"""
Convert cartesian back to spherical
"""
cdef r, dr, dth, dph
r = sqrt(x**2+y**2+z**2)
dr = (x*dx+y*dy*z*dz)/r
dth = ((x*z*dx+y*z*dy)/r**2/sqrt(x**2+y**2)-sqrt(x**2+y**2)*dz)/r**2
dph = -y/(x**2+y**2)*dx+x/(x**2+y**2)*dy
return (dr, dth, dph)
These formulas were obtained automatically by creating another chart and asking for the change of frames. Example:
def print_formulas(): # enclosed in a function to prevent altering the namespace
alpha, beta = var('alpha, beta')
spher.<r, theta, phi> = E.chart()
x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta) # normal Spherical->Cartesian transformation
y, z = y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta) # first rotation
x, z = x*cos(alpha)-z*sin(alpha), z*cos(alpha)+x*sin(alpha) # second rotation
spher.transition_map(E.default_chart(), [x, y, z])
print(list(E.changes_of_frame().values())[0][spher.frame(),:, spher])
print_formulas()
[ -cos(beta)*cos(theta)*sin(alpha) - (sin(alpha)*sin(beta)*sin(phi) - cos(alpha)*cos(phi))*sin(theta) r*cos(beta)*sin(alpha)*sin(theta) - (sin(alpha)*sin(beta)*sin(phi) - cos(alpha)*cos(phi))*r*cos(theta) -(cos(phi)*sin(alpha)*sin(beta) + cos(alpha)*sin(phi))*r*sin(theta)] [ cos(beta)*sin(phi)*sin(theta) - cos(theta)*sin(beta) r*cos(beta)*cos(theta)*sin(phi) + r*sin(beta)*sin(theta) r*cos(beta)*cos(phi)*sin(theta)] [ cos(alpha)*cos(beta)*cos(theta) + (cos(alpha)*sin(beta)*sin(phi) + cos(phi)*sin(alpha))*sin(theta) -r*cos(alpha)*cos(beta)*sin(theta) + (cos(alpha)*sin(beta)*sin(phi) + cos(phi)*sin(alpha))*r*cos(theta) (cos(alpha)*cos(phi)*sin(beta) - sin(alpha)*sin(phi))*r*sin(theta)]
The thickness of the disk will follow this profile, obtained from my high level understanding of black holes mechanics (i.e. it's mostly random, but looks nice)
%%cython
from libc.math cimport exp, erf
cpdef float profile(float x, float disk_min, float disk_max):
cdef float y
# we really don't want negative values
if x<disk_min or x>disk_max:
return 0
y = exp(-(disk_min-20-x)**2/400)*(x-disk_min)**2*(disk_max-x)**2/10000 +\
exp(-(32-x)**2/70)/2*(x-disk_min)**2*(disk_max-x)**2/150000
return max(y, 0)
plot(lambda x, d1=disk_min, d2=disk_max: profile(x, d1, d2),
xmin=0, xmax=60, ymin=0, ymax=1.2)
We can know rewrite the function intersection
to take everything into account. This time it directly returns an RGB value.
%%cython
from libc.math cimport cos, sin, acos, sqrt, abs, atan2
cimport cython
from __main__ import profile
from __main__ import xyz_to_spherical
from __main__ import spherical_to_xyz
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef tuple intersection(list curve, float m, float alpha, float beta,
float dmin, float dmax):
cdef float x, y, z
cdef float x2, y2, z2
cdef float r, theta, phi
cdef int n, i
cdef float t
cdef float sinalpha, cosalpha
cdef float sinbeta, cosbeta
cdef float R, G, B
cdef float dr, dtheta, dphi
cdef float dx, dy, dz
cdef float th
R, G, B = 0., 0., 0. # return values
sinalpha = sin(alpha)
cosalpha = cos(alpha)
sinbeta = sin(beta)
cosbeta = cos(beta)
n = len(curve)
r, theta, phi = curve[0][2:5]
x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)
z = z*cos(alpha)+x*sin(alpha)
for i in range(1, n):
r = curve[i][2]
theta = curve[i][3]
phi = curve[i][4]
# rotations
x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
y2, z2 = y2*cosbeta-z2*sinbeta, z2*cosbeta+y2*sinbeta
x2, z2 = x2*cosalpha-z2*sinalpha, z2*cosalpha+x2*sinalpha
if z!=z2:
t = z/(z-z2)
if t>=0 and t<1 and curve[i][2]>dmin and curve[i][2]<dmax:
# 4-velocity
dt = curve[i+1][1]-curve[i][1]
dr = curve[i+1][2]-curve[i][2]
dtheta = curve[i+1][3]-curve[i][3]
dphi = curve[i+1][4]-curve[i][4]
# rotation of the 4-velocity
dx, dy, dz = spherical_to_xyz(dr, dtheta, dphi, r, theta, phi,
alpha, beta)
dr, dtheta, dphi = xyz_to_spherical(dx, dy, dz, x2, y2, z2)
# change of frame to orthonorma frame:
dr = 2*m/(sqrt(r*(2*m+r)))*dt+sqrt((2*m+r)/r)*dr
dtheta = dtheta*r
dphi = dphi*r # theta = pi/2
dt = sqrt(r/(2*m+r))*dt
# angle theta
th = atan2(dtheta, (sqrt(dr**2+dphi**2)))
# thickness
R += 40*profile(r, dmin, dmax)/sin(abs(th)) # arbitrary factor
x, y, z = x2, y2, z2
# clipping
if R >= 255:
R = 255
return R, G, B
We also have to rewrite render_row
to accept an RGB value.
def render_row(x):
"""
Render a single row of the image
"""
res = np.zeros((ny,3)) # result row in RGB format
for y in range(ny):
beta = atan2(y-ny/2,x) # beta angle
r = sqrt(x**2+(y-ny/2)**2) # pixel distance to the center of the image
ind_geo = int(r/400*n_geod*720/nx) # index of the geodesic to use. values are obtained by trial and error.
if ind_geo<n_geod: # don't bother if the index is too big.
res[ny-1-y,:] = intersection(geo[ind_geo], 2, float(pi/20), float(beta),
float(disk_min), float(disk_max))
return x, res # also return the argument, which is not remembered otherwise
Let's try it:
data = np.zeros( (ny, nx, 3), dtype=float )
render()
img2 = smp.toimage(data)
img2
FloatProgress(value=0.0, max=720.0)
/home/eric/sage/8.7/local/lib/python2.7/site-packages/sage/repl/ipython_kernel/__main__.py:3: DeprecationWarning: `toimage` is deprecated! `toimage` is deprecated in SciPy 1.0.0, and will be removed in 1.2.0. Use Pillow's ``Image.fromarray`` directly instead. IPKernelApp.launch_instance(kernel_class=SageKernel)
Much more beautiful!
But enough of red black holes, let's add real spectra into the mix.
It's quite easy to associate a spectrum to each light ray, but it's also important to be able to convert it back to displayable colors.
This part shows how to do that, using the CIE standard XYZ function. This function should contain every information about a spectrum that are relevant to the human eye. They can then be converted to RGB for example. What's more, these function are simply obtained by integrating the spectrum against a function. They are the linearly dependant of the spectrum.
To compute the XYZ functions, we first need the function that defines them. I got mine from here :
http://www.cvrl.org/database/data/cmfs/ciexyzjv.csv
But because we are planning for the Doppler effect, I added a lot of zeros at the beginning and the end to encompass a more wavelengths (I went from 5 to 3000 nm, instead of 380-780 for visible light).
from six.moves.urllib.request import urlretrieve # valid for both Python 2 and Python 3
urlretrieve("http://www.cvrl.org/database/data/cmfs/ciexyzjv.csv",
"ciexyzjv.csv")
ciexyz = np.genfromtxt("ciexyzjv.csv", delimiter=",")
Here is what they look like. Y can be seen as the luminosity of a normalized spectrum (e.g. for the same light intensity, green will appear much brighter to the human eye than purple). The other two don't really have an interpretation.
plt.plot(ciexyz[:,0],ciexyz[:,3], label='Z')
plt.plot(ciexyz[:,0],ciexyz[:,2], label='Y')
plt.plot(ciexyz[:,0],ciexyz[:,1], label='X')
plt.legend(loc='best')
plt.show()
Next we need to compute the black body spectrum for a given temperature, using Planck's law.
I already include doppler effect in the formula.
%%cython
from libc.math cimport exp
cpdef float blackbody(float nu, float T, float doppler):
"""
Spectral power emmited at frequency nu by a black-body at
temperature T per square meter par steradian
"""
cdef float h = 6.62e-34
cdef float k = 1.38e-23
cdef float c = 3e8
cdef float h_sur_k = 4.79710144927536e-11
return (2*h)*nu/c*nu/c*nu/(exp(h_sur_k*nu/doppler/T)-1)
In fact there is a small trick in the previous cell. Because of the exponent in the constants, it's fairly easy to overflow. That's why the $\nu^3/c^2$ is written this way.
Let's try at 10,000 K:
plot(spline([(l, blackbody(3e8/(l/1e9), 10000, 1)) for l in ciexyz[:,0]]),
xmin=5, xmax=3000)
verbose 0 (3629: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 170 points. verbose 0 (3629: plot.py, generate_plot_points) Last error message: 'Unable to compute f(3000.0)'
Once again, we'll use cython to convert a temperature into XYZ. Because XYZ depend linearly on the intensity, and not RGB, the conversion will be done in post-processing. Everything is redeclared in Cython for optimization (hence all the type declarations). Even the XYZ arrays are reloaded in this Cython environment.
%%cython
from __main__ import blackbody
from libc.math cimport exp
import numpy as np
cimport numpy as np
cimport cython
DTYPE = np.float
ctypedef np.float_t DTYPE_t
cdef np.float_t[:] cielamb
cdef np.float_t[:] ciex
cdef np.float_t[:] ciey
cdef np.float_t[:] ciez
def init_arrays():
global cielamb, ciex, ciey, ciez
cdef np.ndarray[np.float_t, ndim=2] ciexyz = np.genfromtxt ('ciexyzjv.csv', delimiter=",")
cielamb = ciexyz[:, 0]
ciex = ciexyz[:, 1]
ciey = ciexyz[:, 2]
ciez = ciexyz[:, 3]
init_arrays()
cpdef tuple temp_to_XYZ(float T, float doppler):
global cielamb, ciex, ciey, ciez
cdef int nl = len(cielamb)
cdef np.ndarray[np.float_t, ndim=1] sp = np.zeros(nl)
cdef int i
cdef float x, y, z
for i in range(nl):
sp[i] = blackbody(3e8/(cielamb[i]/1e9), T, doppler)
x = np.dot(sp, ciex)
y = np.dot(sp, ciey)
z = np.dot(sp, ciez)
return (x, y, z)
Let's convert it to RGB, using one of the many formulas, just to see how it looks:
def xyz_to_rgb(*args):
# constants
fact = 3e-8 # arbitrary
gamma = 1/2.2
mat = [[3.24047, -1.53715, -0.498835],
[-0.96256, 1.8752, 0.041556],
[0.055648, -0.204043, 1.057311]]
# conversion
r, g, b = np.dot(mat, np.transpose(args)/fact).tolist()
# gamma correction and clipping
r = min(1,max(0, r)**gamma)
g = min(1,max(0, g)**gamma)
b = min(1,max(0, b)**gamma)
return r, g, b
from matplotlib.patches import Rectangle
axes = plt.gca()
axes.set_xlim([1000, 6000])
for T in range(1000, 6000, 200):
# plot a rectangle at a color obtained from the temperature
axes.add_patch(Rectangle((T, 0), 200, 1, facecolor=xyz_to_rgb(*temp_to_XYZ(T, 1))))
With this formula, it becomes realy bright really soon. Adjusting the $\gamma$ correction could help.
We can now add it to our black hole ray tracer:
%%cython
from libc.math cimport cos, sin, acos, sqrt, abs, atan2
cimport cython
from __main__ import profile
from __main__ import xyz_to_spherical
from __main__ import spherical_to_xyz
from __main__ import temp_to_XYZ
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef tuple intersection(list curve, float m, float alpha, float beta,
float dmin, float dmax):
cdef float x, y, z
cdef float x2, y2, z2
cdef float r, theta, phi
cdef int n, i
cdef float t
cdef float sinalpha, cosalpha
cdef float sinbeta, cosbeta
cdef float X, Y, Z
cdef float X0, Y0, Z0
cdef float dr, dtheta, dphi
cdef float dx, dy, dz
cdef float th, doppler, factor,
X, Y, Z = 0., 0., 0. # return values
# 20 percent speed gain
sinalpha = sin(alpha)
cosalpha = cos(alpha)
sinbeta = sin(beta)
cosbeta = cos(beta)
n = len(curve)
r, theta, phi = curve[0][2:5]
# rotations
x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)
z = z*cos(alpha)+x*sin(alpha)
for i in range(1, n):
r = curve[i][2]
theta = curve[i][3]
phi = curve[i][4]
# rotations
x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
y2, z2 = y2*cosbeta-z2*sinbeta, z2*cosbeta+y2*sinbeta
x2, z2 = x2*cosalpha-z2*sinalpha, z2*cosalpha+x2*sinalpha
if z!=z2:
t = z/(z-z2)
if t>=0 and t<1 and curve[i][2]>dmin and curve[i][2]<dmax:
# 4-velocity
dt = curve[i+1][1]-curve[i][1]
dr = curve[i+1][2]-curve[i][2]
dtheta = curve[i+1][3]-curve[i][3]
dphi = curve[i+1][4]-curve[i][4]
# rotation of the 4-velocity
dx, dy, dz = spherical_to_xyz(dr, dtheta, dphi, r, theta, phi,
alpha, beta)
dr, dtheta, dphi = xyz_to_spherical(dx, dy, dz, x2, y2, z2)
# change of frame to orthonormal frame:
dr = 2*m/(sqrt(r*(2*m+r)))*dt+sqrt((2*m+r)/r)*dr
dtheta = dtheta*r
dphi = dphi*r # theta = pi/2
dt = sqrt(r/(2*m+r))*dt
# angle theta
th = atan2(dtheta, (sqrt(dr**2+dphi**2)))
# thickness
factor = profile(r, dmin, dmax)/sin(abs(th))
# unused for now
doppler = 1
# added luminosity
X0, Y0, Z0 = temp_to_XYZ(8000, doppler) # T=8000 because I decided
# updating XYZ
X += X0*factor
Y += Y0*factor
Z += Z0*factor
x, y, z = x2, y2, z2
return X, Y, Z
This time, after the render we still need to convert data
from the XYZ format to RGB.
data = np.zeros( (ny, nx, 3), dtype=float )
render()
FloatProgress(value=0.0, max=720.0)
%%cython
import numpy as np
cimport numpy as np
cpdef np.ndarray[np.int, ndim=3] XYZ_to_RGB(np.ndarray[double, ndim=3] data):
"""
Convert an 2d-array of XYZ values to RGB using a simple formula
"""
cdef int ny = len(data)
cdef int nx = len(data[0])
cdef float x, y, z
cdef int r, g, b
cdef int i, j
cdef np.ndarray[int, ndim=3] res = np.zeros((ny, nx, 3), dtype = np.int32)
cdef float fact = 3e-8
cdef float gamma = 1/1.2
for i in range(ny):
for j in range(nx):
x, y, z = data[i, j, :]
#fact = max(.002*y,1e-7)
x, y, z = x/fact, y/fact, z/fact
r = int(max(0,(3.240479*x -1.537150*y -0.498535*z))**(gamma))
g = int(max(0,(-0.969256*x +1.875992*y +0.041556*z))**(gamma))
b = int(max(0,(0.055648*x -0.204043*y +1.057311*z))**(gamma))
r = min(255, max(0, r))
g = min(255, max(0, g))
b = min(255, max(0, b))
res[i, j, :] = r, g ,b
return res
data_rgb = XYZ_to_RGB(data)
Let's plot the result again:
img3 = smp.toimage(data_rgb)
img3
/home/eric/sage/8.7/local/lib/python2.7/site-packages/sage/repl/ipython_kernel/__main__.py:1: DeprecationWarning: `toimage` is deprecated! `toimage` is deprecated in SciPy 1.0.0, and will be removed in 1.2.0. Use Pillow's ``Image.fromarray`` directly instead. from ipykernel.kernelapp import IPKernelApp
Everything was done to make this part easy. All we have to do is evalute the Doppler factor.
To compute de Doppler effect due to the rotation of the disk, we need to consider the interaction between to disk and the light ray in the orthonormal frame again. The goal is to use the general formula:
$$f_0 = \frac{f_s}{\gamma \left(1+\beta \cos(\theta)\right)}$$To do that, we need to estimate $\beta$ and $\theta$.
Lets first look at the equation of motion of a single particule, and impose circular motion: $dr=0$ and $\theta=\pi/2$
eqs = curve._equations_rhs[C]
dt, dr, dth, dph = curve._velocities
The 4 equations becomes:
# d^2t/dtau^2 =
eqs[0].subs({th: pi/2, dr: 0, dth: 0})
# d^2r/dtau^2 =
eqs[1].subs({th: pi/2, dr: 0, dth: 0})
# d^2th/dtau^2 =
eqs[2].subs({th: pi/2, dr: 0, dth: 0})
# d^2ph/dtau^2 =
eqs[3].subs({th: pi/2, dr: 0, dth: 0})
In particular $d\phi$ is constant. We also solve the second equation for $dt$:
(the sign depend on the rotation directions)
eqs[1].subs({th: pi/2, dr: 0, dth: 0}).solve(dt)
And so by reinjecting in the first equation $dt$ is constant, and so is $d\phi$, like in Newtonian gravity.
Their values are obtained from the normalization condition.
p = M((0, r, pi/2, 0))
Tp = M.tangent_space(p)
v = Tp((r*sqrt(r/m)*dph, 0, 0, dph))
g.at(p)(v,v)
(g.at(p)(v,v)==-1).solve(dph)
As we can see, the smallest circular orbit can only by at $r=3m$, with:
We can see that:
If we switch to the orthonormal frame (just multiply $d\phi$ by r, as seen in the change of frame), we can see that:
$$\beta = \frac{v}{c} = \sqrt{\frac{m}{r-3m}}$$To find $\theta$, we use the formula:
$$\theta = \arccos\left(\frac{\overrightarrow{v_1}\cdot \overrightarrow{v_2}}{||\overrightarrow{v_1}||\cdot||\overrightarrow{v_2}||}\right)$$Finally, the gravitational redshift is easier. It is simply given by the ratio of first component of the 4-velocity in the orthonormal frame between the two end of the geodesic.
%%cython
from libc.math cimport cos, sin, acos, sqrt, abs, atan2
cimport cython
from __main__ import profile
from __main__ import xyz_to_spherical
from __main__ import spherical_to_xyz
from __main__ import temp_to_XYZ
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef tuple intersection(list curve, float m, float alpha, float beta,
float dmin, float dmax):
"""
Return True if the curve intersect the disk comprised between dmin
and dmax tilted of angles alpha and beta
"""
cdef float x, y, z
cdef float x2, y2, z2
cdef float r, theta, phi
cdef int n, i
cdef float t
cdef float sinalpha, cosalpha
cdef float sinbeta, cosbeta
cdef float X, Y, Z
cdef float X0, Y0, Z0
cdef float dt, dt0, dr, dtheta, dphi
cdef float dx, dy, dz
cdef float th, rho, doppler, beta_rel
X, Y, Z = 0., 0., 0. # return values
# 20 percent speed gain
sinalpha = sin(alpha)
cosalpha = cos(alpha)
sinbeta = sin(beta)
cosbeta = cos(beta)
n = len(curve)
r, theta, phi = curve[0][2:5]
dt0 = curve[1][1]-curve[0][1]
# rotations
x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)
z = z*cos(alpha)+x*sin(alpha)
for i in range(1, n):
r = curve[i][2]
theta = curve[i][3]
phi = curve[i][4]
# rotations
x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
y2, z2 = y2*cosbeta-z2*sinbeta, z2*cosbeta+y2*sinbeta
x2, z2 = x2*cosalpha-z2*sinalpha, z2*cosalpha+x2*sinalpha
if z!=z2:
t = z/(z-z2)
if t>=0 and t<1 and curve[i][2]>dmin and curve[i][2]<dmax:
# 4-velocity
dt = curve[i+1][1]-curve[i][1]
dr = curve[i+1][2]-curve[i][2]
dtheta = curve[i+1][3]-curve[i][3]
dphi = curve[i+1][4]-curve[i][4]
# rotation of the 4-velocity
dx, dy, dz = spherical_to_xyz(dr, dtheta, dphi, r, theta, phi,
alpha, beta)
dr, dtheta, dphi = xyz_to_spherical(dx, dy, dz, x2, y2, z2)
# change of frame to orthonormal frame:
dr = 2*m/(sqrt(r*(2*m+r)))*dt+sqrt((2*m+r)/r)*dr
dtheta = dtheta*r
dphi = dphi*r # theta = pi/2
dt = sqrt(r/(2*m+r))*dt
# angle theta (used for thickness)
th = atan2(dtheta, (sqrt(dr**2+dphi**2)))
# thickness
factor = profile(r, dmin, dmax)/sin(abs(th))
# angle rho (used for Doppler effect)
rho = acos(dphi/sqrt(dr**2+dtheta**2+dphi**2))
# beta (used for Doppler effect)
beta_rel = sqrt(m/(r-3*m))
# Doppler factor
doppler = 1
doppler *= dt0/dt # gravitational redshift
doppler *= 1/(sqrt(1-beta_rel**2)*(1+beta_rel*cos(rho))) # rotation of the disk
# added luminosity
X0, Y0, Z0 = temp_to_XYZ(8000, doppler) # T=8000 because I decided
# updating XYZ
X += X0*factor
Y += Y0*factor
Z += Z0*factor
x, y, z = x2, y2, z2
return X, Y, Z
When plotted, this gives a nice asymetry.
data = np.zeros( (ny, nx, 3), dtype=float )
render()
data_rgb = XYZ_to_RGB(data)
img4 = smp.toimage(data_rgb)
img4
FloatProgress(value=0.0, max=720.0)
/home/eric/sage/8.7/local/lib/python2.7/site-packages/sage/repl/ipython_kernel/__main__.py:4: DeprecationWarning: `toimage` is deprecated! `toimage` is deprecated in SciPy 1.0.0, and will be removed in 1.2.0. Use Pillow's ``Image.fromarray`` directly instead.
One last thing we didn't take into account is that a moving object that emit isotropic light at rest doesn't do so when in motion. The classical aberration formula reads:
$$\theta' = \arccos \left( \frac{\cos \theta -\beta}{1-\beta \cos \theta} \right)$$In the rest frame, the number of light rays comming out from the object at an angle $\theta'$ is constant, equal to $\frac{dN}{d\theta'}$.
By using chain derivation, in the moving frame :
$$\frac{dN}{d\theta} = \frac{dN}{d\theta'} \frac{d\theta'}{d\theta}$$As we can see, the number or ray per angle unit is multiplied by $\frac{d\theta'}{d\theta}$. Let's plot this function for multiple values of $\beta$.
beta = var('beta')
thp = acos((cos(th)-beta)/(1-beta*cos(th)))
diff(thp,th)
aberration = abs(diff(thp,th)) # abs needed if we want negative angles
P = plot(aberration.subs(beta=0.), [-pi*.99,pi*.99])
for i in range(1,10):
P += plot(aberration.subs({beta: 0.1*i}), [-pi*.99,pi*.99])
P.plot()
As we can see, the effect does focuses rays toward the direction of motion. One can also check that the integral does not depend on $\beta$.
This is a fast cython implementation:
%%cython
from libc.math cimport cos, sin ,sqrt
cpdef float focalisation_factor(float th, float beta):
# uses an other expression stricly equal, but faster to compute.
return (cos(th)*(cos(th)+beta)+sin(th)**2)/\
(sin(th)**2*sqrt(1-beta**2)+1/sqrt(1-beta**2)*(beta+cos(th))**2)
Finally, this is the final code for the intersection:
%%cython
from libc.math cimport cos, sin, acos, sqrt, abs, atan2
cimport cython
from __main__ import profile
from __main__ import xyz_to_spherical
from __main__ import spherical_to_xyz
from __main__ import temp_to_XYZ
from __main__ import focalisation_factor
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef tuple intersection(list curve, float m, float alpha, float beta,
float dmin, float dmax):
cdef float x, y, z
cdef float x2, y2, z2
cdef float r, theta, phi
cdef int n, i
cdef float t
cdef float sinalpha, cosalpha
cdef float sinbeta, cosbeta
cdef float X, Y, Z
cdef float X0, Y0, Z0
cdef float dt, dt0, dr, dtheta, dphi
cdef float dx, dy, dz
cdef float th, rho, doppler, beta_rel
X, Y, Z = 0., 0., 0. # return values
# 20 percent speed gain
sinalpha = sin(alpha)
cosalpha = cos(alpha)
sinbeta = sin(beta)
cosbeta = cos(beta)
n = len(curve)
r, theta, phi = curve[0][2:5]
dt0 = curve[1][1]-curve[0][1]
# rotations
x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)
z = z*cos(alpha)+x*sin(alpha)
for i in range(1, n):
r = curve[i][2]
theta = curve[i][3]
phi = curve[i][4]
# rotations
x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)
y2, z2 = y2*cosbeta-z2*sinbeta, z2*cosbeta+y2*sinbeta
x2, z2 = x2*cosalpha-z2*sinalpha, z2*cosalpha+x2*sinalpha
if z!=z2:
t = z/(z-z2)
if t>=0 and t<1 and curve[i][2]>dmin and curve[i][2]<dmax:
# 4-velocity
dt = curve[i+1][1]-curve[i][1]
dr = curve[i+1][2]-curve[i][2]
dtheta = curve[i+1][3]-curve[i][3]
dphi = curve[i+1][4]-curve[i][4]
# rotation of the 4-velocity
dx, dy, dz = spherical_to_xyz(dr, dtheta, dphi, r, theta, phi,
alpha, beta)
dr, dtheta, dphi = xyz_to_spherical(dx, dy, dz, x2, y2, z2)
# change of frame to orthonormal frame:
dr = 2*m/(sqrt(r*(2*m+r)))*dt+sqrt((2*m+r)/r)*dr
dtheta = dtheta*r
dphi = dphi*r # theta = pi/2
dt = sqrt(r/(2*m+r))*dt
# angle theta (used for thickness)
th = atan2(dtheta, (sqrt(dr**2+dphi**2)))
# thickness
factor = profile(r, dmin, dmax)/sin(abs(th))
# angle rho (used for Doppler effect)
rho = acos(dphi/sqrt(dr**2+dtheta**2+dphi**2))
# beta (used for Doppler effect and focalisation)
beta_rel = sqrt(m/(r-3*m))
# Doppler factor
doppler = 1
doppler *= dt0/dt # gravitational redshift
doppler *= 1/(sqrt(1-beta_rel**2)*(1+beta_rel*cos(rho))) # rotation of the disk
# added luminosity
X0, Y0, Z0 = temp_to_XYZ(8000, doppler) # T=8000 because I decided
# focalisation
factor *= focalisation_factor(rho, beta_rel)
# updating XYZ
X += X0*factor
Y += Y0*factor
Z += Z0*factor
x, y, z = x2, y2, z2
return X, Y, Z
And this is the final rendering:
data = np.zeros( (ny, nx, 3), dtype=float )
render()
data_rgb = XYZ_to_RGB(data)
img5 = smp.toimage(data_rgb)
# img5.save("my home made black hole","png") uncomment to save
img5
FloatProgress(value=0.0, max=720.0)
/home/eric/sage/8.7/local/lib/python2.7/site-packages/sage/repl/ipython_kernel/__main__.py:4: DeprecationWarning: `toimage` is deprecated! `toimage` is deprecated in SciPy 1.0.0, and will be removed in 1.2.0. Use Pillow's ``Image.fromarray`` directly instead.
The difference with the previous image seems small, but that's only because the transition to RGB intoduces clipping. In fact the white area is twice as bright on the new picture.
img1
img2
img3
img4
img5