%display latex
Using the manifold catalog (https://trac.sagemath.org/ticket/25869):
K.<t, r, th, ph> = manifolds.Kerr(2, -1, coordinates="BL")
Using the metric to normalize the initial 4-velocity.
g = K.metric()
g[:]
$\tau$ is the proper time.
tau = var('tau')
We choose a starting point $p$ and a normalized initial 4-speed $v$ in the tangent space $T_p$.
p = K((0, 14.98, pi/2, 0))
Tp = K.tangent_space(p)
v = Tp((2, 0, 0.005, 0.05))
v = v / sqrt(-g.at(p)(v, v))
Declaration of the geodesic:
c = K.integrated_geodesic(g, (tau, 0, 3000), v)
Fast integration:
sol = c.solve(step = 1, method="ode_int")
Interpolation:
interp = c.interpolate()
An embedding in $R^3$ is needed to plot the result:
E.<x,y,z> = manifolds.Euclidean(3)
phi = K.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
Plotting the result:
P = c.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)
%time sol = c.solve(step=1, method="ode_int")
CPU times: user 248 ms, sys: 8.38 ms, total: 257 ms Wall time: 244 ms
%time sol = c.solve(step=1, method="rk4_maxima")
CPU times: user 2min 28s, sys: 597 ms, total: 2min 29s Wall time: 2min 47s
%timeit K.default_chart().valid_coordinates(1, 1, 1, 1)
10000 loops, best of 3: 87.5 µs per loop
%timeit K.default_chart().valid_coordinates_numerical(1, 1, 1, 1)
The slowest run took 13.03 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 933 ns per loop