Geodesics from initial conditions

[VojtechW]

If I wrote everything correctly, this should add a sub-package InitialConditions to KerrGeodesics that has one key function:

KerrGeoInitOrbit[a,{t0,r0,th0,ph0},u] returns a KerrGeoOrbitFunction[..] (which stores the orbital trajectory and parameters) corresponding to the initial positions and four-velocity u (which does not need to be normalized). One can also optionally provide {[ScriptCapitalE],[ScriptCapitalL], [ScriptCapitalQ]} and {p,e,x} in other arguments.

Even though I made all the commits, about half of the code is actually from the discussion of the issue with Maarten.

[nielsw2]

I wonder if the API should be slightly different. Rather than adding a new public function, all functions could be able to take a,p,e,x and also an association with a different set of initial parameters

So for example you could do:

KerrGeoOrbit[a,p,e,x]

or

KerrGeoOrbit[a,<|"t0"->0, "r0"->10, "th0"->0.3, "ph"->0, "u"-> |>]

This would allow future alternative parameterizations.

[MvdMeent]

One thing that needs to happen is that the Print statements need to be replaced by proper error messages using Message.

[VojtechW]

KerrGeoOrbit[a,<|"t0"->0, "r0"->10, "th0"->0.3, "ph"->0, "u"-> |>]

This seems to be possible in principle, see this post. However, the issue is that I don't know how to do this without breaking the heritage API (i.e. KerrGeoOrbit[a,p,e,x]).

One thing that needs to happen is that the Print statements need to be replaced by proper error messages using Message.

I am not familiar with this too much, but I took a jab at it in the InitialConditions.m. There are many other places in KerrGeodesics where Print[] is still used. (Perhaps a matter for an issue for the whole package...)

[Philip-Lynch]

Sorry for taking so long to get around to this package. It's a nice piece of work and will be a very nice addition once we figure out the best API for it.

I have noticed a bug however when I try to evaluate the trajectory. It seems like the initial polar phase is not passed to the KerrGeoOrbit function correctly, and so I cannot evaluate the functions for t, theta and phi for the orbit's trajectory (r works just fine).

`(Let's give some example orbit) a = 0.9; p = 10; e = 0.4; x = 0.8; [Lambda] = 0;

test = KerrGeoOrbit[a, p, e, x];

t0 = test["Trajectory"][[1]][[Lambda]]; r0 = test["Trajectory"][[2]][[Lambda]]; [Theta]0 = test["Trajectory"][[3]][[Lambda]]; [Phi]0 = test["Trajectory"][[4]][[Lambda]];

u0 = Table[ Values[ KerrGeoFourVelocity[a, p, e, x]][[i]][[Lambda]], {i, 1, 4}];

orbit = KerrGeoInitOrbit[a, {t0, r0, [Theta]0, [Phi]0}, u0]; {t, r, [Theta], [Phi]} = orbit["Trajectory"]; {t[[Lambda]], [Theta][[Lambda]], [Phi][[Lambda]]}`

This should just return the initial positions again, but it instead gives a function in terms of initial polar phase. Same thing happens if evaluated at any other value of Mino time.

[nielsw2]

Just so we don't forget. Following a meeting on the 24th of Feb. we decided on the following API:

KerrGeoOrbit[a,"\[ScriptCapitalE]\[ScriptCapitalL]\[ScriptCapitalQ]"->{0.9,3,1.2},"InitialPhases"->{1,2,4}]
KerrGeoOrbit[a,"\[ScriptCapitalE]\[ScriptCapitalL]\[ScriptCapitalQ]"->{0.9,3,1.2},"OrbitType"->"Plunge"]
KerrGeoOrbit[a,"Subscript[r, ISSO]"->5,"OrbitType"->"Plunge"]
KerrGeoOrbit[a,"InitialPosition"->{1,2,3,4},"FourVelocity"->{1,5,6,7}]

This is the main change needed for before merging this pull request and #44

[VojtechW]

Ok, so now one can access the initial conditions functionality through the KerrGeoOrbit function. Load this branch of the package and this should work (note that also KerrGeoOrbit.m has been modified):

(*Regular bound orbit, you can give the constant of motion instead of orbital elements:*)
KerrGeoOrbit[0.9,"\[ScriptCapitalE]\[ScriptCapitalL]\[ScriptCapitalQ]"->{0.96,3.1`,2.3`},"InitialPhases"->{2,10,15,8}]
(*Specify initial phases:*)
KerrGeoOrbit[0.9,"\[ScriptCapitalE]\[ScriptCapitalL]\[ScriptCapitalQ]"->{0.96,3.1`,2.3`},"InitialPhases"->{2,10,15,8}]
(*Check p,e,x through internal function:*)
elemAssoc= KerrGeodesics`InitialConditions`Private`KerrGeoInit2pex[0.9,{0,0,0,0},{0,0,0,0},{0.96,3.1`,2.3`}]
(*The syntax also allows an alternative to the usual argument pattern*)
{p,e,x} = Values[elemAssoc];
KerrGeoOrbit[0.9,"InitialPhases"->{2,10,15,8},"pex"->{p,e,x}]

You can also specify only coordinate initial conditions:

(*Now let us verify the initial conditions property*)
orb = KerrGeoOrbit[0.9,p,e,x]["Trajectory"];
FV =KerrGeoFourVelocity[0.9,p,e,x];
inVel = {"u^t"[0],"u^r"[0],"u^\[Theta]"[0],"u^\[Phi]"[0]}/.FV
 inPos = {orb[[1]][0],orb[[2]][0],orb[[3]][0],orb[[4]][0]}
(*Check consistency*)
KerrGeodesics`InitialConditions`Private`KerrGeoInit2Constants[0.9,inPos,inVel]
KerrGeoOrbit[0.9,"FourVelocity"->inVel,"InitialPosition"->inPos]

The API currently does not support unbound and plunging/vortical orbits, but stops the computation and provides an error message (previously there was an error message, but the computation kept running nonsensically):

(*Unbound orbit:*)
KerrGeoOrbit[0.9,"\[ScriptCapitalE]\[ScriptCapitalL]\[ScriptCapitalQ]"->{1.1,3.1`,2.3`}]
(*Plunging orbit:*)
KerrGeoOrbit[0.9,"\[ScriptCapitalE]\[ScriptCapitalL]\[ScriptCapitalQ]"->{0.95,3.1`,2.3`}]
(*Vortical orbit:*)
KerrGeoOrbit[0.9,"\[ScriptCapitalE]\[ScriptCapitalL]\[ScriptCapitalQ]"->{0.96,3.1`,-2.3`}]

[VojtechW]

It is not obvious to me why I am failing unit tests, apparently no "TestReport.xml" file is produced. I saw no particular issue when running the package itself. It is true, however, that I did overload the definition of some of the core functions of KerrGeodesics, which may cause issues in some of the unit tests. Any ideas?

[barrywardell]

It is not obvious to me why I am failing unit tests, apparently no "TestReport.xml" file is produced. I saw no particular issue when running the package itself. It is true, however, that I did overload the definition of some of the core functions of KerrGeodesics, which may cause issues in some of the unit tests. Any ideas?

The way they are setup the unit tests only work on the main repository, not on forks. This is because they need access to a Mathematica license.

You can always run the unit tests on your own machine. On the command line, change the the Tests/ directory and run the AllTests.wls script. I have done this with @VojtechW's latest commit and all passed.

[Philip-Lynch]

Only had time to get around to this now. Really nice work on this package. There is just one small thing outstanding. I noticed that orbits with non-zero values of initial phases don't have these taken into account when calculated with Inital Position and FourVelocity. Here's a quick example code to show what I mean:

(*Picking orbital elements*)
{atest,ptest,etest,xtest} = {0.9,10,0.3,0.8};
(*Picking initial phases*)
{t0,qr0,qz0,phi0} = {1,2,3,4};

(*Reference Orbit*)
orbit1 = KerrGeoOrbit[atest,"pex"->{ptest,etest,xtest},"InitialPhases"->{t0,qr0,qz0,phi0}];
{t1,r1,theta1,phi1} = orbit1["Trajectory"];

(*Calculate initial position and velocity from this orbit*)
inPos={t1[0],r1[0],theta1[0],phi1[0]};
FV=KerrGeoFourVelocity[atest,ptest,etest,xtest,{qr0,qz0}];
inVel={Values[FV][[1]][0],Values[FV][[2]][0],Values[FV][[3]][0],Values[FV][[4]][0]};

(*Now calculate the orbit using new code*)
orbit2 = KerrGeoOrbit[atest,"FourVelocity"->inVel,"InitialPosition"->inPos];
{t2,r2,theta2,phi2} = orbit3["Trajectory"];

(*Note that the orbit trajectories do not start at the initial position*)
{{t2[0],r2[0],theta2[0],phi2[0]},inPos}

If the initial phases are zero, they agree, but otherwise they don't. Hopefully this is just a quick fix.