by Florian Hollants
This Jupyter Notebook simulates a rotating Black Hole in SageMath using the Kerr Metric. The Project is based on the work of Florentin Jaffredo. The parameters for mass and angular momentum can be altered in Cell 5 by changing "m_val" and "a_val", while the "Angle" can be changed in Cell 18.
# n_cpu = 4 # 4 Go Ram minimum
# n_geod = 100
# nx, ny = 180, 90
n_cpu = 8 # 8 Go Ram minimum
n_geod = 1000
nx, ny = 720, 360
# n_cpu = 36 # 144 Go Ram minimum
# n_geod = 30000
# nx, ny = 4000, 2000
%display latex
M = Manifold(4, 'M', structure='Lorentzian')
m = var('m')
a = var('a')
m_val = 2
a_val = 1
C.<t, r, th, ph> = M.chart(r't r:(4.2,+oo) th:(0,pi):\theta ph:\phi') #check that minimum radius is bigger than singularity in g[1,1]
C.coord_range()
g = M.metric()
g[0,0] = (2*m*r)/(r^2+(a*cos(th))^2)-1
g[0,3] = -2*m*r*a*sin(th)^2/(r^2+(a*cos(th))^2)
g[1,1] = (r^2+(a*cos(th))^2)/(r^2-2*m*r+a^2)
g[2,2] = r^2+(a*cos(th))^2
g[3,3] = (r^2+a^2+(2*m*r*a^2*sin(th)^2)/(r^2+(a*cos(th))^2))*sin(th)^2
g[:]
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()
Try the following for circular Orbit
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))
tau = var('tau')
curve = M.integrated_geodesic(g, (tau, 0, 5000), v)
sol = curve.solve(step = 1, method="ode_int", parameters_values={m: m_val, a: a_val})
interp = curve.interpolate()
P = curve.plot_integrated(mapping=phi, color="red", thickness=2, plot_points=5000)
P += sage.plot.plot3d.shapes.Sphere(4, color='grey')
P
Try the following for light at different Values for a, aimed directly at the Black Hole
# p = M((0, 20, pi/2, 0))
# Tp = M.tangent_space(p)
# v = Tp((2, -1, 0.000, 0.00))
# v = v / sqrt(-g.at(p)(v, v))
# A=[0,1,2,4] #list of values for a
# Color=[["red"], ["blue"], ["green"],["orange"], ["black"]] #color for different Geodesics
# P = sage.plot.plot3d.shapes.Sphere(4, color='grey')
# tau = var('tau')
# for i in range(4):
# curve = M.integrated_geodesic(g, (tau, 0, 200), initial_tangent_vector=v, across_charts=True)
# curve.solve_across_charts(step=0.2, parameters_values={m:m_val, a:A[i]})
# interp = curve.interpolate()
# P += curve.plot_integrated(mapping=phi, color=Color[i], thickness=2, plot_points=10000, across_charts=True)
# P
import multiprocessing
from ipywidgets import FloatProgress
from IPython.display import display
def chunks(l, n):
"""Yield successive n-sized chunks from l."""
for i in range(0, len(l), n):
yield l[i:i + n]
n_batches_per_cpu = 3
curve = M.integrated_geodesic(g, (tau, 0, 200), v, across_charts=True)
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)
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))
sol = g.at(p)(v, v).solve(dt)
sol
def calc_some_geodesics(args):
"""
Compute nb geodesics starting at index n0
"""
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):
dy = i*0.006/n_geod
v = Tp([sol[0].rhs()(r0=r, y=dy, m=m_val, a=a_val).n(), -1, 0, dy])
curve._initial_tangent_vector = v
curve.solve_across_charts(step=0.2, parameters_values={m:m_val, a:a_val})
res[i] = (p.coord(), curve._solutions.copy())
curve._solutions.clear()
return res
geo = {}
pool = multiprocessing.Pool(n_cpu)
%display plain
f = FloatProgress(min=0, max=n_geod)
display(f)
for i, some_res in enumerate(pool.map(calc_some_geodesics, args)):
f.value += len(some_res)
geo.update(some_res)
pool.close()
pool.join()
FloatProgress(value=0.0, max=1000.0)
P = sage.plot.plot3d.shapes.Sphere(4, color='grey')
for i in range(0, n_geod, 5*n_geod/100):
curve._solutions = geo[i][1]
interp = curve.interpolate()
P += curve.plot_integrated(mapping=phi, color=["red"], thickness=2, plot_points=150,
label_axes=False, across_charts=True)#
P