# Interpreting Mathematica code on black holes

I am trying to understand the code written down on page 7 of this document (code is in Mathematica) I understand pretty much all of the code on the previous page needed to setup the page 7 code (setting up geodesics etc):

https://www.physics.uci.edu/~etolleru/KerrOrbitProject.pdf

Would appreciate if someone could explain the 2 large code blocks (preferably annotate or somehow show what each part means with respect to the final results/graphs). I get vaguely parts what they are doing but not all of it.

From what I understand,

• 1st block: Sets up code for initial conditions and solutions for the geodesics equations as a list of functions. Not 100% sure on there meaning.
• 2nd block: Putting in actual numerical values for these functions for specific initial conditions and then plotting the curve.

Useful information about the document:

• The code is designed to graph particle orbits around a Kerr black hole. However as a build up to looking at the Kerr black hole the author decides to look at Minkowski and Schwarzschild metric first as they require similiar setup. (code on pg 7 is looking at the Schwarzschild solution).
• The description of what the code on pg7 is doing is at the bottom of pg6. I am just unsure of how it is actually being done.
• Did you contact the author for explanation? Or have you decided not to do that? Jan 21, 2019 at 20:21
• Hey, no I haven't contacted the original author primarily because the paper is over 10years old and didn't quite know how long it would take etc. Thought that If I posted on here someone would be able to figure it out Jan 22, 2019 at 13:42
• This is outdated code. I checked that it does not work on version 11.3. We must write a new one. Jan 22, 2019 at 15:08
• Thanks Alex, how did you know how to edit the old code? Was it from understanding the error messages and then translating the maths over into new code? Jan 23, 2019 at 22:08
• @user61882 Yes, the system localizes errors, and I correct them. It is not difficult, but you need to know the language well and understand what the result should be. Jan 25, 2019 at 16:15

## 1 Answer

I managed to debug the code from KerrOrbitGRProject and reproduce the results for the Schwarzschild metric (figures on pp. 7-8)

coords = {t, r, θ, φ}; n =
Length[coords]; a = 0; ρ =
r^2 + a^2 Cos[θ]^2; Δ = r^2 \[Minus] 2*M*r + a^2;
tt = 2 M*r/ρ \[Minus]
1; rr = ρ/Δ; θθ = ρ; \
φφ = (Δ + (2 M*
r*(r^2 +
a^2))/ρ) Sin[θ]^2; tφ = \[Minus]4 a M \
r Sin[θ]^2/ρ; metric = {{tt, 0, 0, tφ}, {0, rr,
0, 0}, {0, 0, θθ, 0}, {tφ, 0,
0, φφ}}; metric // MatrixForm
inversemetric = Simplify[Inverse[metric]]; inversemetric // MatrixForm
christoffel :=
christoffel =
Simplify[Table[
1/2*Sum[inversemetric[[i,
s]]*(D[metric[[s, j]], coords[[k]]] +
D[metric[[s, k]], coords[[j]]] \[Minus]
D[metric[[j, k]], coords[[s]]]), {s, 1, n}], { i, 1, n}, {j,
1, n}, {k, 1, n}]]
listchristoffel :=
Table[If[UnsameQ[christoffel[[i, j, k]],
0 ], {ToString[Γ[i, j, k]],
christoffel[[i, j, k]]}], {i, 1, n}, {j, 1, n}, {k, 1,
j}] TableForm[
Partition[DeleteCases[Flatten[listchristoffel], Null], 2 ],
TableSpacing -> {2, 2}]
geodesic :=
geodesic =
Simplify[Table[\[Minus]Sum[
christoffel[[i, j, k]] coords[[j]]' coords[[k]]', {j, 1, n}, {k,
1, n}], {i, 1, n}]]
listgeodesic :=
Table[{"d/dτ" ToString[coords[[i]]'], " =", geodesic[[i]]}, {i,
1, n}]
TableForm[listgeodesic, TableSpacing -> {2}]
maxτ = 750; ivs = {0, 0, .088}; ics = {0, 6.5, π/2, 0}; M = 1;
computeSoln[maxτi_, ivsi_, icsi_] :=
Block[{ivs, ics, i, χ, tmp, soln}, ics = icsi;
ivs = Join[{χ}, ivsi];
op1 = Table[coords[[i]] -> ics[[i]], {i, 0, n}];
tm = metric /. op1;
tmp = ivs.(tm.ivs); χslv = Solve[tmp == uinvar, χ];
ivs[[1]] = Last[χ /. χslv];

op = {Derivative[1][t] -> Derivative[1][t][τ],
Derivative[1][r] -> Derivative[1][r][τ],
Derivative[1][θ] -> Derivative[1][θ][τ],
Derivative[1][φ] -> Derivative[1][φ][τ],
t -> t[τ],
r -> r[τ], θ -> θ[τ], φ -> \
φ[τ]};
deq = Table[
coords[[i]]''[τ] == Simplify[geodesic[[i]] /. op], {i, 1,
n}]; deq =
Join[deq, Table[coords[[i]]'[0] == ivs[[i]], {i, 1, n}],
Table[coords[[i]][0] == ics[[i]], {i, 1, n}]];
soln = NDSolve[deq, coords, {τ, 0, maxτi}]; soln]
uinvar = \[Minus]1;
sphslnToCartsln[soln_] :=
Block[{xs, ys, zs},
xs = r[τ] Sin[θ[τ]] Cos[φ[τ]] /.
soln; ys =
r[τ] Sin[θ[τ]] Sin[φ[τ]] /. soln;
zs = r[τ] Cos[θ[τ]] /. soln; {xs, ys, zs}]
udotu[solni_, τval_] :=
Block[{xα, uα},
xα =
Table[coords[[i]][τ] /. solni, {i, 1, n}] // Flatten;
uα = D[xα, τ];
xα = xα /. τ -> τval;
uα = uα /. τ -> τval;
uα.((metric /.
Table[coords[[i]] -> xα[[i]], {i, 1, n}]).uα)]
coordlist[τin_] :=
Table[ToString[coords[[i]]] <>
" = " <> {ToString[coords[[i]][τin] /. soln // First]}, {i, 1,
n}]

soln = computeSoln[maxτ, ivs, ics];

xyzsoln = sphslnToCartsln[soln];
Join[{"Final Coordinates:"}, coordlist[maxτ]] // TableForm
Join[{{"", "", "", "u.u values"}},
Table[{"τ=", ToString[i], "->", udotu[soln, i]}, {i, 0,
maxτ, maxτ/5}]] // TableForm
{Plot[Evaluate[
Table[coords[[i]][τ] /. soln, {i, 1, n}]], {τ, 0,
maxτ}, AxesLabel -> {"τ", "Coordinate"},
PlotLegends -> {"t", "r", "θ", "φ"},
PlotRange -> {0, 30}],
Show[ParametricPlot[
Evaluate[{xyzsoln[[1]], xyzsoln[[2]]} // Flatten], { τ, 0,
maxτ}, AspectRatio -> 1, PlotStyle -> Gray],
Graphics[{Red, Circle[{0, 0}, 2]}]]}


To get the pictures on pp. 23-24, we need to change the settings

a=.9; maxτ = 200; ivs = {-0.15, 0, .002}; ics = {0, 10, π/2,
Pi/4}; M = 1; uinvar = 0;


and print options

{Plot[Evaluate[
Table[coords[[i]][τ] /. soln, {i, 1, n}]], {τ, 0,
maxτ}, AxesLabel -> {"τ", "Coordinate"},
PlotLegends -> {"t", "r", "θ", "φ"},
PlotRange -> {0, 30}],
Show[ParametricPlot[
Evaluate[{xyzsoln[[1]], xyzsoln[[2]]} // Flatten], { τ, 0,
maxτ}, AspectRatio -> 1, PlotStyle -> Gray],
Graphics[{Red, Circle[{0, 0}, M + Sqrt[M^2 - a^2]]}]]}


Consider the scattering of light rays in the Schwarzschild metric

a=0; max\[Tau] = 150; ivs = {-0.15, 0, .02}; ics = {0, 10, \[Pi]/4,
Pi/4}; M = 1;
Show[ParametricPlot[
Evaluate[{xyzsoln[[1]], xyzsoln[[2]]} // Flatten], { \[Tau], 0,
max\[Tau]}, AspectRatio -> 1, PlotStyle -> Gray,
PlotRange -> {{-10, 5}, {-7, 8}}],
Graphics[{Red, Circle[{0, 0}, M + Sqrt[M^2 - a^2]]}]]}


• This is great, thanks! Just one thing, if I wanted to do a similar thing but get the trajectory of a light ray in a Schwarzschild orbit rather than a massive particle (what the code in this document is doing) would it be correct just to edit the "univar" function and set it equal to 0 rather than -1? It has a statement about this at the bottom of page 3 (how they are considering massive particles) and at the bottom of page 6 about setting up the initial conditions. Jan 23, 2019 at 19:24
• If you want to get pictures on pages 23-24, then you need to put uinvar = 0; a = 0.9 and change other parameters. See update to my answer. Jan 23, 2019 at 20:32
• @BobHanlon Thank you! Jan 23, 2019 at 20:56
• How could I fiddle the other parameters like maxt, ivs and ics (for Schwarzschild light orbit its obvious that Univar = 0 and a=0) to get a meaningful plot for the Schwarzschild light orbit curves (Not the Kerr as is done in the booklet)? Jan 23, 2019 at 22:05
• @Ordinaryusers68 I have tested code on 12.0.0 and 12.1.1 for Windows 10. In the last picture we should use AspectRatio -> Automatic. Jul 17, 2020 at 10:47