diff --git a/Kernel/KerrGeoOrbit.m b/Kernel/KerrGeoOrbit.m index d67143b..4810091 100644 --- a/Kernel/KerrGeoOrbit.m +++ b/Kernel/KerrGeoOrbit.m @@ -12,14 +12,12 @@ {"KerrGeodesics`ConstantsOfMotion`", "KerrGeodesics`OrbitalFrequencies`", "KerrGeodesics`FourVelocity`"}]; - KerrGeoOrbit::usage = "KerrGeoOrbit[a,p,e,x] returns a KerrGeoOrbitFunction[..] which stores the orbital trajectory and parameters."; KerrGeoOrbitFunction::usage = "KerrGeoOrbitFunction[a,p,e,x,assoc] an object for storing the trajectory and orbital parameters in the assoc Association."; - Begin["`Private`"]; -(* ::Section::Closed:: *) +(* ::Section:: *) (*Schwarzschild*) @@ -85,11 +83,11 @@ ] -(* ::Section::Closed:: *) +(* ::Section:: *) (*Kerr*) -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*Equatorial (Darwin)*) @@ -161,14 +159,14 @@ -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*Equatorial (Fast Spec - Darwin)*) (* Hopper, Forseth, Osburn, and Evans, PRD 92 (2015)*) -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Subroutines that checks for the number of samples necessary for spectral integration*) @@ -234,7 +232,7 @@ ]; -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Main file that calculates geodesics using spectral integration*) @@ -311,14 +309,14 @@ ] -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*Circular (Fast Spec - Darwin)*) (* Hopper, Forseth, Osburn, and Evans, PRD 92 (2015)*) -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Main file that calculates geodesics using spectral integration*) @@ -401,7 +399,7 @@ ] -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*Circular (Mino)*) @@ -450,7 +448,7 @@ -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*Generic (Mino)*) @@ -565,14 +563,14 @@ ] -(* ::Subsection::Closed:: *) +(* ::Subsection:: *) (*Generic (Fast Spec - Mino)*) (* Hopper, Forseth, Osburn, and Evans, PRD 92 (2015)*) -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Subroutines for calculating \[Lambda](\[Psi]) and \[Lambda](\[Chi])*) @@ -622,7 +620,7 @@ -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Subroutine for calculating \[Psi](\[Lambda]) and \[Chi](\[Lambda])*) @@ -635,7 +633,7 @@ ]; -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Subroutines for calculating Subscript[\[CapitalDelta]\[Phi], r](Subscript[q, r]), Subscript[\[CapitalDelta]\[Phi], \[Theta]](Subscript[q, \[Theta]]), Subscript[\[CapitalDelta]t, r](Subscript[q, r]), Subscript[\[CapitalDelta]t, \[Theta]](Subscript[q, \[Theta]])*) @@ -739,7 +737,7 @@ ]; -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Generic subroutines that transform functions from *) (*\[Lambda] dependence to \[Psi] or \[Chi] dependence*) @@ -791,7 +789,7 @@ ]; -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Subroutines that checks for the number of samples necessary for spectral integration of an even function*) @@ -868,7 +866,7 @@ ]; -(* ::Subsubsection::Closed:: *) +(* ::Subsubsection:: *) (*Subroutine that performs spectral integration on even functions*) @@ -1070,7 +1068,7 @@ ] -(* ::Section::Closed:: *) +(* ::Section:: *) (*KerrGeoOrbit and KerrGeoOrbitFuction*) diff --git a/Kernel/KerrGeoPlunge.m b/Kernel/KerrGeoPlunge.m new file mode 100644 index 0000000..6604d2b --- /dev/null +++ b/Kernel/KerrGeoPlunge.m @@ -0,0 +1,563 @@ +(* ::Package:: *) + +(* ::Title:: *) +(*KerrGeoOrbit subpackage of KerrGeodesics*) + + +(* ::Chapter:: *) +(*Define usage for public functions*) + + +BeginPackage["KerrGeodesics`KerrGeoPlunge`", + {"KerrGeodesics`ConstantsOfMotion`", + "KerrGeodesics`OrbitalFrequencies`", + "KerrGeodesics`SpecialOrbits`", + "KerrGeodesics`FourVelocity`"}]; + + +KerrGeoPlunge::usage = "Takes either KerrGeoPlunge[a, {En,\[ScriptCapitalL],\[ScriptCapitalQ]}] or KerrGeoPlunge[a, ISSO , RI] for generic Plunges or ISSO plunges respectively and returns a KerrGeoPlungeFunction[..] which stores the orbital trajectory and parameters. Here the ISSO plunges are paramaterised in terms of the choice of the radius of the ISSO for a given a there are range of allowed RI's which correspond to differeing inclinations in the prograde and retrograde directions, if an RI outside this range is given the used is provided with the range of allowed values to re-run the funciton."; +KerrGeoPlungeFunction::usage = "KerrGeoPlungeFunction[a,rI,assoc] an object for storing the trajectory and orbital parameters in the assoc Association."; +Begin["`Private`"]; + + +(* ::Section:: *) +(*Kerr*) + + +(* ::Subsection:: *) +(*ISSO Plunges*) + + +KerrGeoPlunge::r0outofbounds = "Intial radius `1` is not between ISSO radius `2`, and inner radial root `3`."; +KerrGeoPlunge::\[Theta]0outofbounds = "Intial polar angle `1` must be between the ranges (`2`,`3`) or (-`2`,-`3`)"; + + +KerrGeoISSOPlunge[a_, PlungeType_ ,Arg_, initCoords_] := Module[ + {consts, \[Theta]0, z, r0,Mino, z0, \[Phi]0, initConditions, t0, assoc,M=1,\[Chi]INC, Minoz, MinozFunc, MinoRFunc,\[CapitalUpsilon]z,MCMP,MCPP,MCPM,MCMM,DN,RISSOMIN,RISSOMAX ,t ,RI , r, \[Theta], \[Phi], \[ScriptCapitalE], \[ScriptCapitalL], \[ScriptCapitalQ], R4, RM, RP,t\[Delta], kz, Z1, Z2, J,velocity,MinoR,tr,tz,\[Phi]r,\[Phi]z,LRoot}, + + RM = 1-Sqrt[1-a^2]; + RP = 1+Sqrt[1-a^2]; + + DN = (27-45 a^2+17 a^4+a^6+8 a^3 (1-a^2))^(1/3); + RISSOMIN= 3+Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]-1/2 \[Sqrt](72+8 (-6+a^2)-(4 (9-10 a^2+a^4))/DN-4 DN+(64 a^2)/Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]); + RISSOMAX=3+Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]+1/2 \[Sqrt](72+8 (-6+a^2)-(4 (9-10 a^2+a^4))/DN-4 DN+(64 a^2)/Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]); + + (*This finds the correct sign of the angular momentum by determining the root of where the non-analytic form of the equation touches the axis on a \[ScriptCapitalL](RI) plot*) + If[PlungeType == "ISSORadialParam", + RI = Arg; + If[a==0,\[ScriptCapitalQ]=0,\[ScriptCapitalQ] = -M (RI)^(5/2) ((Sqrt[(RI-RP)(RI-RM)]-2Sqrt[RI])^2-4a^2)/(4a^2 (RI^(3/2)-Sqrt[RI]-Sqrt[(RI-RP)(RI-RM)])); + ]; + \[ScriptCapitalE] = Sqrt[a^2 \[ScriptCapitalQ]-2 RI^3+3 RI^4]/(Sqrt[3] RI^2); + \[ScriptCapitalL] = Sqrt[3 a^2 \[ScriptCapitalQ]-a^2 RI^2-\[ScriptCapitalQ] RI^2+3 RI^4+a^2 RI^2 \[ScriptCapitalE]^2-3 RI^4 \[ScriptCapitalE]^2]/RI; + RootCheckFunc[y_]:= (a^6-3 a^2 y^(5/2) (3 y^(3/2)+6 Sqrt[a^2+(-2+y) y]-2 y Sqrt[a^2+(-2+y) y])+y^(9/2) (20 Sqrt[y]-11 y^(3/2)+5 Sqrt[a^2+(-2+y) y]+3 y Sqrt[a^2+(-2+y) y])+a^4 (-4 y+3 y^2+Sqrt[y] Sqrt[a^2+(-2+y) y]+3 y^(3/2) Sqrt[a^2+(-2+y) y])); + If[a==0,LRoot=6,y/.FindRoot[RootCheckFunc[y],{y,6}]]; + If[a!=0,If[RI>LRoot ,\[ScriptCapitalL]=-\[ScriptCapitalL]]]; + ]; + + If[PlungeType == "ISSOIncParam", + \[Chi]INC = Cos[Arg - \[Pi]/2]; + RI = KerrGeoISSO[a,\[Chi]INC]; + {\[ScriptCapitalE],\[ScriptCapitalL],\[ScriptCapitalQ]} = Values[KerrGeoConstantsOfMotion[a,RI,0,\[Chi]INC]][[1;;3]] + ]; + J = 1-\[ScriptCapitalE]^2; + If[a==0, {Z1,Z2}={0,\[ScriptCapitalL]}, {Z1,Z2}= {Sqrt[1/2 (1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J)-Sqrt[(1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J))^2-(4 \[ScriptCapitalQ])/(a^2 J)])],Sqrt[ (a^2 J)/2 (1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J)+Sqrt[(1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J))^2-(4 \[ScriptCapitalQ])/(a^2 J)])]}]; + R4 = (a^2 \[ScriptCapitalQ])/(J*RI^3); + kz = a*Sqrt[J](Z1/Z2); + \[CapitalUpsilon]z= (\[Pi]*Z2)/(2*EllipticK[kz^2]); + If[Arg==-(\[Pi]/2), \[ScriptCapitalL]= -\[ScriptCapitalL]]; + If[a<0,\[ScriptCapitalL]=-\[ScriptCapitalL]]; + consts = <|"\[ScriptCapitalE]"-> \[ScriptCapitalE],"\[ScriptCapitalL]"-> \[ScriptCapitalL], "\[ScriptCapitalQ]"->\[ScriptCapitalQ]|>; + If[initCoords===Automatic, + {t0,r0,\[Theta]0,\[Phi]0}={0,R4,\[Pi]/2,0}, + {t0,r0,\[Theta]0,\[Phi]0}=initCoords + ]; + If[r0>RI||r0\[Pi] - Abs[ArcCos[Z1]], + Message[KerrGeoPlunge::\[Theta]0outofbounds,\[Theta]0,Abs[ArcCos[Z1]],\[Pi]-Abs[ArcCos[Z1]] ]; + Return[$Failed] + ]; + z0 = Cos[\[Theta]0]; + Minoz[z_] :=InverseJacobiSN[z/Z1,kz^2]/Z2; + If[a==0, Minoz[z_]:=0]; + If[PlungeType == "ISSOIncParam", + If[Abs[Arg]==\[Pi]/2, Minoz[z_] :=0]; + MinozFunc=Function[{Global`z}, Evaluate[Minoz[Global`z]-Minoz[z0] ],Listable];]; + + If[PlungeType == "ISSORadialParam", + If[Round[Arg - RISSOMIN,10^-14]==0, Minoz[z_] :=0]; + MinozFunc=Function[{Global`z}, Evaluate[Minoz[Global`z]-Minoz[z0] ],Listable]; + If[Round[Arg - RISSOMAX,10^-14]==0,Minoz[z_] :=0]; + MinozFunc=Function[{Global`z}, Evaluate[Minoz[Global`z]-Minoz[z0] ],Listable];]; + MinoR[r_] :=-(2 Sqrt[r-R4])/Sqrt[J (RI-r) (R4-RI)^2]; + MinoRFunc=Function[{Global`r}, Evaluate[MinoR[Global`r]-MinoR[r0] ],Listable]; + + r[\[Lambda]_] := ((RI (RI-R4)^2 J*\[Lambda]^2+4*R4)/((RI-R4)^2 J*\[Lambda]^2+4)); + z[\[Lambda]_]:= Z1*JacobiSN[Z2*\[Lambda] ,kz^2]; + If[a!=0,tr[\[Lambda]_]:= ((a^2+RI^2) (-a \[ScriptCapitalL]+(a^2+RI^2) \[ScriptCapitalE])(\[Lambda]) )/((RI-RM) (RI-RP))+(2 ((R4-RI)^2) \[ScriptCapitalE] (\[Lambda]) )/(4+J (R4-RI)^2 (\[Lambda]) ^2)-(R4+3 RI+2 (RM+RP))/Sqrt[J] \[ScriptCapitalE] ArcTan[((\[Lambda]) (RI-R4)Sqrt[J ] )/2 ]+( (a^2+RM^2) (-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RM^2 \[ScriptCapitalE])Log[Sqrt[(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RM]+2 Sqrt[RM-R4])^2/(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RM]-2 Sqrt[RM-R4])^2]])/(Sqrt[RM-R4] (RI-RM)^(3/2) (-RM+RP)Sqrt[J])+( (a^2+RP^2) (-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RP^2 \[ScriptCapitalE]) Log[Sqrt[(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]+2 Sqrt[RP-R4])^2/(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]-2 Sqrt[RP-R4])^2]])/(Sqrt[J]Sqrt[RP-R4] (RI-RP)^(3/2) (RM-RP))]; + If[a==0,tr[\[Lambda]_]:= ((a^2+RI^2) (-a \[ScriptCapitalL]+(a^2+RI^2) \[ScriptCapitalE])(\[Lambda]) )/((RI-RM) (RI-RP))+(2 ((R4-RI)^2) \[ScriptCapitalE] (\[Lambda]) )/(4+J (R4-RI)^2 (\[Lambda]) ^2)-(R4+3 RI+2 (RM+RP))/Sqrt[J] \[ScriptCapitalE] ArcTan[((\[Lambda]) (RI-R4)Sqrt[J ] )/2 ]+( (a^2+RP^2) (-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RP^2 \[ScriptCapitalE]) Log[Sqrt[(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]+2 Sqrt[RP-R4])^2/(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]-2 Sqrt[RP-R4])^2]])/(Sqrt[J]Sqrt[RP-R4] (RI-RP)^(3/2) (RM-RP))]; + If[a!=0,\[Phi]r[\[Lambda]_]:= 1/Sqrt[J] 2 a (((\[Lambda]) (RI-R4)Sqrt[J ] (-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RI^2 \[ScriptCapitalE]))/( 2(-R4+RI) (RI-RM) (RI-RP))+((-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RM^2 \[ScriptCapitalE])Log[Sqrt[(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RM]+2 Sqrt[RM-R4])^2/(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RM]-2 Sqrt[RM-R4])^2]])/(2Sqrt[RM-R4] (RI-RM)^(3/2) (-RM+RP))+((-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RP^2 \[ScriptCapitalE])Log[Sqrt[(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]+2 Sqrt[RP-R4])^2/(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]-2 Sqrt[RP-R4])^2]])/(2Sqrt[RP-R4] (RI-RP)^(3/2) (RM-RP))) ] ; + If[a==0,\[Phi]r[\[Lambda]_]:= 1/Sqrt[J] 2 a (((\[Lambda]) (RI-R4)Sqrt[J ] (-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RI^2 \[ScriptCapitalE]))/( 2(-R4+RI) (RI-RM) (RI-RP))+((-a \[ScriptCapitalL]+a^2 \[ScriptCapitalE]+RP^2 \[ScriptCapitalE])Log[Sqrt[(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]+2 Sqrt[RP-R4])^2/(\[Lambda] (RI-R4)Sqrt[J ] Sqrt[RI-RP]-2 Sqrt[RP-R4])^2]])/(2Sqrt[RP-R4] (RI-RP)^(3/2) (RM-RP))) ] ; + tz[\[Lambda]_]:= 1/J \[ScriptCapitalE] (-Z2 EllipticE[JacobiAmplitude[Z2 \[Lambda] ,kz^2],kz^2]+(Z2-a^2 J/Z2) Z2 \[Lambda]); + \[Phi]z[\[Lambda]_]:= \[ScriptCapitalL]/Z2 EllipticPi[Z1^2,JacobiAmplitude[Z2 \[Lambda] ,kz^2],kz^2]; + + t=Function[{Global`\[Lambda]}, Evaluate[ a*\[ScriptCapitalL] Global`\[Lambda] + tr[Global`\[Lambda]+ MinoR[r0]] + tz[Global`\[Lambda]+ Minoz[z0]]-tr[MinoR[r0]]-tz[Minoz[z0]] + t0], Listable]; + r=Function[{Global`\[Lambda]}, Evaluate[ r[Global`\[Lambda] + MinoR[r0]] ], Listable]; + \[Theta]=Function[{Global`\[Lambda]}, Evaluate[ ArcCos[z[Global`\[Lambda] + Minoz[z0]]]] , Listable]; + \[Phi]=Function[{Global`\[Lambda]}, Evaluate[-a*\[ScriptCapitalE] Global`\[Lambda] + \[Phi]r[Global`\[Lambda]+ MinoR[r0]] + \[Phi]z[Global`\[Lambda]+ Minoz[z0]] - \[Phi]r[MinoR[r0]] - \[Phi]z[Minoz[z0]] + \[Phi]0], Listable]; + MCPP = -Abs[MinoR[RP]] -MinoR[r0]; + MCMP = -Abs[MinoR[RM]] -MinoR[r0]; + MCMM = Abs[MinoR[RM]] - MinoR[r0]; + MCPM = Abs[MinoR[RP]] - MinoR[r0]; + assoc = Association[ + "a" -> a, + "rI" -> RI, + "Parametrization"->"Mino", + "Energy" -> \[ScriptCapitalE], + "AngularMomentum" -> \[ScriptCapitalL], + "CarterConstant" -> \[ScriptCapitalQ], + "ConstantsOfMotion" -> consts, + "RadialFrequency"-> 0, + "PolarFrequency"-> \[CapitalUpsilon]z, + "Frequencies"-> <|"\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(r\)]\)"-> 0, "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Theta]\)]\)" -> \[CapitalUpsilon]z|>, + "RadialPeriod"-> Infinity, + "PolarPeriod"-> 2\[Pi]/\[CapitalUpsilon]z, + "Periods"-> <|"\!\(\*SubscriptBox[\(\[Lambda]\), \(r\)]\)"-> Infinity, "\!\(\*SubscriptBox[\(\[Lambda]\), \(\[Theta]\)]\)" -> 2\[Pi]/\[CapitalUpsilon]z|>, + "RadialRoots"-> {RI,RI,RI,R4}, + "EllipticBasis" -> <|"\!\(\*SubscriptBox[\(k\), \(z\)]\)"-> kz|>, + "RadialMinoTime"-> MinoR, + "AllowedISSORangeForGivenSpin"-> <|"\!\(\*SubscriptBox[\(r\), \(IMIN\)]\)"-> RISSOMIN, "\!\(\*SubscriptBox[\(r\), \(IMAX\)]\)" -> RISSOMAX|>, + "RadialRootCrossingTimeMino"-> <| + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(I\)], \(\)]\)"-> \[PlusMinus] Infinity(* Time to ISSO *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(inner\)], \(\)]\)"-> Abs[MinoR[R4]] - MinoR[r0](* Time to Inner*)|>, + "HorizonCrossingTimeMino"-> <| + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(+\)], \(+\)]\)"-> MCPP(* Time ingoing branch crosses outer horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(-\)], \(+\)]\)"-> MCMP(* Time ingoing branch crosses inner horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(-\)], \(-\)]\)"-> MCMM(* Time outgoing branch crosses inner horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(+\)], \(-\)]\)"-> MCPM(* Time outgoing branch crosses outer horizon *) + |>, + "PolarRoots"-> {Z1,Z2}, + "Trajectory" -> {t,r,\[Theta],\[Phi]} + ]; + + KerrGeoPlungeFunction[a, \[ScriptCapitalE], \[ScriptCapitalL], \[ScriptCapitalQ], assoc] +] + + + +(* ::Subsection:: *) +(*Real Roots Generic Plunge (Mino)*) + + +KerrGeoPlunge::r0outofboundsGen = "Intial radius `1` is not between the outer radial root `2`, and inner radial root `3`."; +KerrGeoPlunge::\[Theta]0outofbounds = "Intial polar angle `1` must be between the ranges (`2`,`3`) or (-`2`,-`3`)"; + + +Options[KerrGeoRealPlungeMino] = {"Roots"-> Automatic}; +KerrGeoRealPlungeMino[a_, \[ScriptCapitalE]_, \[ScriptCapitalL]_, \[ScriptCapitalQ]_ , initCoords_, OptionsPattern[Options[KerrGeoRealPlungeMino]]] := Module[ + {consts,assoc,M,Mino,t0,r0,\[Phi]0,kr,hR,hP,hM,\[CapitalUpsilon]t,\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]tr,\[CapitalUpsilon]tz,\[CapitalUpsilon]z,\[CapitalUpsilon]r,\[CapitalUpsilon]\[Phi]r,\[CapitalUpsilon]\[Phi]z,t,r,z,z0,\[Theta]0, \[Theta], J,\[Phi], R1,R2,R3,R4,MCMP,MCPP,MCPM,MCMM, Minoz, MinozFunc, MinoRFunc , RM, RP,ROOTS,RealRoots,ComplexRoots, kz, Z1, Z2, MinoR}, + + M=1; + J = 1-\[ScriptCapitalE]^2; + RM = M-Sqrt[M^2-a^2]; + RP = M+Sqrt[M^2-a^2]; + consts = <|"\[ScriptCapitalE]"-> \[ScriptCapitalE],"\[ScriptCapitalL]"-> \[ScriptCapitalL], "\[ScriptCapitalQ]"->\[ScriptCapitalQ]|>; + + If[ + OptionValue["Roots"]===Automatic, + ROOTS = Sort[r/.Solve[(\[ScriptCapitalE](r^2+a^2)-a*\[ScriptCapitalL])^2-(r^2-2M*r+a^2)(r^2+(a*\[ScriptCapitalE]-\[ScriptCapitalL])^2+\[ScriptCapitalQ])==0,r]], + ROOTS = OptionValue["Roots"] + ]; + If[Length[Select[ROOTS,#>=RP&]]>=3, + R1= ROOTS[[4]]; + R2= ROOTS[[3]]; + R3= ROOTS[[2]]; + R4= ROOTS[[1]];]; + + If[Length[Select[ROOTS,#>=RP&]]<3, + R1= ROOTS[[2]]; + R2= ROOTS[[1]]; + R3= ROOTS[[4]]; + R4= ROOTS[[3]];]; + + + + If[initCoords===Automatic, + {t0,r0,\[Theta]0,\[Phi]0}={0,R4,\[Pi]/2,0}, + {t0,r0,\[Theta]0,\[Phi]0}=initCoords + ]; + If[r0R3, + Message[KerrGeoPlunge::r0outofboundsGen,r0,R3,R4]; + Return[$Failed] + ]; + + + z0 = Cos[\[Theta]0]; + + Z1=Sqrt[1/2 (1+\[ScriptCapitalL]^2/(a^2 (1-\[ScriptCapitalE]^2))+\[ScriptCapitalQ]/(a^2 (1-\[ScriptCapitalE]^2))-Sqrt[-((4 \[ScriptCapitalQ])/(a^2 (1-\[ScriptCapitalE]^2)))+(-1-\[ScriptCapitalL]^2/(a^2 (1-\[ScriptCapitalE]^2))-\[ScriptCapitalQ]/(a^2 (1-\[ScriptCapitalE]^2)))^2])]; + Z2 =Sqrt[ (a^2 (1-\[ScriptCapitalE]^2))/2 (1+\[ScriptCapitalL]^2/(a^2 (1-\[ScriptCapitalE]^2))+\[ScriptCapitalQ]/(a^2 (1-\[ScriptCapitalE]^2))+Sqrt[-((4 \[ScriptCapitalQ])/(a^2 (1-\[ScriptCapitalE]^2)))+(-1-\[ScriptCapitalL]^2/(a^2 (1-\[ScriptCapitalE]^2))-\[ScriptCapitalQ]/(a^2 (1-\[ScriptCapitalE]^2)))^2])]; + kz = a^2*J(Z1^2/Z2^2); + kr = ((R3-R4)(R1-R2))/((R3-R1)(R4-R2)); + + If[Abs[\[Theta]0]\[Pi] - Abs[ArcCos[Z1]], + Message[KerrGeoPlunge::\[Theta]0outofbounds,\[Theta]0,Abs[ArcCos[Z1]],\[Pi]-Abs[ArcCos[Z1]] ]; + Return[$Failed] + ]; + + + + hR = (R3-R4)/(R3-R1); + hP = hR (R1-RP)/(R4-RP); + hM = hR (R1-RM)/(R4-RM); + + Minoz[z_] :=InverseJacobiSN[z/Z1,kz^2]/Z2; + MinozFunc=Function[{Global`z}, Evaluate[Minoz[Global`z]-Minoz[z0] ],Listable]; + MinoR[r_]:= 2 InverseJacobiSN[Sqrt[(-R1+R3) (r-R4)]/Sqrt[(r-R1) (R3-R4)] ,kr]/Sqrt[J (R1-R3) (R2-R4)]; + MinoRFunc=Function[{Global`r}, Evaluate[MinoR[Global`r]- MinoR[r0]],Listable]; + + If[Abs[Arg]==\[Pi]/2, Minoz[z_] :=0]; + MinozFunc=Function[{Global`z}, Evaluate[Minoz[Global`z]-Minoz[z0] ],Listable]; + + If[\[ScriptCapitalQ]==0, Minoz[z_] :=0]; + MinozFunc=Function[{Global`z}, Evaluate[Minoz[Global`z]-Minoz[z0] ],Listable]; + + + \[CapitalUpsilon]r= \[Pi]/(2*EllipticK[kr]) Sqrt[J(R3-R1)(R4-R2)]; + \[CapitalUpsilon]z= (\[Pi]*Z2)/(2*EllipticK[kz]); + + \[CapitalUpsilon]\[Phi]r = a/(RP-RM) ((2\[ScriptCapitalE]*RP-a*\[ScriptCapitalL])/(R1-RP) (1-(R4-R1)/(R4-RP) (EllipticPi[hP,kr]/EllipticK[kr]) )- (2\[ScriptCapitalE]*RM-a*\[ScriptCapitalL])/(R1-RM) (1-(R4-R1)/(R4-RM) (EllipticPi[hM,kr]/EllipticK[kr])) ); + \[CapitalUpsilon]\[Phi]z = \[ScriptCapitalL]/EllipticK[kz] EllipticPi[Z1^2,kz]; + \[CapitalUpsilon]\[Phi]= \[CapitalUpsilon]\[Phi]r+\[CapitalUpsilon]\[Phi]z; + + \[Phi]Tr[\[Xi]_]:= (-2*a*\[ScriptCapitalE]*(R4-R1))/((RP-RM)Sqrt[J*(R3-R1)(R4-R2)])*((2*RP-a*\[ScriptCapitalL]/\[ScriptCapitalE])/((R4-RP)(R1-RP)) EllipticPi[hP,\[Xi],kr]-(2*RM-a*\[ScriptCapitalL]/\[ScriptCapitalE])/((R4-RM)(R1-RM)) EllipticPi[hM,\[Xi],kr]); + \[Phi]Tz[\[Xi]_]:= -\[ScriptCapitalL]/Z2 EllipticPi[Z1^2,\[Xi],kz]; + + \[Phi]r[\[Lambda]_] := \[Phi]Tr[JacobiAmplitude[EllipticK[kr]*(\[CapitalUpsilon]r*\[Lambda])/\[Pi],kr]] - \[Phi]Tr[\[Pi]]/(2*\[Pi]) (\[CapitalUpsilon]r*\[Lambda]); + \[Phi]z[\[Lambda]_]:= -\[Phi]Tz[JacobiAmplitude[EllipticK[kz]*(2*(\[CapitalUpsilon]z*\[Lambda]))/\[Pi],kz]] + \[Phi]Tz[\[Pi]]/\[Pi] (\[CapitalUpsilon]z*\[Lambda]); + + + \[CapitalUpsilon]tr = (4+a^2)\[ScriptCapitalE]+\[ScriptCapitalE](1/2 ((4+R3+R4+R1)R1-R3*R4+(R3-R1)(R4-R2) EllipticE[kr]/EllipticK[kr]+(4+R3+R4+R1+R2)(R4-R1) EllipticPi[hR,kr]/EllipticK[kr])+2/(RP-RM) (((4-a*\[ScriptCapitalL]/\[ScriptCapitalE])RP-2*a^2)/(R1-RP) (1-(R4-R1)/(R4-RP) EllipticPi[hP,kr]/EllipticK[kr])-((4-a*\[ScriptCapitalL]/\[ScriptCapitalE])RM-2*a^2)/(R1-RM) (1-(R4-R1)/(R4-RM) EllipticPi[hM,kr]/EllipticK[kr]))); + \[CapitalUpsilon]tz = -a^2*\[ScriptCapitalE]+(\[ScriptCapitalE]*\[ScriptCapitalQ])/(J(Z1^2)) (1-EllipticE[kz]/EllipticK[kz]); + If[\[ScriptCapitalQ]==0, \[CapitalUpsilon]tz = -a^2*\[ScriptCapitalE]]; + \[CapitalUpsilon]t= \[CapitalUpsilon]tr +\[CapitalUpsilon]tz; + + + tTr[\[Xi]_]:= (\[ScriptCapitalE](R4-R1))/Sqrt[J(R3-R1)(R4-R2)] ((4+R3+R4+R1+R2)EllipticPi[hR,\[Xi],kr]-4/(RP-RM) ((RP(4-a*\[ScriptCapitalL]/\[ScriptCapitalE])-2*a^2)/((R4-RP)(R1-RP)) EllipticPi[hP,\[Xi],kr]-(RM(4-a*\[ScriptCapitalL]/\[ScriptCapitalE])-2*a^2)/((R4-RM)(R1-RM)) EllipticPi[hM,\[Xi],kr])+((R3-R1)(R4-R2))/(R4-R1) (EllipticE[\[Xi],kr]-(hR*Sin[\[Xi]]*Cos[\[Xi]]Sqrt[1-kr*(Sin[\[Xi]])^2])/(1-hR*(Sin[\[Xi]])^2))); + tTz[\[Xi]_]:= -\[ScriptCapitalE]/J Z2*EllipticE[\[Xi],kz]; + + r[\[Lambda]_]:=(R1(R3-R4)(JacobiSN[EllipticK[kr]/\[Pi]* (\[CapitalUpsilon]r*\[Lambda]),kr])^2-R4(R3-R1))/((R3-R4)(JacobiSN[EllipticK[kr]/\[Pi]* (\[CapitalUpsilon]r*\[Lambda]),kr])^2-(R3-R1)); + z[\[Lambda]_]:= Z1*JacobiSN[EllipticK[kz]*2 (\[CapitalUpsilon]z*\[Lambda])/\[Pi],kz] ; + + + tr[\[Lambda]_] := tTr[JacobiAmplitude[EllipticK[kr]*(\[CapitalUpsilon]r*\[Lambda])/\[Pi],kr]] - tTr[\[Pi]]/(2*\[Pi]) (\[CapitalUpsilon]r*\[Lambda]); + tz[\[Lambda]_]:= tTz[JacobiAmplitude[EllipticK[kz]*(2*(\[CapitalUpsilon]z*\[Lambda]))/\[Pi],kz]] - tTz[\[Pi]]/\[Pi] ( \[CapitalUpsilon]z*\[Lambda]); + + t=Function[{Global`\[Lambda]}, Evaluate[\[CapitalUpsilon]t*Global`\[Lambda] + tr[Global`\[Lambda]+ MinoR[r0]] + tz[Global`\[Lambda]+ Minoz[z0]]-tr[MinoR[r0]]-tz[Minoz[z0]] + t0], Listable]; + r=Function[{Global`\[Lambda]}, Evaluate[ r[Global`\[Lambda] + MinoR[r0]] ], Listable]; + \[Theta]=Function[{Global`\[Lambda]}, Evaluate[ ArcCos[z[Global`\[Lambda] + Minoz[z0]]]] , Listable]; + \[Phi]=Function[{Global`\[Lambda]}, Evaluate[\[CapitalUpsilon]\[Phi]*Global`\[Lambda]+ \[Phi]r[Global`\[Lambda]+ MinoR[r0]] + \[Phi]z[Global`\[Lambda]+ Minoz[z0]] - \[Phi]r[MinoR[r0]] - \[Phi]z[Minoz[z0]] + \[Phi]0], Listable]; + + MCPP = -Abs[MinoR[RP]] -MinoR[r0]+(2\[Pi])/\[CapitalUpsilon]r; + MCMP = -Abs[MinoR[RM]] -MinoR[r0]+(2\[Pi])/\[CapitalUpsilon]r; + MCMM = Abs[MinoR[RM]] - MinoR[r0]; + MCPM = Abs[MinoR[RP]] - MinoR[r0]; + + assoc = Association[ + "a" -> a, + "Parametrization"->"Mino", + "Energy" -> \[ScriptCapitalE], + "AngularMomentum" -> \[ScriptCapitalL], + "CarterConstant" -> \[ScriptCapitalQ], + "ConstantsOfMotion" -> consts, + "RadialFrequency"-> \[CapitalUpsilon]r, + "PolarFrequency"-> \[CapitalUpsilon]z, + "Frequencies"-> <| "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(r\)]\)"->\[CapitalUpsilon]r, "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Theta]\)]\)" -> \[CapitalUpsilon]z |>, + "RadialPeriod"-> 2\[Pi]/\[CapitalUpsilon]r, + "PolarPeriod"-> 2\[Pi]/\[CapitalUpsilon]z, + "Periods"-> <|"\!\(\*SubscriptBox[\(\[Lambda]\), \(r\)]\)"-> 2\[Pi]/\[CapitalUpsilon]r, "\!\(\*SubscriptBox[\(\[Lambda]\), \(\[Theta]\)]\)" -> 2\[Pi]/\[CapitalUpsilon]z|>, + "EllipticBasis" -> <|"\!\(\*SubscriptBox[\(k\), \(r\)]\)"-> kr,"\!\(\*SubscriptBox[\(k\), \(z\)]\)"-> kz|>, + "RadialRoots"-> {R4,R3,R2,R1}, + "RadialMinoTime"-> MinoRFunc, + "RadialRootCrossingTimeMino"-> <| + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(Inner\)], \(\)]\)"-> Abs[MinoR[R4]] - MinoR[r0](* Time outgoing to inner root*), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(Outer\)], \(\)]\)"-> Abs[MinoR[R3]] - MinoR[r0](* Time outgoing to outer root*)|>, + "HorizonCrossingTimeMino"-> <| + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(+\)], \(+\)]\)"-> MCPP(* Time ingoing branch crosses outer horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(-\)], \(+\)]\)"-> MCMP(* Time ingoing branch crosses inner horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(-\)], \(-\)]\)"-> MCMM(* Time outgoing branch crosses inner horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(+\)], \(-\)]\)"-> MCPM(* Time outgoing branch crosses outer horizon *) + |>, + "PolarRoots"-> {Z1,Z2}, + "Trajectory" -> {t,r,\[Theta],\[Phi]} + ]; + KerrGeoPlungeFunction[a, \[ScriptCapitalE], \[ScriptCapitalL], \[ScriptCapitalQ], assoc] +] + + + +(* ::Subsection:: *) +(*Complex Roots Generic Plunge (Mino)*) + + +(* ::Text:: *) +(*FIXME: make the initial phases work in this case*) + + +KerrGeoPlunge::r0outofboundsGen = "Intial radius `1` is not between the outer radial root `2`, and inner radial root `3`."; +KerrGeoPlunge::\[Theta]0outofbounds = "Intial polar angle `1` must be between the ranges (`2`,`3`) or (-`2`,-`3`)"; + + +Options[KerrGeoComplexPlungeMino]= {"Roots"->Automatic}; +KerrGeoComplexPlungeMino[a_, \[ScriptCapitalE]_, \[ScriptCapitalL]_, \[ScriptCapitalQ]_ , initCoords_, OptionsPattern[Options[KerrGeoComplexPlungeMino]]] := Module[{consts,MCPP, MCPM, MCMP , MCMM , assoc,MinoR,\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],t0,r0,\[Phi]0,M,D1M,D1P,D2M,D2P,Minoz,MinozFunc,MinoRFunc,\[Theta]0,z0,e,b,c,d,A,B,chi,kr,p2,f, t, r,z, \[Theta], J,\[Phi], R1,R2,R3,R4, RM, RP,ROOTS,RealRoots,ComplexRoots,kz,Z1, Z2, AMR,CNR,SNR,DNR,AMZ}, + + M=1; + J = 1-\[ScriptCapitalE]^2; + RM = M-Sqrt[M^2-a^2]; + RP = M+Sqrt[M^2-a^2]; + + ROOTS = r/.Solve[(\[ScriptCapitalE](r^2+a^2)-a*\[ScriptCapitalL])^2-(r^2-2M*r+a^2)(r^2+(a*\[ScriptCapitalE]-\[ScriptCapitalL])^2+\[ScriptCapitalQ])==0,r]; + RealRoots = Select[ROOTS,Im[#]==0&]; + ComplexRoots = Select[ROOTS,Im[#]!=0&]; + consts = <|"\[ScriptCapitalE]"-> \[ScriptCapitalE],"\[ScriptCapitalL]"-> \[ScriptCapitalL], "\[ScriptCapitalQ]"->\[ScriptCapitalQ]|>; + R1= RealRoots[[1]]; + R2= RealRoots[[2]]; + R3=ComplexRoots[[1]]; + R4= ComplexRoots[[2]]; + + If[a==0, {Z1,Z2}={0,\[ScriptCapitalL]}, {Z1,Z2}= {Sqrt[1/2 (1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J)-Sqrt[(1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J))^2-(4 \[ScriptCapitalQ])/(a^2 J)])],Sqrt[ (a^2 J)/2 (1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J)+Sqrt[(1+(\[ScriptCapitalL]^2+\[ScriptCapitalQ])/(a^2 J))^2-(4 \[ScriptCapitalQ])/(a^2 J)])]}]; + + e = R2; + b = R1; + c = Re[R3]; + d = Abs[Im[R4]]; + A = Sqrt[(e-c)^2+d^2]; + B = Sqrt[(b-c)^2+d^2]; + chi = (A*b+e*B)/(A+B); + kr = Sqrt[((e-b)^2-(A-B)^2)/(4*A*B)]; + p2 = b*A^2+e*B^2-(e+b)*A*B; + f = (4 A B)/(A-B)^2; + kz = a*Sqrt[J](Z1/Z2); + + \[CapitalUpsilon]r = 2 \[Pi] Sqrt[A B (1-\[ScriptCapitalE]^2) ]/(4 EllipticK[kr^2]); + \[CapitalUpsilon]\[Theta] = \[Pi]/2 EllipticK[kz^2]/Z2; + + + D1M=Sqrt[ -f]; + D2M=Sqrt[(4 A B (RM-b) (e-RM))]/(A( RM-b)+B(e- RM)); + + D1P=Sqrt[ -f]; + D2P=Sqrt[(4 A B (RP-b) (e-RP))]/(A(RP- b)+B( e- RP)); + + + AMR[\[Lambda]_]:=JacobiAmplitude[Sqrt[A B J] \[Lambda],kr^2]; + SNR[\[Lambda]_]:=JacobiSN[Sqrt[A B J] \[Lambda],kr^2]; + CNR[\[Lambda]_]:=JacobiCN[Sqrt[A B J] \[Lambda],kr^2]; + DNR[\[Lambda]_]:=JacobiDN[Sqrt[A B J] \[Lambda],kr^2]; + + AMZ[\[Lambda]_]:=JacobiAmplitude[Z2 \[Lambda],kz^2]; + + If[initCoords===Automatic, + {t0,r0,\[Theta]0,\[Phi]0}={0,R1,\[Pi]/2,0}, + {t0,r0,\[Theta]0,\[Phi]0}=initCoords + ]; + If[r0R2, + Message[KerrGeoPlunge::r0outofboundsGen,r0,R2,R1]; + Return[$Failed] + ]; + If[Abs[\[Theta]0]\[Pi] - Abs[ArcCos[Z1]], + Message[KerrGeoPlunge::\[Theta]0outofbounds,\[Theta]0,Abs[ArcCos[Z1]],\[Pi]-Abs[ArcCos[Z1]] ]; + Return[$Failed] + ]; + + + z0 = Cos[\[Theta]0]; + + Minoz[z_] :=InverseJacobiSN[z/Z1,kz^2]/Z2; + If[Abs[Arg]==\[Pi]/2, Minoz[z_] :=0]; + If[a==0, Minoz[z_] :=0]; + If[\[ScriptCapitalQ]==0, Minoz[z_] :=0]; + MinozFunc=Function[{Global`z}, Evaluate[Minoz[Global`z]-Minoz[z0] ],Listable]; + + + MinoR[x_]:= -1/Sqrt[J A B] EllipticF[(\[Pi]/2 -ArcSin[(B(e-x)-A(x-b))/( B(e-x)+A(x-b) ) ]),kr^2]; + MinoRFunc=Function[{Global`r}, Evaluate[MinoR[Global`r] -MinoR[r0] ],Listable]; + + + r[\[Lambda]_] := ((A-B) (A b-B e) SNR[\[Lambda]]^2+2 (A B (b+e)+A B (b-e) CNR[\[Lambda]]))/(4 A B+(A-B)^2 SNR[\[Lambda]]^2); + + z[\[Lambda]_]:= Z1*JacobiSN[Z2 \[Lambda],kz^2]; + If[kz==0, z[\[Lambda]_]:= Z1*Sin[Z2*\[Lambda]]]; + +(*Integrals*) + + RINT\[Lambda][\[Lambda]_] := ((A b-B e)/(A-B) \[Lambda]-1/ Sqrt[ J] ArcTan[(e-b)/(2 Sqrt[A B]) SNR[\[Lambda]]/Sqrt[1-kr^2 (SNR[\[Lambda]])^2]]+((A+B) (e-b))/(2 (A-B) Sqrt[A B J]) EllipticPi[-(1/f),AMR[\[Lambda]],kr^2]); + R2INT\[Lambda][\[Lambda]_]:=\[Lambda] /(A-B) (A b^2-B e^2)+ Sqrt[A B ]/Sqrt[ J] (EllipticE[AMR[\[Lambda]],kr^2])-((A+B) (A^2+2 b^2-B^2-2 e^2))/(4 (A-B) Sqrt[A B J]) EllipticPi[-(1/f),AMR[\[Lambda]],kr^2]-(Sqrt[A B ] (A+B-(A-B)CNR[\[Lambda]]))/((A-B) Sqrt[ J]) (SNR[\[Lambda]] DNR[\[Lambda]])/(f+(SNR[\[Lambda]])^2)+ (A^2+2 b^2-B^2-2 e^2)/(4 (e-b) Sqrt[ J])ArcTan[(f-(1+2 f kr^2) SNR[\[Lambda]]^2),2 SNR[\[Lambda]] DNR[\[Lambda]]Sqrt[f (1+f kr^2)]](*ArcTan[(2 SNR[\[Lambda]] DNR[\[Lambda]]Sqrt[f (1+f kr^2)])/(f-(1+2 f kr^2) SNR[\[Lambda]]^2)]*); + If[a!=0,RMINT\[Lambda][\[Lambda]_] := ((A-B) \[Lambda])/(A (b-RM)-B (e-RM))+((e-b) (A (b-RM)+B (e-RM)))/(2 Sqrt[A B J] (b-RM) (-e+RM) (A (b-RM)-B (e-RM))) EllipticPi[1/D2M^2,JacobiAmplitude[Sqrt[A B J] \[Lambda],kr^2],kr^2]- Sqrt[(e-b)]/(Sqrt[ J] Sqrt[ (RM-b) (e-RM)] Sqrt[ (A^2 (RM-b)-(e-RM) (b^2-B^2+e RM-b (e+RM)))]) 1/4 (Log[((D2M Sqrt[1-D2M^2 kr^2]+Sqrt[1-kr^2 SNR[\[Lambda]]^2]SNR[\[Lambda]])^2+( kr (D2M^2-SNR[\[Lambda]]^2))^2)/((D2M Sqrt[1-D2M^2 kr^2]- Sqrt[1-kr^2 SNR[\[Lambda]]^2]SNR[\[Lambda]])^2+( kr (D2M^2-SNR[\[Lambda]]^2))^2)])]; + If[a==0,RMINT\[Lambda][\[Lambda]_] := ((A-B) \[Lambda])/(A (b-RM)-B (e-RM))]; + RPINT\[Lambda][\[Lambda]_] := ((A-B) \[Lambda])/(A (b-RP)+B (-e+RP))+((e-b) (A (b-RP)+B (e-RP)))/(2 Sqrt[A B J] (b-RP) (-e+RP) (A (b-RP)+B (-e+RP))) EllipticPi[1/D2P^2,JacobiAmplitude[Sqrt[A B J] \[Lambda],kr^2],kr^2]- Sqrt[(e-b)]/(Sqrt[ J] Sqrt[(RP-b) (e-RP)] Sqrt[ (A^2 (RP-b)-(e-RP) (b^2-B^2+e RP-b (e+RP)))]) 1/4 (Log[((D2P Sqrt[1-D2P^2 kr^2]+Sqrt[1-kr^2 SNR[\[Lambda]]^2]SNR[\[Lambda]])^2+( kr (D2P^2-SNR[\[Lambda]]^2))^2)/((D2P Sqrt[1-D2P^2 kr^2]- Sqrt[1-kr^2 SNR[\[Lambda]]^2]SNR[\[Lambda]])^2+( kr (D2P^2-SNR[\[Lambda]]^2))^2)]); + + + tr[\[Lambda]_]:= ((2 a^2 +RM^2+RM RP+RP^2)\[ScriptCapitalE])\[Lambda]+(R2INT\[Lambda][\[Lambda]]+RINT\[Lambda][\[Lambda]](RM+RP)) \[ScriptCapitalE] + ((RM^2+a^2)(\[ScriptCapitalE](RM^2+a^2)-a*\[ScriptCapitalL]))/(RM-RP) RMINT\[Lambda][\[Lambda]]+ ((RP^2+a^2)(\[ScriptCapitalE](RP^2+a^2)-a*\[ScriptCapitalL]))/(RP-RM) RPINT\[Lambda][\[Lambda]]; + \[Phi]r[\[Lambda]_]:= a(((\[ScriptCapitalE](RM^2+a^2)-a*\[ScriptCapitalL])/(RM-RP))RMINT\[Lambda][\[Lambda]]+ (\[ScriptCapitalE](RP^2+a^2)-a*\[ScriptCapitalL])/(RP-RM) RPINT\[Lambda][\[Lambda]]); + tz[\[Lambda]_]:= \[ScriptCapitalE]/ J ((Z2-a^2 J/Z2) Z2 \[Lambda] -Z2 EllipticE[AMZ[\[Lambda]],kz^2]); + \[Phi]z[\[Lambda]_]:= (\[ScriptCapitalL] EllipticPi[Z1^2,AMZ[\[Lambda]],(a^2 J Z1^2)/Z2^2])/Z2; + + + +(*Note: If including Precission Error May occur in the case of evaluating R2INT\[Lambda][0] due to a bug with mathematica in + evaluating terms of the form EllipticE[0, num`n], num is any number and n is less thna machine precision*) + Print[1]; + t=Function[{Global`\[Lambda]}, Evaluate[ tr[Global`\[Lambda]+ MinoR[r0]] + tz[Global`\[Lambda]+ Minoz[z0]]-tr[MinoR[r0]]-tz[Minoz[z0]] + t0], Listable]; + Print[2]; + r=Function[{Global`\[Lambda]}, Evaluate[ r[Global`\[Lambda]+ MinoR[r0]]], Listable]; + Print[3]; + \[Theta]=Function[{Global`\[Lambda]}, Evaluate[ ArcCos[z[Global`\[Lambda] + Minoz[z0]]]] , Listable]; + Print[4]; + \[Phi]=Function[{Global`\[Lambda]}, Evaluate[ \[Phi]r[Global`\[Lambda]+ MinoR[r0]] + \[Phi]z[Global`\[Lambda]+ Minoz[z0]] - \[Phi]r[MinoR[r0]] - \[Phi]z[Minoz[z0]] + \[Phi]0], Listable]; + Print[5]; + MCPP = -Abs[MinoR[RP]] -MinoR[r0]+(2\[Pi])/\[CapitalUpsilon]r; + MCMP = -Abs[MinoR[RM]] -MinoR[r0]+(2\[Pi])/\[CapitalUpsilon]r; + MCMM = Abs[MinoR[RM]] - MinoR[r0]; + MCPM = Abs[MinoR[RP]] - MinoR[r0]; + assoc = Association[ + "a" -> a, + "Parametrization"->"Mino", + "Energy" -> \[ScriptCapitalE], + "AngularMomentum" -> \[ScriptCapitalL], + "CarterConstant" -> \[ScriptCapitalQ], + "ConstantsOfMotion" -> consts, + "RadialFrequency"-> \[CapitalUpsilon]r, + "PolarFrequency"-> \[CapitalUpsilon]\[Theta], + "Frequencies"-> <|"\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(r\)]\)"-> \[CapitalUpsilon]r, "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Theta]\)]\)" -> \[CapitalUpsilon]\[Theta]|>, + "RadialPeriod"-> 2\[Pi]/\[CapitalUpsilon]r, + "PolarPeriod"-> 2\[Pi]/\[CapitalUpsilon]\[Theta], + "Periods"-> <|"\!\(\*SubscriptBox[\(\[Lambda]\), \(r\)]\)"-> 2\[Pi]/\[CapitalUpsilon]r, "\!\(\*SubscriptBox[\(\[Lambda]\), \(\[Theta]\)]\)" -> 2\[Pi]/\[CapitalUpsilon]\[Theta]|>, + "RadialRoots"-> {R1,R2,R3,R4}, + "EllipticBasis" -> <|"\!\(\*SubscriptBox[\(k\), \(r\)]\)"-> kr,"\!\(\*SubscriptBox[\(k\), \(z\)]\)"-> kz|>, + "RadialMinoTime"-> MinoRFunc, + "RadialRootCrossingTimeMino"-> <| + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(Inner\)], \(\)]\)"-> Abs[MinoR[R1]] - MinoR[r0](* Time outgoing to inner root*), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(Outer\)], \(\)]\)"-> Abs[MinoR[R2]] - MinoR[r0](* Time outgoing to outer root*)|>, + "HorizonCrossingTimeMino"-> <| + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(-\)], \(-\)]\)"-> MCMM(* Time outgoing branch crosses inner horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(+\)], \(-\)]\)"-> MCPM(* Time outgoing branch crosses outer horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(+\)], \(+\)]\)"-> MCPP(* Time ingoing branch crosses outer horizon *), + "\!\(\*SubsuperscriptBox[\(\[CapitalLambda]\), SubscriptBox[\(r\), \(-\)], \(+\)]\)"-> MCMP(* Time ingoing branch crosses inner horizon *) + |>, + "PolarRoots"-> {Z1,Z2}, + "Trajectory" -> {t,r,\[Theta],\[Phi]} + ]; + KerrGeoPlungeFunction[a, \[ScriptCapitalE], \[ScriptCapitalL], \[ScriptCapitalQ], assoc] +] + + +(* ::Section:: *) +(*KerrGeoOrbit and KerrGeoOrbitFuction*) + + +KerrGeoPlunge::fourcomplexroots = "This function can only currently solve plunging geodesics in the paramater space where we have four finite real roots or two finite real roots and two finite complex roots."; +KerrGeoPlunge::highenergy = "This code does not currently support plunging solutions with \[ScriptCapitalE]>1, \[ScriptCapitalE]=`1` has been selected "; +KerrGeoPlunge::negativecarter = "This code does not currently support plunging solutions with Q<0, Q=`1` has been selected "; +KerrGeoPlunge::rIoutbounds = "Given ISSO radius `1` is not in the allowed range between `2` and `3` for spin `4`."; +KerrGeoPlunge::Inclinationoutofbounds = "Given inclination angle `1` is not between -\!\(\*FractionBox[\(\[Pi]\), \(2\)]\) and \!\(\*FractionBox[\(\[Pi]\), \(2\)]\)." +KerrGeoPlunge::aoutofbounds1 = "The plunges packages does not currently support the extremal a=1 case" +KerrGeoPlunge::aoutofbounds0 = "The plunges packages does not currently support the zero spin a=0 case" + + +KerrGeoPlunge[a_:0.9,{\[ScriptCapitalE]_:0.8,\[ScriptCapitalL]_:0.3,\[ScriptCapitalQ]_:3}, initCoords:{_,_,_,_}:Automatic,OptionsPattern[]]:=Module[ + {param, method, RI, DN, RISSOMIN, RISSOMAX, ROOTS, RealRoots,ComplexRoots}, + + + If[a==1, + Message[KerrGeoPlunge::aoutofbounds1]; + Return[$Failed] + ]; + + If[\[ScriptCapitalE]>=1, + Message[KerrGeoPlunge::highenergy,\[ScriptCapitalE]]; + Return[$Failed] + ]; + + If[\[ScriptCapitalQ]<0, + Message[KerrGeoPlunge::negativecarter,\[ScriptCapitalQ]]; + Return[$Failed] + ]; + + ROOTS = r/.NSolve[(\[ScriptCapitalE](r^2+a^2)-a*\[ScriptCapitalL])^2-(r^2-2r+a^2)(r^2+(a*\[ScriptCapitalE]-\[ScriptCapitalL])^2+\[ScriptCapitalQ])==0,r]; + + RealRoots = Select[ROOTS,PossibleZeroQ@Im[#]&]; + ComplexRoots = Select[ROOTS,Not@PossibleZeroQ@Im[#]&]; + + If[Length[ComplexRoots]== 2, Return[KerrGeoComplexPlungeMino[a, \[ScriptCapitalE], \[ScriptCapitalL], \[ScriptCapitalQ], initCoords,"Roots"->ROOTS]]]; + If[Length[ComplexRoots]== 0, Return[KerrGeoRealPlungeMino[a, \[ScriptCapitalE], \[ScriptCapitalL], \[ScriptCapitalQ], initCoords,"Roots"->ROOTS]]]; + ]; + + +KerrGeoPlunge[a_:0.9,"ISSORadialParam",Arg1_:0.8, initCoords:{_,_,_,_}:Automatic,OptionsPattern[]]:= + Module[{RI,DN,RISSOMIN,RISSOMAX}, + If[a==1, + Message[KerrGeoPlunge::aoutofbounds1,a]; + Return[$Failed] + ]; + RI = Arg1; + DN = (27-45 a^2+17 a^4+a^6+8 a^3 (1-a^2))^(1/3); + RISSOMIN = 3+Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]-1/2 \[Sqrt](72+8 (-6+a^2)-(4 (9-10 a^2+a^4))/DN-4 DN+(64 a^2)/Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]); + RISSOMAX = 3+Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]+1/2 \[Sqrt](72+8 (-6+a^2)-(4 (9-10 a^2+a^4))/DN-4 DN+(64 a^2)/Sqrt[3+a^2+(9-10 a^2+a^4)/DN+DN]); + If[Not[RISSOMIN<=RI<=RISSOMAX], + Message[KerrGeoPlunge::rIoutbounds,RI,RISSOMIN,RISSOMAX,a]; + Return[$Failed] + ]; + KerrGeoISSOPlunge[a,"ISSORadialParam" ,RI, initCoords] + ] + + +KerrGeoPlunge[a_:0.9,"ISSOIncParam",Arg1_:0.8, initCoords:{_,_,_,_}:Automatic,OptionsPattern[]]:= + Module[{}, + If[a==1, + Message[KerrGeoPlunge::aoutofbounds1,a]; + Return[$Failed] + ]; + (* Other inclinations should be fine as well, right?*) + If[Not[-\[Pi]/2<=Arg1<=\[Pi]/2], + Message[KerrGeoPlunge::Inclinationoutofbounds, Arg1]; + Return[$Failed] + ]; + KerrGeoISSOPlunge[a, "ISSOIncParam" ,Arg1 ,initCoords] + ] + + +KerrGeoPlungeFunction /: + MakeBoxes[kgof:KerrGeoPlungeFunction[a_,\[ScriptCapitalE]_,\[ScriptCapitalL]_,\[ScriptCapitalQ]_,assoc_], form:(StandardForm|TraditionalForm)] := + Module[{summary, extended}, + summary = {Row[{BoxForm`SummaryItem[{"a: ", a}], " ", + BoxForm`SummaryItem[{"\[ScriptCapitalE]: ", \[ScriptCapitalE]}], " ", + BoxForm`SummaryItem[{"\[ScriptCapitalL]: ", \[ScriptCapitalL]}], " ", + BoxForm`SummaryItem[{"\[ScriptCapitalQ]: ", \[ScriptCapitalQ]}]}], + BoxForm`SummaryItem[{"Parametrization: ", assoc["Parametrization"]}]}; + extended = {BoxForm`SummaryItem[{"Energy: ", assoc["Energy"]}], + BoxForm`SummaryItem[{"Angular Momentum: ", assoc["AngularMomentum"]}], + BoxForm`SummaryItem[{"Carter Constant: ", assoc["CarterConstant"]}]}; + BoxForm`ArrangeSummaryBox[ + KerrGeoPlungeFunction, + kgof, + None, + summary, + extended, + form] +]; + + +KerrGeoPlungeFunction[a_, \[ScriptCapitalE]_,\[ScriptCapitalL]_,\[ScriptCapitalQ]_,assoc_][\[Lambda]_/;StringQ[\[Lambda]] == False] := Through[assoc["Trajectory"][\[Lambda]]] +KerrGeoPlungeFunction[a_,\[ScriptCapitalE]_,\[ScriptCapitalL]_,\[ScriptCapitalQ]_,assoc_][y_?StringQ] := assoc[y] +Keys[g_KerrGeoPlungeFunction]^:=Keys[g[[5]]] + + +(* ::Section:: *) +(*Close the package*) + + +End[]; + +EndPackage[]; diff --git a/Kernel/KerrGeodesics.m b/Kernel/KerrGeodesics.m index f214e36..55bfd86 100644 --- a/Kernel/KerrGeodesics.m +++ b/Kernel/KerrGeodesics.m @@ -5,6 +5,9 @@ "KerrGeodesics`OrbitalFrequencies`", "KerrGeodesics`SpecialOrbits`", "KerrGeodesics`KerrGeoOrbit`", + "KerrGeodesics`KerrGeoPlunge`", "KerrGeodesics`FourVelocity`"}] - EndPackage[] + + + diff --git a/Source/Documentation/English/Guides/KerrGeodesics.md b/Source/Documentation/English/Guides/KerrGeodesics.md index e86b32c..850cab9 100644 --- a/Source/Documentation/English/Guides/KerrGeodesics.md +++ b/Source/Documentation/English/Guides/KerrGeodesics.md @@ -25,6 +25,7 @@ {"KerrGeoISSO", "Computes the location of the inner-most stable spherical orbit (ISSO)"}, {"KerrGeoSeparatrix", "Computes the value p at the separatrix between stable and plunging/scattered orbits"}, {"KerrGeoOrbit", "Computes the orbit trajectory in Boyer-Lindquist coordinates"}, + {"KerrGeoPlunges", "Computes the Plunging trajectory in Boyer-Lindquist coordinates"}, {"KerrGeoFourVelocity", "Computes the components of the body's four-velocity"}, {"KerrGeoFindResonances", "Computes the location of orbital resonances"} } diff --git a/Tests/KerrGeoPlunge.wlt b/Tests/KerrGeoPlunge.wlt new file mode 100644 index 0000000..2f2933d --- /dev/null +++ b/Tests/KerrGeoPlunge.wlt @@ -0,0 +1,228 @@ +BeginTestSection["PlungesTest"] + +VerificationTest[(* 1 *) + Block[{t, r, \[Theta], \[Phi]}, orbit=KerrGeoPlunge[0.05, "ISSORadialParam", 5.9 ]; +{t, r, \[Theta], \[Phi]} =orbit["Trajectory"]; +{t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]}/.\[Lambda]->5] + , + {166.93605140000858, 5.669341902321057, 2.4860551451282125, 17.222045750314397} + , + TestID -> "KerrGeoPlunge - ISSORadialParam" +] + +VerificationTest[(* 2 *) + orbit=KerrGeoPlunge[0.05, "ISSORadialParam", 200.9 ] + , + $Failed + , + {KerrGeoPlunge::rIoutbounds} + , + TestID -> "KerrGeoPlunge -ISSORadialParam - ROutOfBoundsUpper" +] + +VerificationTest[(* 3 *) + orbit=KerrGeoPlunge[0.05, "ISSORadialParam", 0 ] + , + $Failed + , + {KerrGeoPlunge::rIoutbounds} + , + TestID -> "KerrGeoPlunge -ISSORadialParam - ROutOfBoundsUnder" +] + +VerificationTest[(* 4 *) + Block[{t, r, \[Theta], \[Phi]}, orbit=KerrGeoPlunge[0.08, "ISSORadialParam", 6.06, {63 , 1.1,2.4, -140}]; +{t, r, \[Theta], \[Phi]} =orbit["Trajectory"]; +{t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]}/.\[Lambda]->8.9] + , + {456.9800421782238, 5.99235372581383, 2.771355583007127, -170.597372084777} + , + TestID -> "KerrGeoPlunge - ISSORadialParamICs" +] + +VerificationTest[(* 5 *) + KerrGeoPlunge[0.08, "ISSORadialParam", 6.06, {63 , 78.1,2.4, -140}] + , + $Failed + , + {KerrGeoPlunge::r0outofbounds} + , + TestID -> "KerrGeoPlunge - ISSORadialParamICs - RadialOutOfBounds" +] + +VerificationTest[(* 6 *) + KerrGeoPlunge[0.08, "ISSORadialParam", 6.06, {63 , 2.1,12.4, -140}] + , + $Failed + , + {KerrGeoPlunge::\[Theta]0outofbounds} + , + TestID -> "KerrGeoPlunge - ISSORadialParamICs - PolarOutOfBounds" +] + +VerificationTest[(* 7 *) + Block[{t, r, \[Theta], \[Phi]}, orbit=KerrGeoPlunge[0.08, "ISSORadialParam", 6.06, {63` ,(11/10),(12/5), -140`}]; +{t, r, \[Theta], \[Phi]} =orbit["Trajectory"]; +Round[{t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]}/.\[Lambda]->0, 10^(-8)]] + , + {63, 11/10, 12/5, -140} + , + TestID -> "KerrGeoPlunge - ISSORadialParamICs - TimeZero" +] + +VerificationTest[(* 8 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.99, "ISSORadialParam", 8.717352279606489]; + {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> 8.99] + , + {709.2968761641665, 8.643214446769269, 1.6120407986455219, -34.4819996184413} + , + TestID -> "KerrGeoPlunge - ISSORadialParam MaximalRI" +] + +VerificationTest[(* 9 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.7, "ISSOIncParam", Pi/2.5]; {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; + {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> -7.8] + , + {-136.21547228952554, 3.389631843209988, 1.8758140096570342, -24.6802217009732} + , + TestID -> "KerrGeoPlunge - ISSOIncParam" +] + +VerificationTest[(* 10 *) + KerrGeoPlunge[0.7, "ISSOIncParam", 15] + , + $Failed + , + {KerrGeoPlunge::Inclinationoutofbounds} + , + TestID -> "KerrGeoPlunge - ISSOIncParam MaxIncOutOfBounds" +] + +VerificationTest[(* 11 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.77, "ISSOIncParam", -(Pi/4.2), {123, 4.3, Pi/2.2, 12.8}]; + {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> -2.09] + , + {95.35960150614794, 4.282702895995593, 2.391128629110148, 21.014943131587707} + , + TestID -> "KerrGeoPlunge - ISSOIncParam MaxIncICs" +] + +VerificationTest[(* 12 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.77, "ISSOIncParam", -(Pi/2), {123, 4.3, Pi/2, 12.8}]; + {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> -2.09] + , + {87.79211403209042, 5.614689072605751, 1.5707963267948966, 21.38125127100541} + , + TestID -> "KerrGeoPlunge - ISSOIncParamICs - ExtremeInc Negative" +] + +VerificationTest[(* 13 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.977, "ISSOIncParam", Pi/2, {123, 1.3, Pi/2, 12.8}]; + {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> -2.09] + , + {105.04498996779, 0.6649464562686522, 1.5707963267948966, 4.598376768609648} + , + TestID -> "KerrGeoPlunge - ISSOIncParamICs - ExtremeInc Positive" +] + +VerificationTest[(* 14 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.92, {0.8, 10, 6}]; {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; + {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> 3.66] + , + {20.661343459257957, 0.15028185959212795, 1.5630908171743572, 38.368616042198404} + , + TestID -> "KerrGeoPlunge - Generic Plunge - Complex Roots" +] + +VerificationTest[(* 15 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.987, {0.101, 0.39, 0}]; {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; + {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> 3.66] + , + {2.4976443942888236, 1.0611065184188126, 1.5707963267948966, -0.259629783302759} + , + {Power::infy, Power::infy, Infinity::indet, Power::infy, General::stop, Infinity::indet, Infinity::indet, General::stop} + , + TestID -> "KerrGeoPlunge - Generic Plunge - Real Roots 3 Inside rm" +] + +VerificationTest[(* 16 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.1, {0.956, 2.55, 6.54}]; {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; + {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> 3.66] + , + {48.654713642077695, 0.8239191438889363, 1.126915405459793, 13.12077316222189} + , + TestID -> "KerrGeoPlunge - Generic Plunge - Real Roots 3 Outside rp" +] + +VerificationTest[(* 17 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.92, {0.8, 10, 0}]; {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; + {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> 3.66] + , + {18.729898545961788, 0.9316108659154707, 1.5707963267948966, 37.204816126657136} + , + {Power::infy, Power::infy, Infinity::indet} + , + TestID -> " KerrGeoPlunge - Generic - Zero Carter" +] + +VerificationTest[(* 18 *) + KerrGeoPlunge[0.92, {10, 10, 2}] + , + $Failed + , + {KerrGeoPlunge::highenergy} + , + TestID -> "KerrGeoPlunge - Generic - High Energy" +] + +VerificationTest[(* 19 *) + KerrGeoPlunge[0.92, {0.6, 2.4, -5}] + , + $Failed + , + {KerrGeoPlunge::negativecarter} + , + TestID -> "KerrGeoPlunge - Generic - Negative Carter" +] + +VerificationTest[(* 20 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.9, {0.3, 4, 10}, {-221, 1, Pi/2.9, 42}]; {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; + {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> 8.66] + , + {-186.53361138627804, 1.430951644567102, 0.9346228824493619, 91.0672074418539} + , + TestID -> "KerrGeoPlunge - Generic - ICs" +] + +VerificationTest[(* 21 *) + Block[{t, r, \[Theta], \[Phi]}, orbit = KerrGeoPlunge[0.987, {0.101, 0.39, 0}]; {t, r, \[Theta], \[Phi]} = orbit["Trajectory"]; + {t[\[Lambda]], r[\[Lambda]], \[Theta][\[Lambda]], \[Phi][\[Lambda]]} /. \[Lambda] -> 3.66] + , + {2.4976443942888236, 1.0611065184188126, 1.5707963267948966, -0.259629783302759} + , + {Power::infy, Power::infy, Infinity::indet, Power::infy, General::stop, Infinity::indet, Infinity::indet, General::stop} + , + TestID -> "KerrGeoPlunge - RealRoot QZero" +] + +VerificationTest[(* 22 *) + KerrGeoPlunge[0.9, {0.3, 4, 10}, {-221, 1101, Pi/2.9, 42}] + , + $Failed + , + {KerrGeoPlunge::r0outofboundsGen} + , + TestID -> "KerrGeoPlunge - Generic - IC radius out of bounds" +] + +VerificationTest[(* 23 *) + KerrGeoPlunge[0.9, {0.3, 4, 10}, {-221, 1, Pi/0.9, 42}] + , + $Failed + , + {KerrGeoPlunge::\[Theta]0outofbounds} + , + TestID -> "KerrGeoPlunge - Generic - IC polar out of bounds" +] + +EndTestSection[]