This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.
Click here to download the notebook file (ipynb format). To run it, you must start SageMath with sage -n jupyter
.
A version of Sage at least equal to 9.0 is required to run this notebook:
version()
'SageMath version 9.1.beta0, Release Date: 2020-01-10'
%display latex
M = Manifold(4, 'M', structure='Lorentzian')
X.<t, r, th, ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):periodic:\phi')
X
g = M.metric()
m = 1
g[0,0] = -(1-2*m/r)
g[1,1] = 1/(1-2*m/r)
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()
Xcart.<t, x, y, z> = M.chart()
X_to_Xcart = X.transition_map(Xcart, [t, r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
X_to_Xcart.display()
M.top_charts()
M._top_charts = [X]
M.top_charts()
M.identity_map().coord_functions(X, Xcart)
def initial_vector(r0, b, ph0=0, E=1, inward=False):
t0, th0 = 0, pi/2
L = b*E
vt0 = 1/(1-2*m/r0)
vr0 = sqrt(E^2 - L^2/r0^2*(1-2*m/r0))
if inward:
vr0 = - vr0
vth0 = 0
vph0 = L / r0^2
p0 = M((t0, r0, th0, ph0), name='p_0')
return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')
r0 = 3
ph0 = 0
b = 3*sqrt(3)
E = 1
print("b = {} m".format(float(b)))
b = 5.196152422706632 m
v0 = initial_vector(r0, b)
v0.display()
p0 = v0.parent().base_point()
g.at(p0)(v0, v0)
s = var('s')
geod = M.integrated_geodesic(g, (s, 0, 100), v0)
geod
sol = geod.solve(step=0.01, method="ode_int") # numerical integration
interp = geod.interpolate() # interpolation of the solution for the plot
graph = circle((0, 0), 2*m, edgecolor='black', fill=True, facecolor='grey', alpha=0.5)
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='green')
graph
U(r) = (1 - 2/r)/r^2
U1 = 0.02
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 10
v0 = initial_vector(r0, b, inward=True)
v0.display()
b = 7.07106781186548 m
geod = M.integrated_geodesic(g, (s, 0, 50), v0)
sol = geod.solve(step=0.01, method="ode_int")
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='red', thickness=1.5,
display_tangent=True, color_tangent='red',
plot_points_tangent=4, scale=1)
graph
U1 = 0.04
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 10
v0 = initial_vector(r0, b, inward=True)
v0.display()
b = 5.00000000000000 m
p0 = v0.parent().base_point()
g.at(p0)(v0, v0)
geod = M.integrated_geodesic(g, (s, 0, 13.25), v0)
sol = geod.solve(step=0.001, method="ode_int")
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='brown', thickness=1.5,
display_tangent=True, color_tangent='brown',
plot_points_tangent=3, scale=0.2)
graph
show(graph, xmin=-10)
U1 = 0.046
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 2.1
v0 = initial_vector(r0, b)
v0.display()
b = 4.66252404120157 m
p0 = v0.parent().base_point()
g.at(p0)(v0, v0)
geod = M.integrated_geodesic(g, (s, 0, 20), v0)
sol = geod.solve(step=0.001, method="ode_int")
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='orange', thickness=1.5,
display_tangent=True, color_tangent='orange',
plot_points_tangent=3, scale=0.1)
graph
show(graph, xmin=-10)
U1 = 0.035
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 2.4
v0 = initial_vector(r0, b, ph0=-pi/2)
v0.display()
b = 5.34522483824849 m
geod = M.integrated_geodesic(g, (s, 0, 3.5), v0)
sol = geod.solve(step=0.001, method="ode_int")
interp = geod.interpolate()
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='grey', thickness=1.5,
display_tangent=True, color_tangent='grey',
plot_points_tangent=3, scale=0.1)
graph
graph += text('1', (5, 4.5), color='red' ) \
+ text('2', (3, 3.8), color='brown' ) \
+ text('3', (-6, 3.5), color='orange' ) \
+ text('4', (0.7, -2.7), color='grey' )
show(graph, xmin=-10, ymin=-4, axes_labels=[r'$x/m$', r'$y/m$'])
graph.save("ges_null_geod.pdf", xmin=-10, ymin=-4, axes_labels=[r'$x/m$', r'$y/m$'])
r0 = 100 # distance of the source
b_sel = [12., 8., 6., 5.355, 5.23, 5.2025, 5.19643, 5.196155]
#b_sel = [15., 10., 7., 6., 5.5, 5.23, 5.22, 5.2025, 5.198, 5.197, 5.1965, 5.19632]
b_pec = [5.355, 5.2025, 5.19632, 5.196155]
#b_sel = [5.19643]
for b in b_sel:
print('b = {:.6f} m'.format(float(b)))
graph = circle((0, 0), 2*m, edgecolor='black', fill=True, facecolor='grey', alpha=0.5)
graph += circle((0, 0), 3*m, color='lightgreen', thickness=1., linestyle='--', zorder=1)
graph += text(r'$b = {:.6f} \, m$'.format(float(b)), (-8, 12), fontsize=14,
color='black')
ph0 = asin(b/r0)
v0 = initial_vector(r0, b, ph0=ph0, inward=True)
geod = M.integrated_geodesic(g, (s, 0, 150), v0)
sol = geod.solve(step=0.001, method="odeint", rtol=1.e-12, atol=1.e-12)
interp = geod.interpolate()
first_point = geod(0)
print("s=0: {}".format(first_point.coord()))
last_point = geod(150)
print("s=150: {}".format(last_point.coord()))
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500,
aspect_ratio=1, color='green', thickness=1.5, zorder=2,
display_tangent=True, color_tangent='green',
plot_points_tangent=12, scale=0.1)
show(graph, xmin=-12, xmax=12, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'])
filename = 'ges_null_b_{:.6f}'.format(float(b))
filename = filename.replace('.', '_') + '.pdf'
graph.save(filename, xmin=-12, xmax=12, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'])
b = 12.000000 m s=0: (0.0, 100.0, 1.5707963267948966, 0.12028988239478806) s=150: (161.7386305994338, 49.599207321700575, 1.5707963267948966, 3.3499140843834243)
b = 8.000000 m s=0: (0.0, 100.0, 1.5707963267948966, 0.08008558003365901) s=150: (165.0262559276518, 47.87567237981162, 1.5707963267948966, 3.832461957010324)
b = 6.000000 m s=0: (0.0, 100.0, 1.5707963267948966, 0.06003605844527842) s=150: (169.51340790400954, 46.21354032584246, 1.5707963267948966, 4.730793903988559)
b = 5.355000 m s=0: (0.0, 100.0, 1.5707963267948966, 0.05357562643499787) s=150: (175.13467942245063, 43.837643428812676, 1.5707963267948966, 6.172293011040655)
b = 5.230000 m s=0: (0.0, 100.0, 1.5707963267948966, 0.052323872006442895) s=150: (180.39625091484754, 41.32000814165003, 1.5707963267948966, 7.665165130714838)
b = 5.202500 m s=0: (0.0, 100.0, 1.5707963267948966, 0.052048497112969744) s=150: (186.05222271416335, 38.48872884422706, 1.5707963267948966, 9.317135314376204)
b = 5.196430 m s=0: (0.0, 100.0, 1.5707963267948966, 0.05198771489670833) s=150: (196.57550821738923, 33.13574656928606, 1.5707963267948966, 12.42140053051524)
b = 5.196155 m s=0: (0.0, 100.0, 1.5707963267948966, 0.05198496117647258) s=150: (212.1623795944796, 25.165821003889214, 1.5707963267948966, 17.043633460770014)
b_c = n(3*sqrt(3))
a = n(648*(7*sqrt(3) - 12))
b = 5.230000
Dphi = 7.665165130714838 - 0.052323872006442895
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)
b = 5.202500
Dphi = 9.317135314376204 - 0.052048497112969744
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)
b = 5.196430
Dphi = 12.42140053051524 - 0.05198771489670833
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)
b = 5.196155
Dphi = 17.043633460770014 - 0.05198496117647258
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)