Skip to content

Commit

Permalink
FourVelocity.m now uses functions from OrbitalFrequencies.m #50
Browse files Browse the repository at this point in the history
  • Loading branch information
Philip-Lynch committed Jul 17, 2023
1 parent ab16e47 commit 92e1100
Show file tree
Hide file tree
Showing 4 changed files with 29 additions and 46 deletions.
43 changes: 13 additions & 30 deletions Kernel/FourVelocity.m
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,8 @@


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

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

Expand All @@ -27,7 +28,7 @@
(*Schwarzschild*)


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


Expand Down Expand Up @@ -66,45 +67,31 @@
]


(* ::Subsection:: *)
(*Eccentric*)


(* ::Section:: *)
(*Kerr*)


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


(* ::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,
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];

(*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;

(*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 Down Expand Up @@ -153,7 +140,7 @@
(*Equatorial (Darwin)*)


(* ::Subsubsection:: *)
(* ::Subsubsection::Closed:: *)
(*Circular Case*)


Expand All @@ -176,7 +163,7 @@
]


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


Expand All @@ -187,16 +174,12 @@
{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];

(*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));

(*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 Down Expand Up @@ -232,7 +215,7 @@
]


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


Expand Down
10 changes: 5 additions & 5 deletions Kernel/KerrGeoOrbit.m
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
KerrGeoOrbitFunction::usage = "KerrGeoOrbitFunction[a,p,e,x,assoc] an object for storing the trajectory and orbital parameters in the assoc Association.";


(* ::Subsection:: *)
(* ::Subsection::Closed:: *)
(*Error messages*)


Expand All @@ -36,7 +36,7 @@
Begin["`Private`"];


(* ::Subsection:: *)
(* ::Subsection::Closed:: *)
(*Error messages*)


Expand Down Expand Up @@ -515,7 +515,7 @@



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


Expand Down Expand Up @@ -723,7 +723,7 @@



(* ::Subsection:: *)
(* ::Subsection::Closed:: *)
(*Generic (Fast Spec - Mino)*)


Expand Down Expand Up @@ -1231,7 +1231,7 @@
]


(* ::Section:: *)
(* ::Section::Closed:: *)
(*KerrGeoOrbit and KerrGeoOrbitFuction*)


Expand Down
8 changes: 4 additions & 4 deletions Kernel/OrbitalFrequencies.m
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
Begin["`Private`"];


(* ::Section::Closed:: *)
(* ::Section:: *)
(*Roots of the radial and polar equations*)


Expand Down Expand Up @@ -121,11 +121,11 @@
KerrGeoProperFrequencyFactor[a_?PossibleZeroQ ,p_,e_,x_]:=(p^2 ((1+e) (28+4 e^2+(-12+p) p)-((1+e) (-4+p) (-6+2 e+p) EllipticE[(4 e)/(-6+2 e+p)]+2 (6+2 e-p) (3+e^2-p) EllipticPi[(2 e (-4+p))/((1+e) (-6+2 e+p)),(4 e)/(-6+2 e+p)])/EllipticK[(4 e)/(-6+2 e+p)]))/(2 (-1+e) (1+e)^2 (-4+p)^2)


(* ::Subsection::Closed:: *)
(* ::Subsection:: *)
(*Kerr*)


(* ::Subsubsection::Closed:: *)
(* ::Subsubsection:: *)
(*KerrGeoMinoFrequencyr*)


Expand All @@ -147,7 +147,7 @@



(* ::Subsubsection::Closed:: *)
(* ::Subsubsection:: *)
(*KerrGeoMinoFrequency\[Theta]*)


Expand Down
14 changes: 7 additions & 7 deletions Kernel/SpecialOrbits.m
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
{"KerrGeodesics`ConstantsOfMotion`"}];


(* ::Subsection:: *)
(* ::Subsection::Closed:: *)
(*Usage messages*)


Expand Down Expand Up @@ -81,7 +81,7 @@
];


(* ::Section:: *)
(* ::Section::Closed:: *)
(*Photon Sphere*)


Expand Down Expand Up @@ -254,7 +254,7 @@
KerrGeoISSO[a_,x_]:=KerrGeoSeparatrix[a,0,x]


(* ::Section:: *)
(* ::Section::Closed:: *)
(*Bound Orbit Q*)


Expand All @@ -264,7 +264,7 @@
]


(* ::Section:: *)
(* ::Section::Closed:: *)
(*Scatter Orbit Q*)


Expand All @@ -275,7 +275,7 @@
KerrGeoScatterOrbitQ[a_?NumericQ, p_?NumericQ, e_?NumericQ, x_?NumericQ] := If[p >= KerrGeoSeparatrix[a,e,x] && e >= 1, True, False]


(* ::Section:: *)
(* ::Section::Closed:: *)
(*Plunge Orbit Q*)


Expand All @@ -287,7 +287,7 @@
If[KerrGeoBoundOrbitQ[0,p,e,1] == KerrGeoScatterOrbitQ[0,p,e,1] == False, True, False]


(* ::Section:: *)
(* ::Section::Closed:: *)
(*OnSeparatrixQ*)


Expand Down Expand Up @@ -356,7 +356,7 @@
]


(* ::Section:: *)
(* ::Section::Closed:: *)
(*Resonances*)


Expand Down

0 comments on commit 92e1100

Please sign in to comment.