Interpreting Mathematica code on black holes

Question

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.

Answer 1

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]]}]]}