From initial conditions to geodesic

VojtechW commented 2 years ago

Almost anyone who comes to the KerrGeodesics package without knowing that much about Kerr geodesics expects it to do the following:

You pick initial conditions in terms of positions r,theta,t,phi and time derivatives of the coordinates. (This may be wrt Mino time but more likely coordinate time or proper time.)
You get your geodesic starting from that point with that velocity.

Currently, this is actually not implemented in KerrGeodesics, since KerrGeoOrbit gives you a geodesic based on p,e,x and initial orbital phases.

The easiest way to implement the initial-condition approach would be to

Compute E,L,Q from the initial data. (Check normalization/self-consistency of initial data?)
Compute p,e,x from E,L,Q. (Check boundedness/unboundedness of geodesic first?)
Compute phases from coordinate positions and p,e,x. Use this to get corresponding KerrGeoOrbit.

Point 1. and point 3. are easy to do in closed form.

However, point 2. cannot be done fully analytically. The x parameter is a quadratic root. The p,e parameters are then obtained by finding the two "correct" roots of a quartic. For this, one could trust the Mathematica function Root[]. The issue is with indexing the roots, since

The root indexing representation Root[f,k] applies to polynomial functions f only. The indexing of roots takes the real roots first, in increasing order. For polynomials with rational coefficients, the complex conjugate pairs of roots have consecutive indices.

Given numerical initial data, this could make the algorithm make discrete jumps due to indexing changes near transitions between bound/unbound/plunging. But I guess such an implementation is a good start.

Any other thoughts?

MvdMeent commented 2 years ago

Number 2. can be done in closed form analytically at least for bound/unbound orbits. For plunging orbits the p,e,x parametrization breaks down anyway.

VojtechW commented 2 years ago

Number 2. can be done in closed form analytically at least for bound/unbound orbits.

How exactly? I know only of one method, and that is finding the roots of the R(r;E,L,Q) polynomial - and that one is a quartic. (The Schmidt method gives only p,e,x -> E,L,Q, not the other way around afaik.)

For plunging orbits the p,e,x parametrization breaks down anyway.

Sure. I guess this is a nuance, but given numerical data the roots may jump around in an unexpectedly large band around the separatrix. Having discrete jumps due to root indexing is way worse than just getting a noisy inaccurate result.

MvdMeent commented 2 years ago

The four roots of a quartic can be given in closed form. Using the value of E, Q, and the discriminant of the radial potential, you can determine which of the four roots to use (or whether the orbit is plunging). With a bit of effort, you can get a closed form expression for the two roots that define p and e (for non-plunging orbits). Alternatively, given numerical input, you can evaluate all 4 and pick the two that bracket the initial position.

MvdMeent commented 2 years ago

I've written a prototype function, which takes initial data as input and produces (a,p,e,x).

VojtechW commented 2 years ago

Ok, I see. Here is my jab at it as well. First step:

InitKerrGeoConstants[a_, t0_, \[Phi]0_, r0_, \[Theta]0_, dt0_, d\[Phi]0_, dr0_, d\[Theta]0_, Param_ : "Proper"] := 
 Module[{\[CapitalSigma], \[CapitalDelta], gtt, gt\[Phi], g\[Phi]\[Phi], g\[Theta]\[Theta], grr, \[ScriptCapitalN], utm2, En,L, Q},
  (*Kerr metric:*)
  \[CapitalSigma] = r0^2 + a^2 Cos[\[Theta]0]^2;
  \[CapitalDelta] = r0^2 - 2 r0 + a^2;
  gtt = -(1 - (2 r0)/\[CapitalSigma]);
  gt\[Phi] = -((2 r0 a Sin[\[Theta]0]^2)/\[CapitalSigma]);
  g\[Phi]\[Phi] = 
   Sin[\[Theta]0]^2 (r0^2 + a^2 + (
      2 r0 a^2 Sin[\[Theta]0]^2)/\[CapitalSigma]);
  g\[Theta]\[Theta] = \[CapitalSigma];
  grr = \[CapitalSigma]/\[CapitalDelta];
  (*Normalization factors for other parametrizations:*)
  \[ScriptCapitalN] = 1;
  If[Param == "Mino", 
   \[ScriptCapitalN] = 1/\[CapitalSigma], 
   If[Param == "Coord",
    utm2 = -(gtt + 2 gt\[Phi] d\[Phi]0  + g\[Phi]\[Phi] d\[Phi]0^2 + 
        grr dr0^2 + g\[Theta]\[Theta] d\[Theta]0^2);
    If[utm2 > 0, \[ScriptCapitalN] = 1/Sqrt[utm2], 
     "Geodesic not timelike", \[ScriptCapitalN] = 1/Sqrt[utm2]]
    , Null]
   ];
  (*Result:*)
  En = -\[ScriptCapitalN] (gtt dt0 + gt\[Phi] d\[Phi]0);
  L = \[ScriptCapitalN] (  gt\[Phi] dt0 + g\[Phi]\[Phi] d\[Phi]0);
  Q = \[ScriptCapitalN]^2 ((g\[Theta]\[Theta] d\[Theta]0)^2 + (L/Sin[\[Theta]0] - a En Sin[\[Theta]0])^2 + a^2 Cos[\[Theta]0]^2 - (L - a En)^2);
   <|"\[ScriptCapitalE]" -> En, "\[ScriptCapitalL]" -> L, "\[ScriptCapitalQ]" -> Q |>
  ]

Second step using Root[]:

InitKerrGeoPEX[
   a_, \[ScriptCapitalE]_, \[ScriptCapitalL]_, \[ScriptCapitalQ]_] = 
  Module[{R, r1, r2, zm},
   R = (((#^2 + a^2) \[ScriptCapitalE] - a \[ScriptCapitalL])^2 - (#^2 - 2 # + a^2) (#^2 + (\[ScriptCapitalL] - a \[ScriptCapitalE]^2) + \[ScriptCapitalQ])) &;
   r1 = Root[R, 4];
   r2 = Root[R, 3];
   If[PossibleZeroQ[a],
    zm2 = \[ScriptCapitalQ]/(\[ScriptCapitalL]^2 + \[ScriptCapitalQ]),  
    zm2 = (1/(2 a^2 (1 - \[ScriptCapitalE]^2)))(\[ScriptCapitalQ] + \[ScriptCapitalL]^2 + a^2 (1 - \[ScriptCapitalE]^2) - Sqrt[(\[ScriptCapitalQ] + \[ScriptCapitalL]^2 + a^2 (1 - \[ScriptCapitalE]^2))^2 - 4 a^2 (1 - \[ScriptCapitalE]^2) \[ScriptCapitalQ]])
    ];
    <|"p" -> (2 r1 r2)/(r1 + r2), "e" -> (r1 - r2)/(r1 + r2), 
    "x" -> Sign[\[ScriptCapitalL]] Sqrt[1 - zm2] |>
   ];

I assume that Mathematica recognizes the quartic and uses an analytical formula when using Root[].

duetosymmetry commented 2 years ago

Leaving it as Root[] is best, since the explicit form of the quartic roots are very large expressions. Mma knows the analytic properties of Root[] objects (when evaluating them, it may use root polishing instead of the explicit formulas). I'm aware that the ordering of the different Root[] solutions has changed as the version of Mma has changed; I think it's only guaranteed to be stable if they are all real roots. Ideally we would have some combination of a Root[] within an IsolatingIntervalCopy[], but I don't think that exists in Mma.

VojtechW commented 2 years ago

An alternative would be the nighborhood representation of Root[]:

The root neighborhood representation Root[{f,c}] specifies the equation f[x]==0 as well as the neighborhood rectangle centered at c and width 2 10^-Accuracy[Re(c)] and height 2 10^-Accuracy[Im(c)].

Of course, for that one would first need to obtain sufficiently accurate/stable guesses with FindRoot[] or similar.

MvdMeent commented 2 years ago

My version (works for bound and scattering orbits):

Step 1:

Init2ELQ[{M_, a_}, {t0_, r0_, \[Theta]0_, \[Phi]0_}, u_List] := Module[
  {
   g = g[M, a][t0, r0, \[Theta]0, \[Phi]0],
   \[ScriptCapitalK] = \[ScriptCapitalK][M, a][t0, 
     r0, \[Theta]0, \[Phi]0],
   norm},
  norm = u.g.u;
   If[norm < 0,
   <|
    "E" -> {-1, 0, 0, 0}.g.(u/Sqrt[-norm]),
    "L" -> {0, 0, 0, 1}.g.(u/Sqrt[-norm]),
    "Q" -> (u/Sqrt[-norm]).\[ScriptCapitalK].(u/Sqrt[-norm])
    |>,
   Print["Intial velocity is not timelike."];
   $Failed
   ]
  ]

Step 2:

Init2apex[{M_, a_}, {t0_, r0_, \[Theta]0_, \[Phi]0_}, u_List] := 
 Module[
  {
   E, L, Q,
   rts,
   disc,
   r1, r2
   },
  {E, L, Q} = (Values@
     Init2ELQ[{M, a}, {t0, r0, \[Theta]0, \[Phi]0}, u]);
  Print[{E, L, Q}];
  disc = Disc[a, E, L, Q];
  rts = Sort@Re@roots[a, E, L, Q][[All, 1, 2]];
  Which[
   Q < 0, Print["Vortical orbit detected"],
   disc < 0, Print["Plunge orbit detected"],
   True,
   r1 = If[E < 1, rts[[4]], rts[[1]]];
   r2 = If[E < 1, rts[[3]], rts[[4]]];
   ];
  <|"a" -> a, "p" -> (2 r1 r2)/(r1 + r2), "e" -> (r1 - r2)/(r1 + r2), 
   "x" -> Sign[L] Sqrt[1 - ELQ2zmax[a, E, L, Q]^2] |>
  ]

With a bunch of explicit definitions:

R = Function[
   r, (\[ScriptCapitalE] (r^2 + a^2) - 
      a \[ScriptCapitalL])^2 - (r^2 - 2 r + 
       a^2) (-\[Mu] r^2 + (\[ScriptCapitalL] - 
         a \[ScriptCapitalE])^2 + \[ScriptCapitalQ])];
Z = Function[
   z, \[ScriptCapitalQ] - 
    z^2 (a^2 (1 - \[ScriptCapitalE]^2) (1 - 
          z^2) + \[ScriptCapitalL]^2 + \[ScriptCapitalQ])];

ELQ2zmax[a_, \[ScriptCapitalE]_, \[ScriptCapitalL]_, \
\[ScriptCapitalQ]_] := (
 a^2 - a^2 \[ScriptCapitalE]^2 + \[ScriptCapitalL]^2 + \
\[ScriptCapitalQ] - 
  Sqrt[(-a^2 + 
     a^2 \[ScriptCapitalE]^2 - \[ScriptCapitalL]^2 - \
\[ScriptCapitalQ])^2 - 
   4 (a^2 - a^2 \[ScriptCapitalE]^2) \[ScriptCapitalQ]])/(
 2 (a^2 - a^2 \[ScriptCapitalE]^2))

Disc[a_, \[ScriptCapitalE]_, \[ScriptCapitalL]_, \[ScriptCapitalQ]_] \
:= -16 (-a^8 \[ScriptCapitalE]^4 + a^10 \[ScriptCapitalE]^4 + 
    16 a^6 \[ScriptCapitalE]^6 - 16 a^8 \[ScriptCapitalE]^6 - 
    4 a^10 \[ScriptCapitalE]^6 + 62 a^8 \[ScriptCapitalE]^8 + 
    6 a^10 \[ScriptCapitalE]^8 - 72 a^8 \[ScriptCapitalE]^10 - 
    4 a^10 \[ScriptCapitalE]^10 + 27 a^8 \[ScriptCapitalE]^12 + 
    a^10 \[ScriptCapitalE]^12 + 
    4 a^7 \[ScriptCapitalE]^3 \[ScriptCapitalL] - 
    4 a^9 \[ScriptCapitalE]^3 \[ScriptCapitalL] - 
    96 a^5 \[ScriptCapitalE]^5 \[ScriptCapitalL] + 
    100 a^7 \[ScriptCapitalE]^5 \[ScriptCapitalL] + 
    16 a^9 \[ScriptCapitalE]^5 \[ScriptCapitalL] - 
    428 a^7 \[ScriptCapitalE]^7 \[ScriptCapitalL] - 
    24 a^9 \[ScriptCapitalE]^7 \[ScriptCapitalL] + 
    540 a^7 \[ScriptCapitalE]^9 \[ScriptCapitalL] + 
    16 a^9 \[ScriptCapitalE]^9 \[ScriptCapitalL] - 
    216 a^7 \[ScriptCapitalE]^11 \[ScriptCapitalL] - 
    4 a^9 \[ScriptCapitalE]^11 \[ScriptCapitalL] - 
    6 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 + 
    6 a^8 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 + 
    240 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 - 
    260 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 - 
    21 a^8 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 + 
    1274 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 + 
    27 a^8 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 - 
    1764 a^6 \[ScriptCapitalE]^8 \[ScriptCapitalL]^2 - 
    15 a^8 \[ScriptCapitalE]^8 \[ScriptCapitalL]^2 + 
    756 a^6 \[ScriptCapitalE]^10 \[ScriptCapitalL]^2 + 
    3 a^8 \[ScriptCapitalE]^10 \[ScriptCapitalL]^2 + 
    4 a^5 \[ScriptCapitalE] \[ScriptCapitalL]^3 - 
    4 a^7 \[ScriptCapitalE] \[ScriptCapitalL]^3 - 
    320 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 + 
    360 a^5 \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 + 
    4 a^7 \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 - 
    2128 a^5 \[ScriptCapitalE]^5 \[ScriptCapitalL]^3 + 
    12 a^7 \[ScriptCapitalE]^5 \[ScriptCapitalL]^3 + 
    3276 a^5 \[ScriptCapitalE]^7 \[ScriptCapitalL]^3 - 
    20 a^7 \[ScriptCapitalE]^7 \[ScriptCapitalL]^3 - 
    1512 a^5 \[ScriptCapitalE]^9 \[ScriptCapitalL]^3 + 
    8 a^7 \[ScriptCapitalE]^9 \[ScriptCapitalL]^3 - 
    a^4 \[ScriptCapitalL]^4 + a^6 \[ScriptCapitalL]^4 + 
    240 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 - 
    280 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 + 
    14 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 + 
    2170 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 - 
    45 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 - 
    3780 a^4 \[ScriptCapitalE]^6 \[ScriptCapitalL]^4 + 
    44 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalL]^4 + 
    1890 a^4 \[ScriptCapitalE]^8 \[ScriptCapitalL]^4 - 
    14 a^6 \[ScriptCapitalE]^8 \[ScriptCapitalL]^4 - 
    96 a \[ScriptCapitalE] \[ScriptCapitalL]^5 + 
    116 a^3 \[ScriptCapitalE] \[ScriptCapitalL]^5 - 
    12 a^5 \[ScriptCapitalE] \[ScriptCapitalL]^5 - 
    1372 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL]^5 + 
    24 a^5 \[ScriptCapitalE]^3 \[ScriptCapitalL]^5 + 
    2772 a^3 \[ScriptCapitalE]^5 \[ScriptCapitalL]^5 - 
    12 a^5 \[ScriptCapitalE]^5 \[ScriptCapitalL]^5 - 
    1512 a^3 \[ScriptCapitalE]^7 \[ScriptCapitalL]^5 + 
    16 \[ScriptCapitalL]^6 - 20 a^2 \[ScriptCapitalL]^6 + 
    3 a^4 \[ScriptCapitalL]^6 + 
    518 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^6 + 
    9 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^6 - 
    1260 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalL]^6 - 
    26 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^6 + 
    756 a^2 \[ScriptCapitalE]^6 \[ScriptCapitalL]^6 + 
    14 a^4 \[ScriptCapitalE]^6 \[ScriptCapitalL]^6 - 
    104 a \[ScriptCapitalE] \[ScriptCapitalL]^7 - 
    12 a^3 \[ScriptCapitalE] \[ScriptCapitalL]^7 + 
    324 a \[ScriptCapitalE]^3 \[ScriptCapitalL]^7 + 
    20 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL]^7 - 
    216 a \[ScriptCapitalE]^5 \[ScriptCapitalL]^7 - 
    8 a^3 \[ScriptCapitalE]^5 \[ScriptCapitalL]^7 + 
    8 \[ScriptCapitalL]^8 + 3 a^2 \[ScriptCapitalL]^8 - 
    36 \[ScriptCapitalE]^2 \[ScriptCapitalL]^8 + 
    27 \[ScriptCapitalE]^4 \[ScriptCapitalL]^8 - 
    3 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalL]^8 - 
    4 a \[ScriptCapitalE] \[ScriptCapitalL]^9 + 
    4 a \[ScriptCapitalE]^3 \[ScriptCapitalL]^9 + \
\[ScriptCapitalL]^10 - \[ScriptCapitalE]^2 \[ScriptCapitalL]^10 + 
    a^8 \[ScriptCapitalQ] - a^10 \[ScriptCapitalQ] - 
    20 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalQ] + 
    19 a^8 \[ScriptCapitalE]^2 \[ScriptCapitalQ] + 
    5 a^10 \[ScriptCapitalE]^2 \[ScriptCapitalQ] + 
    48 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalQ] - 
    28 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalQ] - 
    98 a^8 \[ScriptCapitalE]^4 \[ScriptCapitalQ] - 
    10 a^10 \[ScriptCapitalE]^4 \[ScriptCapitalQ] + 
    192 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalQ] + 
    170 a^8 \[ScriptCapitalE]^6 \[ScriptCapitalQ] + 
    10 a^10 \[ScriptCapitalE]^6 \[ScriptCapitalQ] - 
    252 a^6 \[ScriptCapitalE]^8 \[ScriptCapitalQ] - 
    127 a^8 \[ScriptCapitalE]^8 \[ScriptCapitalQ] - 
    5 a^10 \[ScriptCapitalE]^8 \[ScriptCapitalQ] + 
    108 a^6 \[ScriptCapitalE]^10 \[ScriptCapitalQ] + 
    35 a^8 \[ScriptCapitalE]^10 \[ScriptCapitalQ] + 
    a^10 \[ScriptCapitalE]^10 \[ScriptCapitalQ] + 
    40 a^5 \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ] - 
    44 a^7 \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ] - 
    192 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ] + 
    156 a^5 \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ] + 
    268 a^7 \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ] - 
    952 a^5 \[ScriptCapitalE]^5 \[ScriptCapitalL] \[ScriptCapitalQ] - 
    540 a^7 \[ScriptCapitalE]^5 \[ScriptCapitalL] \[ScriptCapitalQ] + 
    1404 a^5 \[ScriptCapitalE]^7 \[ScriptCapitalL] \[ScriptCapitalQ] \
+ 452 a^7 \[ScriptCapitalE]^7 \[ScriptCapitalL] \[ScriptCapitalQ] - 
    648 a^5 \[ScriptCapitalE]^9 \[ScriptCapitalL] \[ScriptCapitalQ] - 
    136 a^7 \[ScriptCapitalE]^9 \[ScriptCapitalL] \[ScriptCapitalQ] - 
    20 a^4 \[ScriptCapitalL]^2 \[ScriptCapitalQ] + 
    25 a^6 \[ScriptCapitalL]^2 \[ScriptCapitalQ] - 
    4 a^8 \[ScriptCapitalL]^2 \[ScriptCapitalQ] + 
    288 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ] \
- 300 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ] - 
    226 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ] \
+ 16 a^8 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ] + 
    1920 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ] + 
    541 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \[ScriptCapitalQ] \
- 24 a^8 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \[ScriptCapitalQ] - 
    3240 a^4 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ] - 
    504 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 \[ScriptCapitalQ] \
+ 16 a^8 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 \[ScriptCapitalQ] + 
    1620 a^4 \[ScriptCapitalE]^8 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ] + 
    164 a^6 \[ScriptCapitalE]^8 \[ScriptCapitalL]^2 \[ScriptCapitalQ] \
- 4 a^8 \[ScriptCapitalE]^8 \[ScriptCapitalL]^2 \[ScriptCapitalQ] - 
    192 a \[ScriptCapitalE] \[ScriptCapitalL]^3 \[ScriptCapitalQ] + 
    244 a^3 \[ScriptCapitalE] \[ScriptCapitalL]^3 \[ScriptCapitalQ] + 
    40 a^5 \[ScriptCapitalE] \[ScriptCapitalL]^3 \[ScriptCapitalQ] - 
    2000 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 \
\[ScriptCapitalQ] - 
    80 a^5 \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 \[ScriptCapitalQ] \
+ 3960 a^3 \[ScriptCapitalE]^5 \[ScriptCapitalL]^3 \[ScriptCapitalQ] \
+ 40 a^5 \[ScriptCapitalE]^5 \[ScriptCapitalL]^3 \[ScriptCapitalQ] - 
    2160 a^3 \[ScriptCapitalE]^7 \[ScriptCapitalL]^3 \
\[ScriptCapitalQ] + 48 \[ScriptCapitalL]^4 \[ScriptCapitalQ] - 
    72 a^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ] + 
    16 a^4 \[ScriptCapitalL]^4 \[ScriptCapitalQ] - 
    6 a^6 \[ScriptCapitalL]^4 \[ScriptCapitalQ] + 
    1120 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \
\[ScriptCapitalQ] - 
    156 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ] \
+ 18 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ] - 
    2700 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 \
\[ScriptCapitalQ] + 
    290 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 \[ScriptCapitalQ] \
- 18 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 \[ScriptCapitalQ] + 
    1620 a^2 \[ScriptCapitalE]^6 \[ScriptCapitalL]^4 \
\[ScriptCapitalQ] - 
    150 a^4 \[ScriptCapitalE]^6 \[ScriptCapitalL]^4 \[ScriptCapitalQ] \
+ 6 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalL]^4 \[ScriptCapitalQ] - 
    312 a \[ScriptCapitalE] \[ScriptCapitalL]^5 \[ScriptCapitalQ] + 
    68 a^3 \[ScriptCapitalE] \[ScriptCapitalL]^5 \[ScriptCapitalQ] + 
    972 a \[ScriptCapitalE]^3 \[ScriptCapitalL]^5 \[ScriptCapitalQ] - 
    188 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL]^5 \[ScriptCapitalQ] \
- 648 a \[ScriptCapitalE]^5 \[ScriptCapitalL]^5 \[ScriptCapitalQ] + 
    120 a^3 \[ScriptCapitalE]^5 \[ScriptCapitalL]^5 \[ScriptCapitalQ] \
+ 32 \[ScriptCapitalL]^6 \[ScriptCapitalQ] - 
    3 a^2 \[ScriptCapitalL]^6 \[ScriptCapitalQ] - 
    4 a^4 \[ScriptCapitalL]^6 \[ScriptCapitalQ] - 
    144 \[ScriptCapitalE]^2 \[ScriptCapitalL]^6 \[ScriptCapitalQ] + 
    48 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^6 \[ScriptCapitalQ] \
+ 8 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^6 \[ScriptCapitalQ] + 
    108 \[ScriptCapitalE]^4 \[ScriptCapitalL]^6 \[ScriptCapitalQ] - 
    44 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalL]^6 \[ScriptCapitalQ] \
- 4 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^6 \[ScriptCapitalQ] - 
    16 a \[ScriptCapitalE] \[ScriptCapitalL]^7 \[ScriptCapitalQ] + 
    16 a \[ScriptCapitalE]^3 \[ScriptCapitalL]^7 \[ScriptCapitalQ] + 
    5 \[ScriptCapitalL]^8 \[ScriptCapitalQ] - 
    a^2 \[ScriptCapitalL]^8 \[ScriptCapitalQ] - 
    5 \[ScriptCapitalE]^2 \[ScriptCapitalL]^8 \[ScriptCapitalQ] + 
    a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^8 \[ScriptCapitalQ] + 
    8 a^4 \[ScriptCapitalQ]^2 - 12 a^6 \[ScriptCapitalQ]^2 + 
    4 a^8 \[ScriptCapitalQ]^2 + 
    48 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^2 - 
    44 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^2 + 
    24 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^2 - 
    16 a^8 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^2 + 
    206 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^2 + 
    22 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^2 + 
    24 a^8 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^2 - 
    324 a^4 \[ScriptCapitalE]^6 \[ScriptCapitalQ]^2 - 
    68 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalQ]^2 - 
    16 a^8 \[ScriptCapitalE]^6 \[ScriptCapitalQ]^2 + 
    162 a^4 \[ScriptCapitalE]^8 \[ScriptCapitalQ]^2 + 
    34 a^6 \[ScriptCapitalE]^8 \[ScriptCapitalQ]^2 + 
    4 a^8 \[ScriptCapitalE]^8 \[ScriptCapitalQ]^2 - 
    96 a \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ]^2 + 
    128 a^3 \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ]^2 - 
    44 a^5 \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ]^2 - 
    628 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ]^2 \
+ 88 a^5 \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ]^2 + 
    1188 a^3 \[ScriptCapitalE]^5 \[ScriptCapitalL] \
\[ScriptCapitalQ]^2 - 
    44 a^5 \[ScriptCapitalE]^5 \[ScriptCapitalL] \[ScriptCapitalQ]^2 \
- 648 a^3 \[ScriptCapitalE]^7 \[ScriptCapitalL] \[ScriptCapitalQ]^2 + 
    48 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^2 - 
    84 a^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^2 + 
    35 a^4 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^2 + 
    4 a^6 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^2 + 
    686 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 - 
    255 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 - 
    12 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 - 
    1620 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 + 
    418 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 + 
    12 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 + 
    972 a^2 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 - 
    198 a^4 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^2 - 
    4 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^2 \
- 312 a \[ScriptCapitalE] \[ScriptCapitalL]^3 \[ScriptCapitalQ]^2 + 
    172 a^3 \[ScriptCapitalE] \[ScriptCapitalL]^3 \[ScriptCapitalQ]^2 \
+ 972 a \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 \[ScriptCapitalQ]^2 - 
    436 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 \
\[ScriptCapitalQ]^2 - 
    648 a \[ScriptCapitalE]^5 \[ScriptCapitalL]^3 \[ScriptCapitalQ]^2 \
+ 264 a^3 \[ScriptCapitalE]^5 \[ScriptCapitalL]^3 \[ScriptCapitalQ]^2 \
+ 48 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^2 - 
    27 a^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^2 - 
    4 a^4 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^2 - 
    216 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^2 + 
    144 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \
\[ScriptCapitalQ]^2 + 
    8 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^2 \
+ 162 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^2 - 
    114 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 \
\[ScriptCapitalQ]^2 - 
    4 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^2 \
- 24 a \[ScriptCapitalE] \[ScriptCapitalL]^5 \[ScriptCapitalQ]^2 + 
    24 a \[ScriptCapitalE]^3 \[ScriptCapitalL]^5 \[ScriptCapitalQ]^2 \
+ 10 \[ScriptCapitalL]^6 \[ScriptCapitalQ]^2 - 
    4 a^2 \[ScriptCapitalL]^6 \[ScriptCapitalQ]^2 - 
    10 \[ScriptCapitalE]^2 \[ScriptCapitalL]^6 \[ScriptCapitalQ]^2 + 
    4 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^6 \[ScriptCapitalQ]^2 \
+ 16 \[ScriptCapitalQ]^3 - 32 a^2 \[ScriptCapitalQ]^3 + 
    22 a^4 \[ScriptCapitalQ]^3 - 6 a^6 \[ScriptCapitalQ]^3 + 
    84 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^3 - 
    90 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^3 + 
    18 a^6 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^3 - 
    180 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^3 + 
    102 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^3 - 
    18 a^6 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^3 + 
    108 a^2 \[ScriptCapitalE]^6 \[ScriptCapitalQ]^3 - 
    34 a^4 \[ScriptCapitalE]^6 \[ScriptCapitalQ]^3 + 
    6 a^6 \[ScriptCapitalE]^6 \[ScriptCapitalQ]^3 - 
    104 a \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ]^3 + 
    92 a^3 \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ]^3 + 
    324 a \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ]^3 - 
    228 a^3 \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ]^3 \
- 216 a \[ScriptCapitalE]^5 \[ScriptCapitalL] \[ScriptCapitalQ]^3 + 
    136 a^3 \[ScriptCapitalE]^5 \[ScriptCapitalL] \[ScriptCapitalQ]^3 \
+ 32 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^3 - 
    33 a^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^3 + 
    4 a^4 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^3 - 
    144 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^3 + 
    144 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^3 - 
    8 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^3 \
+ 108 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^3 - 
    108 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \
\[ScriptCapitalQ]^3 + 
    4 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^3 \
- 16 a \[ScriptCapitalE] \[ScriptCapitalL]^3 \[ScriptCapitalQ]^3 + 
    16 a \[ScriptCapitalE]^3 \[ScriptCapitalL]^3 \[ScriptCapitalQ]^3 \
+ 10 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^3 - 
    6 a^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^3 - 
    10 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^3 + 
    6 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^4 \[ScriptCapitalQ]^3 \
+ 8 \[ScriptCapitalQ]^4 - 12 a^2 \[ScriptCapitalQ]^4 + 
    4 a^4 \[ScriptCapitalQ]^4 - 
    36 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^4 + 
    48 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^4 - 
    8 a^4 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^4 + 
    27 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^4 - 
    35 a^2 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^4 + 
    4 a^4 \[ScriptCapitalE]^4 \[ScriptCapitalQ]^4 - 
    4 a \[ScriptCapitalE] \[ScriptCapitalL] \[ScriptCapitalQ]^4 + 
    4 a \[ScriptCapitalE]^3 \[ScriptCapitalL] \[ScriptCapitalQ]^4 + 
    5 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^4 - 
    4 a^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^4 - 
    5 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^4 + 
    4 a^2 \[ScriptCapitalE]^2 \[ScriptCapitalL]^2 \[ScriptCapitalQ]^4 \
+ \[ScriptCapitalQ]^5 - 
    a^2 \[ScriptCapitalQ]^5 - \[ScriptCapitalE]^2 \[ScriptCapitalQ]^5 \
+ a^2 \[ScriptCapitalE]^2 \[ScriptCapitalQ]^5)

roots[a_, \[ScriptCapitalE]_, \[ScriptCapitalL]_, \[ScriptCapitalQ]_] \
:= Evaluate@Solve[R[r] == 0 /. \[Mu] -> -1, r]

g[M_, a_][t_, 
  r_, \[Theta]_, \[Phi]_] := {{-1 + (2 r)/(r^2 + a^2 Cos[\[Theta]]^2),
    0, 0, -((2 a r Sin[\[Theta]]^2)/(
    r^2 + a^2 Cos[\[Theta]]^2))}, {0, (r^2 + a^2 Cos[\[Theta]]^2)/(
   a^2 + (-2 + r) r), 0, 0}, {0, 0, r^2 + a^2 Cos[\[Theta]]^2, 
   0}, {-((2 a r Sin[\[Theta]]^2)/(r^2 + a^2 Cos[\[Theta]]^2)), 0, 
   0, (Sin[\[Theta]]^2 ((a^2 + r^2) (r^2 + a^2 Cos[\[Theta]]^2) + 
      2 a^2 r Sin[\[Theta]]^2))/(r^2 + a^2 Cos[\[Theta]]^2)}}

\[ScriptCapitalK][M_, a_][t_, 
  r_, \[Theta]_, \[Phi]_] := {{(
   a^2 r Cos[\[Theta]]^2 (a^2 + 
      2 (-1 + r) r + (a^2 - 2 r) Cos[2 \[Theta]]))/(r^2 + 
     a^2 Cos[\[Theta]]^2)^2, 0, 
   0, -((a r (a^4 + 2 r^4 + a^2 r (-2 + 3 r) + 
       a^2 (a^2 + (-2 + r) r) Cos[2 \[Theta]]) Sin[2 \[Theta]]^2)/(
    4 (r^2 + a^2 Cos[\[Theta]]^2)^2))}, {0, -((
    a^2 Cos[\[Theta]]^2 (r^2 + a^2 Cos[\[Theta]]^2))/(
    a^2 + (-2 + r) r)), 0, 0}, {0, 0, r^4 + a^2 r^2 Cos[\[Theta]]^2, 
   0}, {-((a r (a^4 + 2 r^4 + a^2 r (-2 + 3 r) + 
       a^2 (a^2 + (-2 + r) r) Cos[2 \[Theta]]) Sin[2 \[Theta]]^2)/(
    4 (r^2 + a^2 Cos[\[Theta]]^2)^2)), 0, 0, (
   r Cos[\[Theta]]^2 (2 a^6 - 4 a^4 r + 3 a^6 r + 12 a^4 r^2 + 
      11 a^4 r^3 + 16 a^2 r^4 + 16 a^2 r^5 + 8 r^7 + 
      4 a^2 r (a^4 + 2 (-2 + r) r^3 + a^2 r (-2 + 3 r)) Cos[
        2 \[Theta]] + 
      a^4 (-2 + r) (a^2 + (-2 + r) r) Cos[
        4 \[Theta]]) Sin[\[Theta]]^2)/(
   8 (r^2 + a^2 Cos[\[Theta]]^2)^2)}}

VojtechW commented 2 years ago

Cool, I created a branch InitialConditions on a fork that uses this method. It also contains the phase and orbit computations, even though I have some issues with the optional arguments rn.

VojtechW commented 2 years ago

Ok, pull requested, for more see https://github.com/VojtechW/KerrGeodesics/tree/InitialConditions.

BlackHolePerturbationToolkit / KerrGeodesics

From initial conditions to geodesic #37