From ab16e47e0bd336d25ab2f0f4afcc66d45ee8878c Mon Sep 17 00:00:00 2001 From: Philip-Lynch Date: Mon, 17 Jul 2023 15:53:12 +0200 Subject: [PATCH] Created KerrGeoOnSeparatrixQ and added condition to KerrGeoOrbit to abort when on the separatrix --- Kernel/KerrGeoOrbit.m | 89 +++++++++++++++++++------------------ Kernel/OrbitalFrequencies.m | 2 +- Kernel/SpecialOrbits.m | 44 +++++++++++++----- 3 files changed, 79 insertions(+), 56 deletions(-) diff --git a/Kernel/KerrGeoOrbit.m b/Kernel/KerrGeoOrbit.m index 1d1731b..1123ef9 100644 --- a/Kernel/KerrGeoOrbit.m +++ b/Kernel/KerrGeoOrbit.m @@ -19,13 +19,16 @@ 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*) @@ -33,7 +36,7 @@ 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,7 +1231,7 @@ ] -(* ::Section::Closed:: *) +(* ::Section:: *) (*KerrGeoOrbit and KerrGeoOrbitFuction*) @@ -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}"]; + + ] KerrGeoOrbitFunction /: diff --git a/Kernel/OrbitalFrequencies.m b/Kernel/OrbitalFrequencies.m index d0e6ee6..f6f51c9 100644 --- a/Kernel/OrbitalFrequencies.m +++ b/Kernel/OrbitalFrequencies.m @@ -70,7 +70,7 @@ ] -(* ::Section::Closed:: *) +(* ::Section:: *) (*Orbital Frequencies*) diff --git a/Kernel/SpecialOrbits.m b/Kernel/SpecialOrbits.m index aa8eb5a..1c4ce16 100644 --- a/Kernel/SpecialOrbits.m +++ b/Kernel/SpecialOrbits.m @@ -8,7 +8,7 @@ (*Define usage for public functions*) -(* ::Section::Closed:: *) +(* ::Section:: *) (*Create Package*) @@ -16,7 +16,7 @@ {"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*)