This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.
version()
'SageMath version 9.3.beta5, Release Date: 2020-12-27'
%display latex
We use $m=1$ and denote $r_0$ simply by $r$.
a, r = var('a r')
lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph
qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph
The radii of the two horizons:
rp(a) = 1 + sqrt(1 - a^2)
rm(a) = 1 - sqrt(1 - a^2)
rph_ss(a) = 1/2 + cos(2/3*asin(a) + 2*pi/3)
rph_ss
rph_s(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2
rph_s
rph_p(a) = 4*cos(acos(-a)/3)^2
rph_p
rph_m(a) = 4*cos(acos(a)/3)^2
rph_m
We add the radius of the marginally stable orbit:
rph_ms(a) = 1 - (1 - a^2)^(1/3)
rph_ms
as well as the radius of outer and inner polar orbits:
rph_pol(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2)))
rph_pol
rph_pol_in(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2)) + 2*pi/3)
rph_pol_in
a0 = 0.95
show(LatexExpr(r'r_{\rm ph}^{**} = '), n(rph_ss(a0)))
show(LatexExpr(r'r_{\rm ph}^{\rm ms} = '), n(rph_ms(a0)))
show(LatexExpr(r'r_{\rm ph}^{*} = '), n(rph_s(a0)))
show(LatexExpr(r'r_- = '), n(rm(a0)))
show(LatexExpr(r'r_+ = '), n(rp(a0)))
show(LatexExpr(r'r_{\rm ph}^+ = '), n(rph_p(a0)))
show(LatexExpr(r'r_{\rm ph}^- = '), n(rph_m(a0)))
gq = plot(qsph(a0, r), (r, -2, 0.9), thickness=2, legend_label=r'$q$',
axes_labels=[r'$r_0/m$', r'$\ell/m,\ q/m^2$'],
frame=True, gridlines=True) \
+ plot(qsph(a0, r), (r, 1.1, 5), thickness=2)
show(gq, ymin=-4)
gl = plot(lsph(a0, r), (r, -2, 0.99), color='red', thickness=2, legend_label=r'$\ell$') \
+ plot(lsph(a0, r), (r, 1.002, 5), color='red', thickness=2)
show(gq+gl, ymin=-5, ymax=10)
def qmin(a1, r):
l0 = abs(lsph(a1, r))
if l0 < a1:
return -(a1 - l0)^2
return 0
gqmin = plot(lambda r: qmin(a0, r), (r, -2, 5), color='grey', thickness=2, fill=-8)
hor = line([(rp(a0), -10), (rp(a0), 30)], color='black', thickness=2) \
+ line([(rm(a0), -10), (rm(a0), 30)], color='peru', thickness=2)
llim = line([(-2, a0), (5, a0)], color='red', thickness=1.5,
linestyle=':', legend_label=r'$|\ell| = a$') \
+ line([(-2, -a0), (5, -a0)], color='red', thickness=1.5,
linestyle=':')
graph = gq + gl + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph += text(r'$r_{\rm ph}^{**}$', (rph_ss(a0) + 0.1, 2), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^*$', (rph_s(a0) + 0.23, 2), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^+$', (rph_p(a0) + 0.25, 2), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^-$', (rph_m(a0) + 0.2, 2), fontsize=14, zorder=101)
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)
graph.save('gik_spher_orb_exist.pdf', xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)
Zoom:
gq_stable = plot(qsph(a0, r), (r, 0, rph_ms(a0)), thickness=5)
gl_stable = plot(lsph(a0, r), (r, 0, rph_ms(a0)), thickness=5, color='red')
graph = gq + gq_stable + gl + gl_stable + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph += text(r'$r_{\rm ph}^*$', (rph_s(a0) + 0.12, 0.15), fontsize=14, zorder=101)
graph += text(r'$r_{\rm ph}^+$', (rph_p(a0) + 0.09, 0.15), fontsize=14, zorder=101)
r1 = rph_ss(a0)
graph += text(r'$r_{\rm ph}^{**}$', (r1, 0.15), fontsize=14, zorder=101)
graph += line([(r1, 0), (r1, qsph(a0, r1))], linestyle='--', color='blue')
r1 = rph_ms(a0)
graph += text(r'$r_{\rm ph}^{\rm ms}$', (r1, -0.2), color='black',
fontsize=14, zorder=101)
graph += line([(r1, 0), (r1, lsph(a0, r1))], linestyle='--', color='black')
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1, xmax=1.5, ymin=-1.5, ymax=2)
graph.save('gik_spher_orb_exist_zoom.pdf', xmin=-1, xmax=1.5, ymin=-1.5, ymax=2)
a0 = 0.5
gq = plot(qsph(a0, r), (r, -2, 0.9), thickness=2, legend_label=r'$q$',
axes_labels=[r'$r_0/m$', r'$\ell/m,\ q/m^2$'],
frame=True, gridlines=True) \
+ plot(qsph(a0, r), (r, 1.1, 5), thickness=2)
gl = plot(lsph(a0, r), (r, -2, 0.99), color='red', thickness=2, legend_label=r'$\ell$') \
+ plot(lsph(a0, r), (r, 1.002, 5), color='red', thickness=2)
gqmin = plot(lambda r: qmin(a0, r), (r, -2, 5), color='grey', thickness=2, fill=-8)
hor = line([(rp(a0), -10), (rp(a0), 30)], color='black', thickness=2) \
+ line([(rm(a0), -10), (rm(a0), 30)], color='peru', thickness=2)
llim = line([(-2, a0), (5, a0)], color='red', thickness=1.5,
linestyle=':', legend_label=r'$|\ell| = a$') \
+ line([(-2, -a0), (5, -a0)], color='red', thickness=1.5,
linestyle=':')
graph = gq + gl + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)
gq_stable = plot(qsph(a0, r), (r, 0, rph_ms(a0)), thickness=5)
gl_stable = plot(lsph(a0, r), (r, 0, rph_ms(a0)), thickness=5, color='red')
graph = gq + gq_stable + gl + gl_stable + gqmin + hor + llim
show(graph, xmin=-1, xmax=1.5, ymin=-0.2, ymax=0.6)
show(graph, xmin=-1, xmax=1.5, ymin=-0.2, ymax=0.02)
show(graph, xmin=-1, xmax=1.5, ymin=-0.01, ymax=0.01)
a0 = 0.998
gq = plot(qsph(a0, r), (r, -2, 0.999), thickness=2, legend_label=r'$q$',
axes_labels=[r'$r_0/m$', r'$\ell/m,\ q/m^2$'],
frame=True, gridlines=True) \
+ plot(qsph(a0, r), (r, 1.001, 5), thickness=2)
gl = plot(lsph(a0, r), (r, -2, 0.999), color='red', thickness=2, legend_label=r'$\ell$') \
+ plot(lsph(a0, r), (r, 1.001, 5), color='red', thickness=2)
gqmin = plot(lambda r: qmin(a0, r), (r, -2, 5), color='grey', thickness=2, fill=-8)
hor = line([(rp(a0), -10), (rp(a0), 30)], color='black', thickness=2) \
+ line([(rm(a0), -10), (rm(a0), 30)], color='peru', thickness=2)
llim = line([(-2, a0), (5, a0)], color='red', thickness=1.5,
linestyle=':', legend_label=r'$|\ell| = a$') \
+ line([(-2, -a0), (5, -a0)], color='red', thickness=1.5,
linestyle=':')
graph = gq + gl + gqmin + hor + llim
for r1 in [rph_ss(a0), rph_s(a0), rph_p(a0), rph_m(a0)]:
graph += point((r1, qsph(a0, r1)), size=60, color='blue', zorder=100)
graph.set_legend_options(handlelength=2)
show(graph, xmin=-1.5, xmax=4.5, ymin=-7, ymax=28)
gq_stable = plot(qsph(a0, r), (r, 0, rph_ms(a0)), thickness=5)
gl_stable = plot(lsph(a0, r), (r, 0, rph_ms(a0)), thickness=5, color='red')
graph = gq + gq_stable + gl + gl_stable + gqmin + hor + llim
show(graph, xmin=-1, xmax=1.5, ymin=-1.5, ymax=2.5)
P = 2*r^3 - 3*r^2 + a^2
P
a0 = 0.95
plot(P.subs({a: a0}), (r, -0.7, 1.5))
Pdep = (P/2).subs({r: x + 1/2}).simplify_full()
Pdep
pp = -3/4
qq = 1/2*(a^2 - 1/2)
sqrt(-pp/3)
3*qq/(2*pp)
3*qq/(2*pp)*sqrt(-3/pp)
g = Graphics()
for k in range(3):
rk = 1/2 + cos(2/3*asin(a) + 2*k*pi/3)
g += plot(rk, (a, 0, 1), color=hue((3-k)/3), thickness=2)
g
Checking that $x_1 = \cos\left(\frac{2}{3}\arcsin(a) + \frac{2\pi}{3}\right)$ is a solution:
s = Pdep.subs({x: cos(2/3*asin(a) + 2*pi/3)})
s
s.expand_trig().simplify_trig().reduce_trig().expand_trig()
graph = plot(rph_p, (0, 1), axes_labels=[r'$a/m$', r'$r_0/m$'],
color='green', thickness=1.5, legend_label=r'$r_{\rm ph}^+$',
fill=rph_m, fillcolor='palegreen',
frame=True, gridlines=True, axes=True) \
+ plot(rph_m, (0, 1), linestyle='--', color='green', thickness=1.5,
legend_label=r'$r_{\rm ph}^-$') \
+ plot(rph_s, (0, 1), color='green', linestyle=':', thickness=2,
legend_label=r'$r_{\rm ph}^*$', fill=rph_ss, fillcolor='palegreen') \
+ plot(rph_ss, (0, 1), color='green', linestyle='-.', thickness=1.5,
legend_label=r'$r_{\rm ph}^{**}$') \
+ plot(rph_ms, (a, 0, 1), color='blue', linestyle='-', thickness=1.5,
legend_label=r'$r_{\rm ph}^{\rm ms}$', fill=0, fillcolor='darkseagreen') \
+ plot(rp, (0, 1), color='black', thickness=1.5, legend_label=r'$r_+$') \
+ plot(rm, (0, 1), color='darkgoldenrod', thickness=2, legend_label=r'$r_-$')
graph.set_legend_options(handlelength=2, loc=(1.02, 0.36))
show(graph)
graph.save('gik_spher_orb_range.pdf')
graph = plot(lambda a: lsph(a, rph_ss(a)), (0, 1), color='red', thickness=2,
legend_label=r'$\ell_{\rm ph}^{**}$',
axes_labels=[r'$a/m$', r'$\ell_{\rm ph}^{**}/m,\ \ q_{\rm ph}^{**}/m^2$'],
frame=True, gridlines=True, axes=False) \
+ plot(lambda a: qsph(a, rph_ss(a)), (0.001, 0.999), color='blue', thickness=2,
legend_label=r'$q_{\rm ph}^{**}$')
graph.set_legend_options(handlelength=2)
graph.save("gik_ell_q_rss.pdf")
graph
lsph(1, rph_ss(1))
qsph(1, rph_ss(1))
n(_)
th_s = lambda a: arcsin(sqrt(abs(lsph(a, rph_ss(a)))/a))
graph = plot(th_s, (0.0001, 0.9999), color='blue', thickness=2,
axes_labels=[r'$a/m$', r'$\theta^{**}_{\rm ph}$'],
frame=True, gridlines=True, axes=False)
graph.save("gik_theta_ss.pdf")
graph
n(pi/6)
graph = plot(lambda a: lsph(a, rph_s(a)), (0, 1), color='red', thickness=2,
linestyle=':', legend_label=r'$\ell_{\rm ph}^{*}$',
axes_labels=[r'$a/m$', r'$\ell/m$'],
frame=True, gridlines=True, axes=False) \
+ plot(lambda a: lsph(a, rph_p(a)), (0, 1), color='red', thickness=2,
legend_label=r'$\ell_{\rm ph}^+$') \
+ plot(lambda a: lsph(a, rph_m(a)), (0, 1), color='red', thickness=2,
linestyle='--', legend_label=r'$\ell_{\rm ph}^-$')
graph.set_legend_options(handlelength=2)
graph.save("gik_ell_circ_equat.pdf", ymin=-8, ymax=6)
show(graph, ymin=-8, ymax=6)
Generic function $\mathcal{R}(r)$:
l = var('l', latex_name=r'\ell')
q = var('q')
R(a, l, q, r) = r^4 + (a^2 - l^2 - q)*r^2 + 2*(q + (l - a)^2)*r - a^2*q
R
Function $\mathcal{R}(r)$ for a spherical orbit:
r0 = var('r_0')
R0(a, r0, r) = R(a, lsph(a, r0), qsph(a, r0), r).simplify_full()
R0(a, r0, r)
We check that $r_0$ is a double root of $\mathcal{R}(r)$:
R0(a, r0, r).factor()
s = (R0(a, r0, r)*(r0 - 1)^2/(r - r0)^2).simplify_full()
s
bool(s*(r - r0)^2 == R0(a, r0, r).factor().numerator())
s.collect(r)
b = r0^3 - 6*r0^2 + 9*r0 - 4*a^2
b.factor()
bool(R0(a, r0, r) == (r - r0)^2*(r^2 + 2*r0*r + r0*b/(r0 - 1)^2))
Check of the formula $$ \mathcal{R}(r) = (r - r_0)^2 \left( r^2 + 2 r_0 r - a^2 \frac{q}{r_0^2} \right) $$
bool(R0(a, r0, r) == (r - r0)^2*(r^2 + 2*r0*r - a^2/r0^2*qsph(a, r0)))
Rss(a, r) = R0(a, rph_ss(a), r)
Rs(a, r) = R0(a, rph_s(a), r)
Rp(a, r) = R0(a, rph_p(a), r)
Rm(a, r) = R0(a, rph_m(a), r)
Rms(a, r) = R0(a, rph_ms(a), r)
a0 = 0.95
# show(LatexExpr(r'r_{\rm ph}^{**} = ') + latex(n(rph_ss(a0))) + L)
show(LatexExpr(r'r_{\rm ph}^{**} = ' + '{} m'.format(n(rph_ss(a0)))))
plot(lambda r: Rss(a0, r), (-1, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])
A stablle inner orbit:
r0 = 0.4
show(LatexExpr(r'r_0 = ' + '{} m'.format(float(r0))))
plot(lambda r: R0(a0, r0, r), (-0.2, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])
The marginally stable orbit:
show(LatexExpr(r'r_{\rm ph}^{\rm ms} = ' + '{} m'.format(n(rph_ms(a0)))))
plot(lambda r: Rms(a0, r), (-1, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])
show(LatexExpr(r'r_{\rm ph}^{*} = ' + '{} m'.format(n(rph_s(a0)))))
plot(lambda r: Rs(a0, r), (-1, 1), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])
show(LatexExpr(r'r_{\rm ph}^{+} = ' + '{} m'.format(n(rph_p(a0)))))
plot(lambda r: Rp(a0, r), (-1/2, 2), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])
show(LatexExpr(r'r_{\rm ph}^{-} = ' + '{} m'.format(n(rph_m(a0)))))
plot(lambda r: Rm(a0, r), (-1/2, 5), axes_labels=[r'$r/m$', r'$\mathcal{R}(r)$'])
R2(a, l, q, r) = diff(R(a, l, q, r), r, 2).simplify_full()
R2
R2o(a, r) = R2(a, lsph(a, r), qsph(a, r), r).simplify_full().factor()
R2o
dqdr(a, r) = diff(qsph(a, r), r).simplify_full().factor()
dqdr(a, r)
Let us check that marginally stable orbits correspond to an extremum (actually a maximum) of $q(r_0)$:
dqdr(a, rph_ms(a)).simplify_full()
The maximum at $r_0 = 3 m$ does not depend on $a$:
qsph(a, 3)
while for $r_0 = r_{\rm ph}^{ms}$, we have
qms = qsph(a, rph_ms(a)).simplify_full()
qms
q1 = SR(qms._sympy_().simplify())
q1
(qms - q1).simplify_full()
qms = q1
taylor(qms, a, 0, 10)
qms.subs({a: 1})
plot(qms, (a, 0, 1), axes_labels=[r'$a/m$', r'$q_{\rm ms}/m^2$'],
thickness=2, frame=True, gridlines=True)
dldr(a, r) = diff(lsph(a, r), r).simplify_full().factor()
dldr(a, r)
dldr(a, rph_ms(a)).simplify_full()
lms = lsph(a, rph_ms(a)).simplify_full()
lms
taylor(lms, a, 0, 6)
lms.limit(a=1)
plot(lms, (a, 0, 1), color='red', thickness=2,
axes_labels=[r'$a/m$', r'$\ell_{\rm ms}/m$'],
frame=True, gridlines=True) \
+ plot(a, (a, 0, 1), color='black', linestyle='--')
R2(r) = (r + sqrt(r^2 + 4))/2
R2
R4(r) = (r + (r^4 + 16)^(1/4))/2
R4
n(R2(2)), n(R4(2))
n(R2(6)), n(R4(6))
g = plot(R4, (-8, 8), color='green',
axes_labels=[r'$r/m$', r'$\hat{r}/m$'])
show(g, aspect_ratio=1, frame=True, gridlines=True)
g = plot(R2, (-8, 8), axes_labels=[r'$r/m$', r'$\hat{r}/m$']) \
+ plot(r, (r, 0, 8), linestyle='--') \
+ plot(exp(r), (r, -8, 2.1), color='red') \
+ plot(R4, (-8, 8), color='green')
show(g, aspect_ratio=1, frame=True, gridlines=True)
theta0(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 + sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2))))
theta0
theta1(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 - sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2))))
theta1
a0 = 0.95
R = R2
th = var('th')
# We use parametric_plot() instead of arc() because of legend_label
g = parametric_plot((R(0)*sin(th), R(0)*cos(th)),
(th, 0, pi), color='darkorange', thickness=3,
linestyle=':', legend_label=r'$r=0$',
axes_labels=[r'$\frac{\hat{r}}{m}\sin\theta$',
r'$\frac{\hat{r}}{m}\cos\theta$'])
g += parametric_plot((R(rm(a0))*sin(th), R(rm(a0))*cos(th)),
(th, 0, pi), color='peru', thickness=3,
legend_label=r'$r=r_-$')
g += parametric_plot((R(rp(a0))*sin(th), R(rp(a0))*cos(th)),
(th, 0, pi), color='black', thickness=3,
legend_label=r'$r=r_+$')
sing = point((R(0),0), size=60, color='orangered', zorder=100) # ring singularity
rminf = point((0,0), size=60, marker='o', color='white',
markeredgecolor='black', zorder=100)
g += sing + rminf
g.set_legend_options(handlelength=2.5, loc=(0.9, 0.9),
font_size='large')
ns = 9
rmin, rmax = RR(1.00001*rph_p(a0)), RR(0.99999*rph_m(a0))
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2,
plot_points=200, fill=True, fillcolor='palegreen',
fillalpha=0.2)
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
-R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2,
plot_points=200, fill=True, fillcolor='palegreen',
fillalpha=0.2)
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
r0 = rmin + i*dr
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0),
color='green')
rmin, rmax = RR(0), RR(rph_s(a0))
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2,
fill=True, fillcolor='palegreen', fillalpha=0.2)
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
-R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2,
fill=True, fillcolor='palegreen', fillalpha=0.2)
# outer polar orbits:
r0 = rph_pol(a0)
print(lsph(a0, r0))
q0 = qsph(a0, r0)
th0 = 0
print(theta0(a0, 0, q0))
g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0),
linestyle='--', color='green')
ns = 3
rmin, rmax = RR(rph_ms(a0)), 0.9999*RR(rph_s(a0))
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
r0 = rmin + i*dr
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0),
color='green')
ns = 5
rmin, rmax = RR(0), RR(rph_ms(a0))
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
r0 = rmin + i*dr
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0),
color='green')
r0 = rph_ms(a0)
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
g += arc((0, 0), R(r0), sector=(th0 - pi/2, pi/2 - th0),
thickness=2, color='blue', zorder=100)
rmin, rmax = RR(rph_ss(a0)), -1e-8
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2)
g += parametric_plot([R(r)*sin(theta1(a0, lsph(a0, r), qsph(a0, r))),
R(r)*cos(theta1(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2)
g += parametric_plot([R(r)*sin(theta0(a0, lsph(a0, r), qsph(a0, r))),
-R(r)*cos(theta0(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2)
g += parametric_plot([R(r)*sin(theta1(a0, lsph(a0, r), qsph(a0, r))),
-R(r)*cos(theta1(a0, lsph(a0, r), qsph(a0, r)))],
(r, rmin, rmax), color='darkgreen', thickness=2)
ns = 5
dr = (rmax - rmin)/(ns - 1)
for i in range(ns):
r0 = rmin + i*dr
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
th1 = theta1(a0, l0, q0)
g += arc((0, 0), R(r0), sector=(pi/2 - th1, pi/2 - th0), color='green')
g += arc((0, 0), R(r0), sector=(th0 - pi/2, th1 - pi/2), color='green')
g += arc((0,0), R(rph_ss(a0)), sector=(- pi/2, pi/2), color='grey',
linestyle=':', thickness=1)
# inner polar orbits:
r0 = rph_pol_in(a0)
print(n(lsph(a0, r0)))
q0 = qsph(a0, r0)
th0 = 0
th1 = theta1(a0, 0, q0)
print(n(theta0(a0, 0, q0)))
g += arc((0, 0), R(r0), sector=(pi/2 - th1, pi/2 - th0),
linestyle='--', color='green')
g += arc((0, 0), R(r0), sector=(th0 - pi/2, th1 - pi/2),
linestyle='--', color='green')
g
-6.26333640524599e-16 0.000000000000000 6.68118973136174e-16 0.000000000000000
Adding the ergoregion:
x, y = var('x y')
Rxy = sqrt(x^2 + y^2)
r1 = Rxy - 1/Rxy
costh2 = y^2/(x^2 + y^2)
g += region_plot(r1^2 - 2*r1 + a0^2*costh2 < 0, (0, 3), (-3, 3),
incol='whitesmoke', bordercol='grey')
Coloring the region of vortical spherical orbits:
def vortical_orb(x, y):
R0 = sqrt(x*x + y*y)
r0 = R0 - 1/R0
if r0 < rph_ss(a0) or r0 > 0:
return False
phi = atan2(y, x)
th = pi/2 - phi
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
th1 = theta1(a0, l0, q0)
if (th >= th0 and th <= th1) or (th >= pi - th1 and th <= pi - th0):
return True
return False
g += region_plot(vortical_orb, (0, 1), (-1, 1), incol='palegreen', alpha=0.2)
Coloring the region of stable spherical orbits:
def stable_orb(x, y):
R0 = sqrt(x*x + y*y)
r0 = R0 - 1/R0
if r0 <= 0 or r0 > rph_ms(a0):
return False
phi = atan2(y, x)
th = pi/2 - phi
l0 = lsph(a0, r0)
q0 = qsph(a0, r0)
th0 = theta0(a0, l0, q0)
if th < th0 or th > pi - th0:
return False
return True
g += region_plot(stable_orb, (1, 1.5), (-0.7, 0.7), incol='green', alpha=0.5)
Adding circular orbits:
for r1 in [rph_s(a0), rph_p(a0), rph_m(a0)]:
g += point((R(r1), 0), size=60, color='green', zorder=100)
r1 = rph_ss(a0)
sin_th_ss = sqrt( r1/(r1 - 1)*(1 - r1^2/a0^2) )
cos_th_ss = sqrt( 1 - sin_th_ss^2 )
g += point((R(r1)*sin_th_ss, R(r1)*cos_th_ss), size=60, color='green', zorder=100)
g += point((R(r1)*sin_th_ss, -R(r1)*cos_th_ss), size=60, color='green', zorder=100)
show(g, xmin=0, xmax=4.5, ticks=[1, 1], figsize=8)
g.save("gik_spo_meridional.pdf", xmin=0, xmax=4.5, ticks=[1, 1], figsize=8)
show(g, xmin=0, xmax=2.2, ymin=-1, ymax=1, show_legend=False)
g.save("gik_spo_meridional_zoom.pdf", xmin=0, xmax=2.2, ymin=-1, ymax=1,
show_legend=False)