Fast Kerr geodesics

In [1]:
%display latex

Using the manifold catalog (https://trac.sagemath.org/ticket/25869):

In [2]:
K.<t, r, th, ph> = manifolds.Kerr(2, -1, coordinates="BL")

Using the metric to normalize the initial 4-velocity.

In [3]:
g = K.metric()
g[:]
Out[3]:

$\tau$ is the proper time.

In [4]:
tau = var('tau')

We choose a starting point $p$ and a normalized initial 4-speed $v$ in the tangent space $T_p$.

In [5]:
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:

In [6]:
c = K.integrated_geodesic(g, (tau, 0, 3000), v)

Fast integration:

In [7]:
sol = c.solve(step = 1, method="ode_int")

Interpolation:

In [8]:
interp = c.interpolate()

An embedding in $R^3$ is needed to plot the result:

In [9]:
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:

In [10]:
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)

A few speed tests:

In [11]:
%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
In [12]:
%time sol = c.solve(step=1, method="rk4_maxima")
CPU times: user 2min 28s, sys: 597 ms, total: 2min 29s
Wall time: 2min 47s
In [13]:
%timeit K.default_chart().valid_coordinates(1, 1, 1, 1)
10000 loops, best of 3: 87.5 ┬Ás per loop
In [14]:
%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
In [ ]: