Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Issue50 #51

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Issue50 #51

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
6e9d7bb
Simplified Schwarzschild circ velocities so they are finite on ISCO #50
Philip-Lynch Jul 12, 2023
ab16e47
Created KerrGeoOnSeparatrixQ and added condition to KerrGeoOrbit to a…
Philip-Lynch Jul 17, 2023
92e1100
FourVelocity.m now uses functions from OrbitalFrequencies.m #50
Philip-Lynch Jul 17, 2023
1f6492a
Simplified KerrGeoFourVelocity API by removing index option. All comp…
Philip-Lynch Jul 17, 2023
29ddd77
Using ?PossibleZeroQ and closing secitons #50
Philip-Lynch Aug 2, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
168 changes: 103 additions & 65 deletions Kernel/FourVelocity.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,55 +9,86 @@


BeginPackage["KerrGeodesics`FourVelocity`",
{"KerrGeodesics`ConstantsOfMotion`"}];
{"KerrGeodesics`ConstantsOfMotion`",
"KerrGeodesics`OrbitalFrequencies`"}];

KerrGeoFourVelocity::usage = "KerrGeoVelocity[a,p,e,x] returns the four-velocity components as parametrized functions.";

Begin["`Private`"];


(* ::Subsection::Closed:: *)
barrywardell marked this conversation as resolved.
Show resolved Hide resolved
(* ::Subsection:: *)
(*Error messages*)


KerrGeoFourVelocity::parametrization = "Parameterization error: `1`"


(* ::Section::Closed:: *)
(*Kerr*)
(* ::Section:: *)
barrywardell marked this conversation as resolved.
Show resolved Hide resolved
(*Schwarzschild*)


(* ::Subsection:: *)
(*Generic (Mino)*)
(* ::Subsection::Closed:: *)
(*Circular, Equatorial*)


KerrGeoVelocityMino[(0|0.),p_,(0|0.),x_,initPhases_ ]:= Module[{En,L,Q,r,z,r1,r2,r3,r4,kr,zp,zm,kz, \[CapitalUpsilon]r, \[CapitalUpsilon]z,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use of (0|0.) misses the case where a or e are non-machine precision zero's (e.g. 0.``32) I suggest the following instead

KerrGeoVelocityMino[a_,p_,e_,x_,initPhases_ ] /; a==0 && e==0

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see in the rest of the package ?PossibleZeroQ is used. This should do the same thing, but I also think it's better style to be consistent throughout the package and use

KerrGeoVelocityMino[a_?PossibleZeroQ,p_,e_?PossibleZeroQ,x_,initPhases_ ] := ....

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PossibleZeroQ could be used here, although note the slightly different. PossibleZeroQ prints a message and returns True if a number might be zero but it's not sure, whereas the version with Equal will not return True in that case, which is equivalent to returning False for the purposes of the function definition. There are also circumstances where PossibleZeroQ will get it wrong (see the docs for example). These are all very much corner cases, so I think either PossibleZeroQ or Equal would be fine.

qr, qz, \[Lambda]local ,qr0, qz0, rprime, zprime, \[CapitalDelta], \[CapitalSigma], \[Omega], utContra,urContra,u\[Theta]Contra,uzContra,u\[Phi]Contra, utCo, urCo, u\[Theta]Co, u\[Phi]Co},

(*Constants of Motion*)
{En,L,Q}= {"\[ScriptCapitalE]","\[ScriptCapitalL]","\[ScriptCapitalQ]"}/.KerrGeoConstantsOfMotion[0,p,0,x];

\[CapitalUpsilon]z = p/Sqrt[-3+p];

{qr0,qz0} = initPhases;

qz[\[Lambda]_] := \[Lambda] \[CapitalUpsilon]z + qz0;



utContra= Function[{Global`\[Lambda]},Evaluate[Sqrt[p/(-3+p)] ], Listable];
urContra:= Function[{Global`\[Lambda]},Evaluate[0],Listable];
u\[Theta]Contra = Function[{Global`\[Lambda]}, Evaluate[(Sqrt[((1-x^2)/(-3+p))] Sin[qz[Global`\[Lambda]]] )/(p Sqrt[1+(-1+x^2) Cos[qz[Global`\[Lambda]]]^2])],Listable];
u\[Phi]Contra = Function[{Global`\[Lambda]},Evaluate[x/(Sqrt[-3+p] (p+p (-1+x^2) Cos[qz[Global`\[Lambda]]]^2))],Listable];

utCo = Function[{Global`\[Lambda]},Evaluate[-En], Listable];
urCo= Function[{Global`\[Lambda]},Evaluate[0],Listable];
u\[Theta]Co= Function[{Global`\[Lambda]},Evaluate[(p Sqrt[(1-x^2)/(-3+p)] Sin[qz[Global`\[Lambda]]])/ Sqrt[1+(-1+x^2) Cos[qz[Global`\[Lambda]]]^2]],Listable];
u\[Phi]Co= Function[{Global`\[Lambda]},Evaluate[L],Listable];

<|"\!\(\*SuperscriptBox[\(u\), \(t\)]\)"->utContra, "\!\(\*SuperscriptBox[\(u\), \(r\)]\)"->urContra, "\!\(\*SuperscriptBox[\(u\), \(\[Theta]\)]\)"-> u\[Theta]Contra, "\!\(\*SuperscriptBox[\(u\), \(\[Phi]\)]\)"-> u\[Phi]Contra,
"\!\(\*SubscriptBox[\(u\), \(t\)]\)"->utCo, "\!\(\*SubscriptBox[\(u\), \(r\)]\)"->urCo, "\!\(\*SubscriptBox[\(u\), \(\[Theta]\)]\)"-> u\[Theta]Co, "\!\(\*SubscriptBox[\(u\), \(\[Phi]\)]\)"-> u\[Phi]Co|>

(* ::Text:: *)
(*FixMe: Circular Equatorial Retrograde orbits don't normalize to -1*)


KerrGeoVelocityMino[a_,p_,e_,x_,initPhases_,index_ ]:= Module[{En,L,Q,r,z,r1,r2,r3,r4,kr,zp,zm,kz, \[CapitalUpsilon]r, \[CapitalUpsilon]z,
]


(* ::Section:: *)
barrywardell marked this conversation as resolved.
Show resolved Hide resolved
(*Kerr*)


(* ::Subsection::Closed:: *)
(*Generic (Mino)*)


KerrGeoVelocityMino[a_,p_,e_,x_,initPhases_]:= Module[{En,L,Q,r,z,r1,r2,r3,r4,kr,zp,zm,kz, \[CapitalUpsilon]r, \[CapitalUpsilon]z,
qr, qz, \[Lambda]local ,qr0, qz0, rprime, zprime, \[CapitalDelta], \[CapitalSigma], \[Omega], utContra,urContra,u\[Theta]Contra,uzContra,u\[Phi]Contra, utCo, urCo, u\[Theta]Co, u\[Phi]Co},

(*Constants of Motion*)
{En,L,Q}= {"\[ScriptCapitalE]","\[ScriptCapitalL]","\[ScriptCapitalQ]"}/.KerrGeoConstantsOfMotion[a,p,e,x];

(*Roots*)
r1 = p/(1-e);
r2 = p/(1+e);
zm = Sqrt[1-x^2];
{r1,r2,r3,r4}=KerrGeodesics`OrbitalFrequencies`Private`KerrGeoRadialRoots[a, p, e, x,En,Q];

(*Other Roots*)
r3 = 1/(1-En^2) - (r1 + r2)/2 + Sqrt[(-(1/(1-En^2)) + (r1 + r2)/2 )^2 - (a^2 Q)/(r1 r2 (1 - En^2))];
r4 = (a^2 Q)/(r1 r2 r3 (1-En^2));
{zp,zm}= KerrGeodesics`OrbitalFrequencies`Private`KerrGeoPolarRoots[a, p, e, x];

zp = Sqrt[a^2 (1 - En^2) + (L^2)/(1 - zm^2) ];
kr = ((r1-r2)(r3-r4))/((r1-r3)(r2-r4));

kz = a^2 (1-En^2) zm^2/zp^2;
Comment on lines 85 to 87
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can these be removed now? These are the lines that lead to divide-by-zero at the separatrix.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you mean the expressions for kr and kz?

These aren't singular at the separatrix. kz is finite and kr goes to 1. The problem is that EllipticK[kr] returns complex infinity when kr = 1 and that shows up in the frequencies, the analytic solutions for the coordinates and the four-velocities, with no obvious way (to me at least) for taking a finite limit exactly on the separatrix for generic orbits.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When there is a double root, you get special solutions. These are currently not explicitly implemented in the package (except circular orbits). We should add a check for r2==r3 and produce an error message.

Copy link
Member

@barrywardell barrywardell Aug 2, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I'm suggesting to remove the definitions of kr and kz. The problem with these is that in the case of repeated roots (eg. at the ISCO) they involve 0/0 and this causes problems. I had thought kr and kz are no longer used in this function, in which case we could remove them to avoid warning messages. Unfortunately I now see they are used further down. Could we define functions for these quantities also and then remove kr and kz here? The special handling of repeated roots can then be done in KerrGeoMinoFrequencyr and the other added functions.


(*Frequencies*)
\[CapitalUpsilon]r = \[Pi]/(2 EllipticK[kr]) Sqrt[(1 - En^2)(r1 - r3)(r2 - r4)];
\[CapitalUpsilon]z = (\[Pi] zp)/(2EllipticK[kz] );
\[CapitalUpsilon]r = KerrGeodesics`OrbitalFrequencies`Private`KerrGeoMinoFrequencyr[a,p,e,x,{En,L,Q},{r1,r2,r3,r4}];
\[CapitalUpsilon]z = KerrGeodesics`OrbitalFrequencies`Private`KerrGeoMinoFrequency\[Theta][a,p,e,x,{En,L,Q},{zp,zm}];

(*Action Angle Phases*)
{ qr0, qz0} = {initPhases[[1]], initPhases[[2]]};
Expand All @@ -79,28 +110,25 @@
\[Omega][qr_] := Sqrt[r[qr]^2+ a^2];
\[CapitalSigma][qr_,qz_] := r[qr]^2 + a^2 z[qz]^2;

If[index == "Contravariant",

utContra= Function[{Global`\[Lambda]},Evaluate[1/\[CapitalSigma][qr[Global`\[Lambda]],qz[Global`\[Lambda]]] (\[Omega][qr[Global`\[Lambda]]]^2/\[CapitalDelta][qr[Global`\[Lambda]]] ( \[Omega][qr[Global`\[Lambda] ]]^2 En - a L) - a^2 (1-z[qz[Global`\[Lambda]]]^2)En + a L)], Listable];
urContra:= Function[{Global`\[Lambda]},Evaluate[( rprime[qr[Global`\[Lambda]]] \[CapitalUpsilon]r)/\[CapitalSigma][qr[Global`\[Lambda]],qz[Global`\[Lambda]]]],Listable];
u\[Theta]Contra = Function[{Global`\[Lambda]}, Evaluate[-(\[CapitalUpsilon]z zprime[qz[Global`\[Lambda]]])/(\[CapitalSigma][qr[Global`\[Lambda]],qz[Global`\[Lambda]]]Sqrt[1-z[qz[Global`\[Lambda]]]^2])],Listable];
u\[Phi]Contra = Function[{Global`\[Lambda]},Evaluate[1/\[CapitalSigma][qr[Global`\[Lambda]],qz[Global`\[Lambda]]] (a/\[CapitalDelta][qr[Global`\[Lambda]]] ( \[Omega][qr[Global`\[Lambda]]]^2 En - a L) - a En + L/(1-z[qz[Global`\[Lambda]]]^2))],Listable];

<|"\!\(\*SuperscriptBox[\(u\), \(t\)]\)"->utContra, "\!\(\*SuperscriptBox[\(u\), \(r\)]\)"->urContra, "\!\(\*SuperscriptBox[\(u\), \(\[Theta]\)]\)"-> u\[Theta]Contra, "\!\(\*SuperscriptBox[\(u\), \(\[Phi]\)]\)"-> u\[Phi]Contra|>,

(*Else if Index \[Equal] Covariant*)

utCo = Function[{Global`\[Lambda]},Evaluate[-En], Listable];
urCo= Function[{Global`\[Lambda]},Evaluate[( rprime[qr[Global`\[Lambda]]] \[CapitalUpsilon]r)/\[CapitalDelta][qr[Global`\[Lambda]]]],Listable];
u\[Theta]Co= Function[{Global`\[Lambda]},Evaluate[-((\[CapitalUpsilon]z zprime[qz[Global`\[Lambda]]])/Sqrt[1-z[qz[Global`\[Lambda]]]^2])],Listable];
u\[Phi]Co= Function[{Global`\[Lambda]},Evaluate[L],Listable];

<|"\!\(\*SubscriptBox[\(u\), \(t\)]\)"->utCo, "\!\(\*SubscriptBox[\(u\), \(r\)]\)"->urCo, "\!\(\*SubscriptBox[\(u\), \(\[Theta]\)]\)"-> u\[Theta]Co, "\!\(\*SubscriptBox[\(u\), \(\[Phi]\)]\)"-> u\[Phi]Co|>
<|"\!\(\*SuperscriptBox[\(u\), \(t\)]\)"->utContra, "\!\(\*SuperscriptBox[\(u\), \(r\)]\)"->urContra,
"\!\(\*SuperscriptBox[\(u\), \(\[Theta]\)]\)"-> u\[Theta]Contra, "\!\(\*SuperscriptBox[\(u\), \(\[Phi]\)]\)"-> u\[Phi]Contra,
"\!\(\*SubscriptBox[\(u\), \(t\)]\)"->utCo, "\!\(\*SubscriptBox[\(u\), \(r\)]\)"->urCo,
"\!\(\*SubscriptBox[\(u\), \(\[Theta]\)]\)"-> u\[Theta]Co, "\!\(\*SubscriptBox[\(u\), \(\[Phi]\)]\)"-> u\[Phi]Co|>
]


]


(* ::Subsection:: *)
(*Equatorial (Darwin)*)
Expand All @@ -110,46 +138,53 @@
(*Circular Case*)


KerrGeoVelocityDarwin[a_,p_,(0|0.),x_,initPhases_,index_ ]:= Module[{ut,ur,u\[Theta],u\[Phi], MinoVelocities,ut1,ur1,u\[Theta]1,u\[Phi]1},
KerrGeoVelocityDarwin[a_,p_,(0|0.),x_,initPhases_]:= Module[{ut,ur,u\[Theta],u\[Phi], MinoVelocities,utContra,urContra,u\[Theta]Contra,u\[Phi]Contra,
barrywardell marked this conversation as resolved.
Show resolved Hide resolved
utCo,urCo,u\[Theta]Co,u\[Phi]Co,utUp,urUp,u\[Theta]Up,u\[Phi]Up, utDown,urDown,u\[Theta]Down,u\[Phi]Down},

MinoVelocities = KerrGeoVelocityMino[a,p,0,x,{0,0}, index];
MinoVelocities = KerrGeoVelocityMino[a,p,0,x,{0,0}];

utUp="\!\(\*SuperscriptBox[\(u\), \(t\)]\)"; urUp="\!\(\*SuperscriptBox[\(u\), \(r\)]\)";
u\[Theta]Up="\!\(\*SuperscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]Up="\!\(\*SuperscriptBox[\(u\), \(\[Phi]\)]\)";
utDown="\!\(\*SubscriptBox[\(u\), \(t\)]\)"; urDown="\!\(\*SubscriptBox[\(u\), \(r\)]\)";
u\[Theta]Down="\!\(\*SubscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]Down="\!\(\*SubscriptBox[\(u\), \(\[Phi]\)]\)";

If[index == "Contravariant",
ut1="\!\(\*SuperscriptBox[\(u\), \(t\)]\)"; ur1="\!\(\*SuperscriptBox[\(u\), \(r\)]\)"; u\[Theta]1="\!\(\*SuperscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]1="\!\(\*SuperscriptBox[\(u\), \(\[Phi]\)]\)";,
ut1="\!\(\*SubscriptBox[\(u\), \(t\)]\)"; ur1="\!\(\*SubscriptBox[\(u\), \(r\)]\)"; u\[Theta]1="\!\(\*SubscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]1="\!\(\*SubscriptBox[\(u\), \(\[Phi]\)]\)";
];
(*All components are Constants*)
ut = Function[{Global`\[Chi]},Evaluate[MinoVelocities [ut1][Global`\[Chi]]], Listable];
ur = Function[{Global`\[Chi]},Evaluate[0],Listable];
u\[Theta]= Function[{Global`\[Chi]},Evaluate[0],Listable];
u\[Phi]= Function[{Global`\[Chi]},Evaluate[MinoVelocities [u\[Phi]1][Global`\[Chi]]],Listable];
utContra = Function[{Global`\[Chi]},Evaluate[MinoVelocities [utUp][Global`\[Chi]]],Listable];
urContra = Function[{Global`\[Chi]},Evaluate[0],Listable];
u\[Theta]Contra = Function[{Global`\[Chi]}, Evaluate[0],Listable];
u\[Phi]Contra = Function[{Global`\[Chi]}, Evaluate[MinoVelocities [u\[Phi]Up][Global`\[Chi]]],Listable];

utCo = Function[{Global`\[Chi]},Evaluate[MinoVelocities [utDown][Global`\[Chi]]],Listable];
urCo = Function[{Global`\[Chi]},Evaluate[0],Listable];
u\[Theta]Co = Function[{Global`\[Chi]}, Evaluate[0],Listable];
u\[Phi]Co = Function[{Global`\[Chi]}, Evaluate[MinoVelocities [u\[Phi]Down][Global`\[Chi]]],Listable];

<|ut1-> ut, ur1-> ur, u\[Theta]1-> u\[Theta], u\[Phi]1-> u\[Phi] |>
<|utUp ->utContra, urUp ->urContra,
u\[Theta]Up-> u\[Theta]Contra, u\[Phi]Up-> u\[Phi]Contra,
utDown->utCo, urDown->urCo,
u\[Theta]Down-> u\[Theta]Co, u\[Phi]Down-> u\[Phi]Co|>

]


(* ::Subsubsection:: *)
(* ::Subsubsection::Closed:: *)
(*Eccentric Case*)


KerrGeoVelocityDarwin[a_,p_,e_,x_/;x^2==1,initPhases_,index_ ]:= Module[{En,L,Q,r,z,r1,r2,r3,r4,kr, \[CapitalUpsilon]r, \[CapitalLambda]r,yr,\[Lambda]0r,r01,\[CapitalLambda]r1,\[Lambda],
\[Chi]0,\[Nu], \[Chi]local ,qr0, qz0, rprime, zprime, \[CapitalDelta], \[CapitalSigma], \[Omega], ut,ur,u\[Theta],u\[Phi], MinoVelocities,ut1,ur1,u\[Theta]1,u\[Phi]1},
KerrGeoVelocityDarwin[a_,p_,e_,x_/;x^2==1,initPhases_]:= Module[{En,L,Q,r,z,r1,r2,r3,r4,kr, \[CapitalUpsilon]r, \[CapitalLambda]r,yr,\[Lambda]0r,r01,\[CapitalLambda]r1,\[Lambda],
\[Chi]0,\[Nu], \[Chi]local ,qr0, qz0, rprime, zprime, \[CapitalDelta], \[CapitalSigma], \[Omega], ut,ur,u\[Theta],u\[Phi], MinoVelocities,utContra,urContra,u\[Theta]Contra,u\[Phi]Contra,
utCo,urCo,u\[Theta]Co,u\[Phi]Co,utUp,urUp,u\[Theta]Up,u\[Phi]Up, utDown,urDown,u\[Theta]Down,u\[Phi]Down},

(*Constants of Motion*)
{En,L,Q}= {"\[ScriptCapitalE]","\[ScriptCapitalL]","\[ScriptCapitalQ]"}/.KerrGeoConstantsOfMotion[a,p,e,x];

(*Roots*)
r1 = p/(1-e);
r2 = p/(1+e);
{r1,r2,r3,r4}=KerrGeodesics`OrbitalFrequencies`Private`KerrGeoRadialRoots[a, p, e, x,En,Q];

(*Other Roots*)
r3 = 1/(1-En^2) - (r1 + r2)/2 + Sqrt[(-(1/(1-En^2)) + (r1 + r2)/2 )^2 - (a^2 Q)/(r1 r2 (1 - En^2))];
r4 = (a^2 Q)/(r1 r2 r3 (1-En^2));
kr = ((r1-r2)(r3-r4))/((r1-r3)(r2-r4));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove this line?

Copy link
Collaborator Author

@Philip-Lynch Philip-Lynch Aug 2, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need the expression for kr in the function definition on line 198, so we can't just get rid of it.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the same issue as above.


(*Frequencies*)
\[CapitalUpsilon]r = \[Pi]/(2 EllipticK[kr]) Sqrt[(1 - En^2)(r1 - r3)(r2 - r4)];
\[CapitalUpsilon]r = KerrGeodesics`OrbitalFrequencies`Private`KerrGeoMinoFrequencyr[a,p,e,x,{En,L,Q},{r1,r2,r3,r4}];

(*Initial Phase*)
\[Chi]0 = initPhases[[1]];
Expand All @@ -167,47 +202,50 @@
\[Lambda][\[Nu]_]:=\[CapitalLambda]r Floor[\[Nu]/(2\[Pi])]+If[Mod[\[Nu],2\[Pi]]<=\[Pi], \[Lambda]0r[r[\[Nu]]]-\[CapitalLambda]r1,\[CapitalLambda]r-\[Lambda]0r[r[\[Nu]]]];


MinoVelocities = KerrGeoVelocityMino[a,p,e,x,{0,0}, index];
MinoVelocities = KerrGeoVelocityMino[a,p,e,x,{0,0}];

If[index == "Contravariant",
ut1="\!\(\*SuperscriptBox[\(u\), \(t\)]\)"; ur1="\!\(\*SuperscriptBox[\(u\), \(r\)]\)"; u\[Theta]1="\!\(\*SuperscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]1="\!\(\*SuperscriptBox[\(u\), \(\[Phi]\)]\)";,
ut1="\!\(\*SubscriptBox[\(u\), \(t\)]\)"; ur1="\!\(\*SubscriptBox[\(u\), \(r\)]\)"; u\[Theta]1="\!\(\*SubscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]1="\!\(\*SubscriptBox[\(u\), \(\[Phi]\)]\)";
];
utUp="\!\(\*SuperscriptBox[\(u\), \(t\)]\)"; urUp="\!\(\*SuperscriptBox[\(u\), \(r\)]\)";
u\[Theta]Up="\!\(\*SuperscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]Up="\!\(\*SuperscriptBox[\(u\), \(\[Phi]\)]\)";
utDown="\!\(\*SubscriptBox[\(u\), \(t\)]\)"; urDown="\!\(\*SubscriptBox[\(u\), \(r\)]\)";
u\[Theta]Down="\!\(\*SubscriptBox[\(u\), \(\[Theta]\)]\)"; u\[Phi]Down="\!\(\*SubscriptBox[\(u\), \(\[Phi]\)]\)";

ut = Function[{Global`\[Chi]},Evaluate[MinoVelocities [ut1][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];
ur = Function[{Global`\[Chi]},Evaluate[MinoVelocities [ur1][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];
u\[Theta] = Function[{Global`\[Chi]}, Evaluate[0],Listable];
u\[Phi] = Function[{Global`\[Chi]}, Evaluate[MinoVelocities [u\[Phi]1][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];
utContra = Function[{Global`\[Chi]},Evaluate[MinoVelocities [utUp][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];
urContra = Function[{Global`\[Chi]},Evaluate[MinoVelocities [urUp][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];
u\[Theta]Contra = Function[{Global`\[Chi]}, Evaluate[0],Listable];
u\[Phi]Contra = Function[{Global`\[Chi]}, Evaluate[MinoVelocities [u\[Phi]Up][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];

<|ut1-> ut, ur1-> ur, u\[Theta]1-> u\[Theta], u\[Phi]1-> u\[Phi] |>
utCo = Function[{Global`\[Chi]},Evaluate[MinoVelocities [utDown][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];
urCo = Function[{Global`\[Chi]},Evaluate[MinoVelocities [urDown][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];
u\[Theta]Co = Function[{Global`\[Chi]}, Evaluate[0],Listable];
u\[Phi]Co = Function[{Global`\[Chi]}, Evaluate[MinoVelocities [u\[Phi]Down][\[Lambda][Global`\[Chi]-\[Chi]0]]],Listable];

<|utUp ->utContra, urUp ->urContra,
u\[Theta]Up-> u\[Theta]Contra, u\[Phi]Up-> u\[Phi]Contra,
utDown->utCo, urDown->urCo,
u\[Theta]Down-> u\[Theta]Co, u\[Phi]Down-> u\[Phi]Co|>

]


(* ::Section::Closed:: *)
(* ::Section:: *)
(*KerrGeoFourVelocity Wrapper*)


Options[KerrGeoFourVelocity] = {"Covariant" -> False, "Parametrization"-> "Mino"}
Options[KerrGeoFourVelocity] = {"Parametrization"-> "Mino"}
SyntaxInformation[KerrGeoFourVelocity] = {"ArgumentsPattern"->{_,_,_,_,OptionsPattern[]}};


KerrGeoFourVelocity[a_,p_,e_,x_,initPhases:{_,_}:{0,0}, OptionsPattern[]]:= Module[{param, index},
KerrGeoFourVelocity[a_,p_,e_,x_,initPhases:{_,_}:{0,0}, OptionsPattern[]]:= Module[{param},
param = OptionValue["Parametrization"];

If[OptionValue["Covariant"], index = "Covariant" , index="Contravariant", Message[KerrGeoFourVelocity::opttf,"Covariant",OptionValue["Covariant"]]; Return[] ];


If[param == "Darwin",

If[ Abs[x]!=1,
Message[KerrGeoFourVelocity::parametrization, "Darwin parameterization only valid for equatorial motion"];
Return[];,
Return[KerrGeoVelocityDarwin[a,p,e,x,initPhases, index]]]];
Return[KerrGeoVelocityDarwin[a,p,e,x,initPhases]]]];


If[param == "Mino", Return[KerrGeoVelocityMino[a,p,e,x,initPhases, index]]];
If[param == "Mino", Return[KerrGeoVelocityMino[a,p,e,x,initPhases]]];

Message[KerrGeoFourVelocity::parametrization, "Unrecognized Paramaterization: " <> param];

Expand Down
79 changes: 41 additions & 38 deletions Kernel/KerrGeoOrbit.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@
KerrGeoOrbit::OutOfBounds = "For this hyperbolic orbit the Darwin parameter \[Chi] must be between `1` and `2`"


KerrGeoOrbit::OnSeparatrix = "This orbit is on the separatrix where many expressions are numerically singular. Aborting..."


(* ::Subsection::Closed:: *)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be better if we could fix the KerrGeoMinoFrequency* functions to work in this case.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, but from my preliminary look at this, there doesn't seem to be an easy finite limit that can be taken. These frequencies may need to be completely re-derived for this case. Hence the next best thing of stopping the user from using the functions on the separatrix until these are known.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are there common sub-cases (e.g. ISCO) where we could make it work, but leave the generic case for another time?

(*Being the private context*)

Expand Down Expand Up @@ -133,7 +136,7 @@
]


(* ::Section::Closed:: *)
(* ::Section:: *)
barrywardell marked this conversation as resolved.
Show resolved Hide resolved
(*Kerr*)


Expand Down Expand Up @@ -374,7 +377,7 @@
(* Hopper, Forseth, Osburn, and Evans, PRD 92 (2015)*)


(* ::Subsubsection:: *)
(* ::Subsubsection::Closed:: *)
(*Main file that calculates geodesics using spectral integration*)


Expand Down Expand Up @@ -597,7 +600,7 @@



(* ::Subsubsection:: *)
(* ::Subsubsection::Closed:: *)
(*Scattering orbit (e > 1)*)


Expand Down Expand Up @@ -1236,41 +1239,41 @@
SyntaxInformation[KerrGeoOrbit] = {"ArgumentsPattern"->{_,_,OptionsPattern[]}};


KerrGeoOrbit[a_,p_,e_,x_, initPhases:{_,_,_,_}:{0,0,0,0},OptionsPattern[]]:=Module[{param, method},
(*FIXME: add stability check but make it possible to turn it off*)

method = OptionValue["Method"];
param = OptionValue["Parametrization"];

If[param == "Darwin" && Abs[x]!=1, Message[KerrGeoOrbit::parametrization, "Darwin parameterization only valid for equatorial motion"]; Return[];];

If[Precision[{a,p,e,x}] > 30, method = "Analytic"];
If[e > 1, method = "Analytic"];

If[method == "FastSpec",

If[param == "Mino", If[PossibleZeroQ[a] || PossibleZeroQ[e], Return[KerrGeoOrbitMino[a, p, e, x, initPhases]], Return[KerrGeoOrbitFastSpec[a, p, e, x, initPhases]]]];
If[param == "Darwin",
If[PossibleZeroQ[a], Return[KerrGeoOrbitSchwarzDarwin[p, e]], Return[KerrGeoOrbitFastSpecDarwin[a,p,e,x,initPhases]]]
];
Message[KerrGeoOrbit::parametrization, "Unrecognized parametrization: " <> OptionValue["Parametrization"]];

];

If[method == "Analytic",
(*Changed "KerrGeoOrbitDarwin" to "KerrGeoOrbitEquatorialDarwin"*)
If[param == "Mino", Return[KerrGeoOrbitMino[a, p, e, x, initPhases]]];
If[param == "Phases", Return[KerrGeoOrbitPhases[a, p, e, x]]];
If[param == "Darwin",
If[PossibleZeroQ[a], Return[KerrGeoOrbitSchwarzDarwin[p, e]], Return[KerrGeoOrbitEquatorialDarwin[a,p,e,x,initPhases]]]
];
Message[KerrGeoOrbit::parametrization, "Unrecognized parametrization: " <> OptionValue["Parametrization"]];

];

Message[KerrGeoOrbit::general, "Method " <> method <> " is not one of {FastSpec, Analytic}"];

]
KerrGeoOrbit[a_, p_, e_, x_, initPhases : {_, _, _, _} : {0, 0, 0, 0}, OptionsPattern[]] := Module[{param, method},
(*FIXME: add stability check but make it possible to turn it off*)
If[a!=0 ||e != 0,If[KerrGeodesics`SpecialOrbits`Private`KerrGeoOnSeparatrixQ[a, p, e, x],Message[KerrGeoOrbit::OnSeparatrix];Abort[];]];
method = OptionValue["Method"];
param = OptionValue["Parametrization"];
If[param == "Darwin" && Abs[x] != 1, Message[KerrGeoOrbit::parametrization, "Darwin parameterization only valid for equatorial motion"]; Return[];];
If[Precision[{a, p, e, x}] > 30, method = "Analytic"];
If[e > 1, method = "Analytic"];
If[method == "FastSpec",
If[param == "Mino", If[PossibleZeroQ[a] || PossibleZeroQ[e], Return[KerrGeoOrbitMino[a, p, e, x, initPhases]], Return[KerrGeoOrbitFastSpec[a, p, e, x, initPhases]]]];
If[param == "Darwin",
If[PossibleZeroQ[a], Return[KerrGeoOrbitSchwarzDarwin[p, e]], Return[KerrGeoOrbitFastSpecDarwin[a, p, e, x, initPhases]]]
];
Message[KerrGeoOrbit::parametrization, "Unrecognized parametrization: " <> OptionValue["Parametrization"]];
];
If[method == "Analytic",
(*Changed "KerrGeoOrbitDarwin" to "KerrGeoOrbitEquatorialDarwin"*)
If[param == "Mino", Return[KerrGeoOrbitMino[a, p, e, x, initPhases]]];
If[param == "Phases", Return[KerrGeoOrbitPhases[a, p, e, x]]];
If[param == "Darwin",
If[PossibleZeroQ[a], Return[KerrGeoOrbitSchwarzDarwin[p, e]], Return[KerrGeoOrbitEquatorialDarwin[a, p, e, x, initPhases]]]
];
Message[KerrGeoOrbit::parametrization, "Unrecognized parametrization: " <> OptionValue["Parametrization"]];
];
Message[KerrGeoOrbit::general, "Method " <> method <> " is not one of {FastSpec, Analytic}"];
]
Comment on lines +1242 to +1276
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Lots of unrelated whitespace changes here...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure why these all changed. I didn't manually add these in so it must have been automatic.



KerrGeoOrbitFunction /:
Expand Down
Loading