Created KerrGeoOnSeparatrixQ and added condition to KerrGeoOrbit to a…

…bort when on the separatrix

Kernel/KerrGeoOrbit.m

@@ -19,21 +19,24 @@
 KerrGeoOrbitFunction::usage = "KerrGeoOrbitFunction[a,p,e,x,assoc] an object for storing the trajectory and orbital parameters in the assoc Association.";
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*Error messages*)
 
 
 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:: *)
 (*Being the private context*)
 
 
 Begin["`Private`"];
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*Error messages*)
 
 
@@ -133,7 +136,7 @@
 ]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Kerr*)
 
 
@@ -374,7 +377,7 @@
 (* Hopper, Forseth, Osburn, and Evans, PRD 92 (2015)*)
 
 
-(* ::Subsubsection:: *)
+(* ::Subsubsection::Closed:: *)
 (*Main file that calculates geodesics using spectral integration*)
 
 
@@ -512,7 +515,7 @@
 
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*Generic (Mino)*)
 
 
@@ -597,7 +600,7 @@
 
 
 
-(* ::Subsubsection:: *)
+(* ::Subsubsection::Closed:: *)
 (*Scattering orbit (e > 1)*)
 
 
@@ -720,7 +723,7 @@
 
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*Generic (Fast Spec - Mino)*)
 
 
@@ -1228,49 +1231,49 @@
 ]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*KerrGeoOrbit and KerrGeoOrbitFuction*)
 
 
 Options[KerrGeoOrbit] = {"Parametrization" -> "Mino", "Method" -> "FastSpec"}
 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}"];
+   
+   ]
 
 
 KerrGeoOrbitFunction /:

Kernel/OrbitalFrequencies.m

@@ -70,7 +70,7 @@
 ]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Orbital Frequencies*)
 
 

Kernel/SpecialOrbits.m

@@ -8,15 +8,15 @@
 (*Define usage for public functions*)
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Create Package*)
 
 
 BeginPackage["KerrGeodesics`SpecialOrbits`",
 	{"KerrGeodesics`ConstantsOfMotion`"}];
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*Usage messages*)
 
 
@@ -34,7 +34,9 @@
 
 (*KerrGeoBoundOrbitQ::usage = "KerrGeoBoundOrbitQ[a,p,e,x] tests if the orbital parameters correspond to a bound orbit."
 KerrGeoScatterOrbitQ::usage = "KerrGeoScatterOrbitQ[a,p,e,x] tests if the orbital parameters correspond to a scatter orbit."
-KerrGeoPlungeOrbitQ::usage = "KerrGeoPlungeOrbitQ[a,p,e,x] tests if the orbital parameters correspond to a plunge orbit."*)
+KerrGeoPlungeOrbitQ::usage = "KerrGeoPlungeOrbitQ[a,p,e,x] tests if the orbital parameters correspond to a plunge orbit."
+KerrGeoOnSeparatrixQ::usage = "KerrGeoOnSeparatrixQ[a,p,e,x] tests if the orbital parameters correspond to being on the separatrix."
+*)
 
 
 (* ::Subsection::Closed:: *)
@@ -134,7 +136,7 @@
 ]
 
 
-(* ::Section:: *)
+(* ::Section::Closed:: *)
 (*Innermost bound spherical orbits (IBSO)*)
 
 
@@ -176,7 +178,7 @@
 p/.FindRoot[IBSOPoly/.{a->a1,x->x1},{p,KerrGeoIBSO[a1,0],KerrGeoIBSO[a1,-1]},WorkingPrecision->Max[MachinePrecision,prec-1]]];
 
 
-(* ::Section:: *)
+(* ::Section::Closed:: *)
 (*Separatrix*)
 
 
@@ -252,7 +254,7 @@
 KerrGeoISSO[a_,x_]:=KerrGeoSeparatrix[a,0,x]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Bound Orbit Q*)
 
 
@@ -262,7 +264,7 @@
 ]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Scatter Orbit Q*)
 
 
@@ -273,7 +275,7 @@
 KerrGeoScatterOrbitQ[a_?NumericQ, p_?NumericQ, e_?NumericQ, x_?NumericQ] := If[p >= KerrGeoSeparatrix[a,e,x] && e >= 1, True, False]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Plunge Orbit Q*)
 
 
@@ -285,6 +287,24 @@
 	If[KerrGeoBoundOrbitQ[0,p,e,1] == KerrGeoScatterOrbitQ[0,p,e,1] == False, True, False]
 
 
+(* ::Section:: *)
+(*OnSeparatrixQ*)
+
+
+(* ::Text:: *)
+(*Test to see if the orbit is exactly on the separatrix, where many of our expressions become singular. This is a stopgap solution until we can think of a nice way of taking the separatrix limit of all of these functions. *)
+
+
+KerrGeoOnSeparatrixQ[a_?NumericQ,p_?NumericQ,e_?NumericQ,x_?NumericQ]:= Module[{En,L,Q,r1,r2,r3,r4,prec},
+
+prec = Precision[{a,p,e,x}];
+{En,L,Q} = Values[KerrGeoConstantsOfMotion[a,p,e,x]];
+{r1,r2,r3,r4} = KerrGeodesics`OrbitalFrequencies`Private`KerrGeoRadialRoots[a, p, e, x, En, Q];
+(*Assuming KerrGeoSeparatrix only finds solution to within half of the working precision*)
+Abs[r2 - r3] <=  10^(-prec/2)  
+]
+
+
 (* ::Section::Closed:: *)
 (*Orbit type*)
 
@@ -336,11 +356,11 @@
 ]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Resonances*)
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*r\[Theta]-resonances*)
 
 
@@ -562,7 +582,7 @@
 ];
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*Generic resonance interface*)
 
 
@@ -582,7 +602,7 @@
 ]
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Close the package*)