This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.
version()
'SageMath version 9.1.beta8, Release Date: 2020-03-18'
%display latex
We evaluate $\Phi_b(u)$ via the functions elliptic_f and elliptic_kc of SageMath, taking into account that the second argument of elliptic_f
and the argument of elliptic_kc
is $m = k^2$ and not $k$.
def Phi(u, b):
u = RDF(u)
b = RDF(b)
bc = RDF(3*sqrt(3))
if b > bc:
xi = 2*arcsin(bc/b)/3
un = RDF(1/3*cos(xi + 2*pi/3) + 1/6)
up = RDF(1/3*cos(xi + 4*pi/3) + 1/6)
ua = RDF(1/3*cos(xi) + 1/6)
k2 = (up - un)/(ua - un)
phi = arcsin(sqrt((u - un)/(up - un)))
# NB: elliptic_kc() and elliptic_f() takes m=k^2 as argument
aa = elliptic_kc(k2) - elliptic_f(phi, k2)
return sqrt(2/(ua - un)) * aa
else:
xi = bc/b - sqrt((bc/b)^2 - 1)
un = (1 - xi^(2/3) - xi^(-2/3))/6
us = sqrt(un*(3*un - 1)) + un
k2 = (us - 2.5*un + 0.25)/(2*(us - un))
phi = arccos(abs(us - u)/(us + u - 2*un))
# NB: elliptic_kc() take m=k^2 as argument
aa = elliptic_kc(k2) - elliptic_f(phi, k2)
if u > us:
aa = - aa
return aa / sqrt(2*(us - un))
def phi_infinity(r_em, b, phi_em=0, pre_periastron=True, eps_L=1):
# print("r_em: {}, b: {}".format(r_em, b))
u_em = RDF(1)/RDF(r_em)
b = RDF(b)
bc = RDF(3*sqrt(3))
phi_em = RDF(phi_em)
if b > bc:
if pre_periastron:
return phi_em + eps_L*(Phi(0, b) + Phi(u_em, b))
return phi_em + eps_L*(Phi(0, b) - Phi(u_em, b))
return phi_em + eps_L*(Phi(0, b) - Phi(u_em, b))
phi_infinity(6., 5.2, phi_em=pi/2, pre_periastron=False)
def find_b(phi_inf, r_em, phi_em, pre_periastron=True):
phi_inf = RDF(phi_inf)
r_em = RDF(r_em)
phi_em = RDF(phi_em)
bmax = 0.999*sqrt(r_em**3 / (r_em - RDF(2)))
eps_L = 1 if phi_inf > phi_em else -1
bc1 = 1.0000001*RDF(3*sqrt(3))
def delta_phi(b):
return phi_infinity(r_em, b, phi_em=phi_em,
pre_periastron=pre_periastron,
eps_L=eps_L) - phi_inf
return find_root(delta_phi, bc1, bmax)
find_b(0, 6., pi/2, pre_periastron=False)
find_b(2*pi, 6., pi/2, pre_periastron=True)
find_b(-2*pi, 6., pi/2, pre_periastron=True)
def phi_sol(r, r_em, phi_em, b, em_pre_periastron=True,
pre_periastron=True, eps_L=1):
u = RDF(1)/RDF(r)
u_em = RDF(1)/RDF(r_em)
phi_em = RDF(phi_em)
if em_pre_periastron:
if pre_periastron:
return phi_em + eps_L*(Phi(u_em, b) - Phi(u, b))
return phi_em + eps_L*(Phi(u_em, b) + Phi(u, b))
return phi_em + eps_L*(Phi(u, b) - Phi(u_em, b))
phi_sol(7, 6, pi/2, 6.93)
def r_periastron(b):
return RDF(2/sqrt(3)*b*cos(pi/3 - arccos(3*sqrt(3)/b)/3))
r_periastron(3*sqrt(3))
r_periastron(6), r_periastron(10)
def geod(r_em, phi_em, b, eps_L, em_pre_periastron=True, rmax=25, np=500):
resu = []
if em_pre_periastron:
rp = 1.001*r_periastron(b)
if r_em < rp:
raise ValueError("r_em < r_per(b)")
dr = RDF(r_em - rp)/RDF(np - 1)
for i in range(np):
r = r_em - i*dr
phi = phi_sol(r, r_em, phi_em, b, em_pre_periastron=True,
pre_periastron=True, eps_L=eps_L)
resu.append((RDF(r)*cos(phi), RDF(r)*sin(phi)))
#
np = 2*np
dr = RDF(rmax - rp)/RDF(np - 1)
for i in range(np):
r = rp + i*dr
phi = phi_sol(r, r_em, phi_em, b, em_pre_periastron=True,
pre_periastron=False, eps_L=eps_L)
resu.append((RDF(r)*cos(phi), RDF(r)*sin(phi)))
else:
dr = RDF(rmax - r_em)/RDF(np - 1)
for i in range(np):
r = r_em + i*dr
phi = phi_sol(r, r_em, phi_em, b, em_pre_periastron=False,
pre_periastron=False, eps_L=eps_L)
resu.append((RDF(r)*cos(phi), RDF(r)*sin(phi)))
return resu
Inner and outer radii of the accretion disk:
r_disk_in = 6
r_disk_out = 9
Colors of the various images:
color1 = 'orange'
color2 = 'peru'
color3 = 'red'
phi_em = pi/2
r_em = r_disk_in
em_pre_periastron = False
b1_in = find_b(0, r_em, phi_em, pre_periastron=em_pre_periastron)
b1_in
eps_L = -1
graph = circle((0, 0), 2, edgecolor='black', fill=True, facecolor='grey', alpha=0.5)
graph += circle((0, 0), 3, color='lightgreen', thickness=1., linestyle='--', zorder=1)
graph += line([(0, r_disk_in), (0, r_disk_out)], color='darkorange',
thickness=4, zorder=50)
graph += line(geod(r_em, phi_em, b1_in, eps_L, em_pre_periastron=em_pre_periastron),
color=color1)
graph
em_pre_periastron = True
b2_in = find_b(2*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
b2_in
eps_L = 1
graph += line(geod(r_em, phi_em, b2_in, eps_L, em_pre_periastron=em_pre_periastron),
color=color2)
graph
eps_L = -1
b3_in = find_b(-2*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
b3_in
graph += line(geod(r_em, phi_em, b3_in, eps_L, em_pre_periastron=em_pre_periastron,
np=4000), color=color3)
graph
r_em = r_disk_out
em_pre_periastron = False
b1_out = find_b(0, r_em, phi_em, pre_periastron=em_pre_periastron)
eps_L = -1
graph += line(geod(r_em, phi_em, b1_out, eps_L, em_pre_periastron=em_pre_periastron),
color=color1)
em_pre_periastron = True
eps_L = 1
b2_out = find_b(2*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
graph += line(geod(r_em, phi_em, b2_out, eps_L, em_pre_periastron=em_pre_periastron),
color=color2)
eps_L = -1
b3_out = find_b(-2*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
graph += line(geod(r_em, phi_em, b3_out, eps_L, em_pre_periastron=em_pre_periastron,
np=4000), color=color3)
show(graph, xmax=15, axes_labels=[r'$x/m$', r'$y/m$'])
graph.save("gis_disk_image_detail.pdf", xmax=15,
axes_labels=[r'$x/m$', r'$y/m$'])
Inner radius and width of primary image:
b1_in, b1_out - b1_in
Inner radius and width of secondary image:
b2_in, b2_out - b2_in
Inner radius and width of tertiary image:
b3_in, b3_out - b3_in
phi_em = 3*pi/2
r_em = r_disk_in
em_pre_periastron = False
b = find_b(2*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
eps_L = 1
graph += line([(0, -r_disk_in), (0, -r_disk_out)], color='darkorange',
thickness=4, zorder=50)
graph += line(geod(r_em, phi_em, b, eps_L, em_pre_periastron=em_pre_periastron),
color=color1)
em_pre_periastron = True
eps_L = -1
b = find_b(0, r_em, phi_em, pre_periastron=em_pre_periastron)
graph += line(geod(r_em, phi_em, b, eps_L, em_pre_periastron=em_pre_periastron),
color=color2)
eps_L = 1
b = find_b(4*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
graph += line(geod(r_em, phi_em, b, eps_L, em_pre_periastron=em_pre_periastron,
np=1000), color=color3)
r_em = r_disk_out
em_pre_periastron = False
b = find_b(2*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
eps_L = 1
graph += line(geod(r_em, phi_em, b, eps_L, em_pre_periastron=em_pre_periastron),
color=color1)
em_pre_periastron = True
eps_L = -1
b = find_b(0, r_em, phi_em, pre_periastron=em_pre_periastron)
graph += line(geod(r_em, phi_em, b, eps_L, em_pre_periastron=em_pre_periastron),
color=color2)
eps_L = 1
b = find_b(4*pi, r_em, phi_em, pre_periastron=em_pre_periastron)
graph += line(geod(r_em, phi_em, b, eps_L, em_pre_periastron=em_pre_periastron,
np=1000), color=color3)
graph
show(graph, xmax=20)
graph1 = graph + polygon([(15,-10), (15,10), (25,10), (25,-10)], color='white',
zorder=100, axes=False)
graph1
graph1 += circle((25, 0), b1_out, color=color1, fill=True, zorder=101)
graph1 += circle((25, 0), b1_in, color='white', fill=True, thickness=0.1, zorder=102)
graph1 += circle((25, 0), b2_out, color=color2, fill=True, zorder=103)
graph1 += circle((25, 0), b2_in, color='white', fill=True, thickness=0.1, zorder=104)
graph1 += circle((25, 0), b3_out, color=color3, zorder=105)
graph1
Slight reajustment to take into account that the null geodesics are cut at finite $r$:
a = 0.995
b1_in *= a
b1_out *= a
b2_in *= a
b2_out *= a
b3_in *= a
b3_out *= a
graph1 += line([(15, b1_out), (25, b1_out)], linestyle='dashed',
color=color1, zorder=106)
graph1 += line([(15, b1_in), (25, b1_in)], linestyle='dashed',
color=color1, zorder=107)
graph1 += line([(15, b2_out), (25, b2_out)], linestyle='dashed',
color=color2, zorder=108)
graph1 += line([(15, b2_in), (25, b2_in)], linestyle='dashed',
color=color2, zorder=109)
graph1 += line([(15, b3_out), (25, b3_out)], linestyle='dashed',
color=color3, zorder=110)
graph1 += line([(15, -b1_out), (25, -b1_out)], linestyle='dashed',
color=color1, zorder=106)
graph1 += line([(15, -b1_in), (25, -b1_in)], linestyle='dashed',
color=color1, zorder=107)
graph1 += line([(15, -b2_out), (25, -b2_out)], linestyle='dashed',
color=color2, zorder=108)
graph1 += line([(15, -b2_in), (25, -b2_in)], linestyle='dashed',
color=color2, zorder=109)
graph1 += line([(15, -b3_out), (25, -b3_out)], linestyle='dashed',
color=color3, zorder=110)
graph1
graph1.save("gis_disk_image.pdf")