3
$\begingroup$

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.
$\endgroup$
  • 2
    $\begingroup$ Did you contact the author for explanation? Or have you decided not to do that? $\endgroup$ – Somos Jan 21 at 20:21
  • $\begingroup$ 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 $\endgroup$ – user61882 Jan 22 at 13:42
  • $\begingroup$ This is outdated code. I checked that it does not work on version 11.3. We must write a new one. $\endgroup$ – Alex Trounev Jan 22 at 15:08
  • $\begingroup$ 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? $\endgroup$ – user61882 Jan 23 at 22:08
  • $\begingroup$ @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. $\endgroup$ – Alex Trounev Jan 25 at 16:15
15
$\begingroup$

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

fig1

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

fig2

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

fig3

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – user61882 Jan 23 at 19:24
  • $\begingroup$ 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. $\endgroup$ – Alex Trounev Jan 23 at 20:32
  • $\begingroup$ @BobHanlon Thank you! $\endgroup$ – Alex Trounev Jan 23 at 20:56
  • $\begingroup$ Awesome, so basically you can fiddle around with the parameters and get whatever trajectory you would like. For instance by setting uinvar = 0; a = 0 you get a light trajectory in a Schwarzschild geometry too. $\endgroup$ – user61882 Jan 23 at 21:51
  • $\begingroup$ 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)? $\endgroup$ – user61882 Jan 23 at 22:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.