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

Parallel transport dev #32

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
+190 −91
Open

Parallel transport dev #32

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
aab5522
Expose MarckCarterFrame and PrecessionPhase
Jul 28, 2022
da58399
First pass at adding phase parametrization.
Jul 28, 2022
0f0e95d
Expose Delta's, implement "Phases" parameterization, various polishments
Jul 28, 2022
149e435
Second pass at implementing "Phases" parametrization
Jul 28, 2022
2bc026b
Fix a bug causing unit tests to fail for KerrGeoOrbit
Jul 29, 2022
d902e49
1) implement `Keys` for ParallelTransportFrameFunction
Jul 29, 2022
1657a5c
Added InitialPhases key to association for KerrGeoOrbit and have exte…
znasipak Jul 29, 2022
b7c1059
Merge branch 'ParallelTransport_Dev' of https://github.com/BlackHoleP…
znasipak Jul 29, 2022
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
114 changes: 73 additions & 41 deletions Kernel/KerrGeoOrbit.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@
]


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


Expand Down Expand Up @@ -158,6 +158,8 @@

]



(* ::Subsection::Closed:: *)
(*Equatorial (Fast Spec - Darwin)*)

Expand Down Expand Up @@ -404,7 +406,6 @@
(*FIXME: make the initial phases work in this case*)



KerrGeoOrbitMino[a_, p_, (0|0.), (1|1.), initPhases:{_,_,_,_}:{0,0,0,0}] := Module[{consts, assoc, t, r, \[Theta], \[Phi], En, L, Q, \[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]t, e=0, x=1,r1,r2,r3,r4,velocity},

{\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]t} = Values[KerrGeodesics`OrbitalFrequencies`Private`KerrGeoMinoFrequencies[a,p,e,x]];
Expand Down Expand Up @@ -445,12 +446,11 @@



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

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


KerrGeoOrbitMino[a_,p_,e_,x_,initPhases:{_,_,_,_}:{0,0,0,0}]:=Module[{M=1,consts,En,L,Q,assoc,\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]t,r1,r2,r3,r4,zp,zm,kr,k\[Theta],rp,rm,hr,hp,hm,rq,zq,\[Psi]r,tr,\[Phi]f,\[Psi]z,tz,\[Phi]z,qt0,qr0,qz0,q\[Phi]0,t,r,\[Theta],\[Phi],\[Phi]t,\[Phi]r,Ct,C\[Phi],velocity},
KerrGeoOrbitPhases[a_,p_,e_,x_]:=Module[{M=1,consts,En,L,Q,assoc,\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]t,r1,r2,r3,r4,zp,zm,kr,k\[Theta],rp,rm,hr,hp,hm,rq,zq,\[Psi]r,tr,\[Phi]f,\[Psi]z,tz,\[Phi]z,t,r,\[Theta],\[Phi],\[Phi]t,\[Phi]r,Ct,C\[Phi],velocity,\[CapitalDelta]tr,\[CapitalDelta]t\[Theta],\[CapitalDelta]\[Phi]r,\[CapitalDelta]\[Phi]\[Theta]},
consts = KerrGeoConstantsOfMotion[a,p,e,x];
{En,L,Q} = Values[consts];
{\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]t} = Values[KerrGeodesics`OrbitalFrequencies`Private`KerrGeoMinoFrequencies[a,p,e,x]];
Expand All @@ -460,62 +460,60 @@
kr = (r1-r2)/(r1-r3) (r3-r4)/(r2-r4);
k\[Theta] = a^2 (1-En^2)(zm/zp)^2;

rp=M+Sqrt[M^2-a^2];
rm=M-Sqrt[M^2-a^2];
hr=(r1-r2)/(r1-r3);
hp=((r1-r2)(r3-rp))/((r1-r3)(r2-rp));
hm=((r1-r2)(r3-rm))/((r1-r3)(r2-rm));
rp=M+Sqrt[M^2-a^2];
rm=M-Sqrt[M^2-a^2];
hr=(r1-r2)/(r1-r3);
hp=((r1-r2)(r3-rp))/((r1-r3)(r2-rp));
hm=((r1-r2)(r3-rm))/((r1-r3)(r2-rm));

rq = Function[{qr},(r3(r1 - r2)JacobiSN[EllipticK[kr]/\[Pi] qr,kr]^2-r2(r1-r3))/((r1-r2)JacobiSN[EllipticK[kr]/\[Pi] qr,kr]^2-(r1-r3))];
rq = Function[{qr},(r3(r1 - r2)JacobiSN[EllipticK[kr]/\[Pi] qr,kr]^2-r2(r1-r3))/((r1-r2)JacobiSN[EllipticK[kr]/\[Pi] qr,kr]^2-(r1-r3))];

zq = Function[{qz}, zm JacobiSN[EllipticK[k\[Theta]] 2/\[Pi] (qz+\[Pi]/2),k\[Theta]]];
zq = Function[{qz}, zm JacobiSN[EllipticK[k\[Theta]] 2/\[Pi] (qz+\[Pi]/2),k\[Theta]]];

\[Psi]r[qr_]:=\[Psi]r[qr]= JacobiAmplitude[EllipticK[kr]/\[Pi] qr,kr];
\[Psi]r[qr_]:=\[Psi]r[qr]= JacobiAmplitude[EllipticK[kr]/\[Pi] qr,kr];

tr[qr_]:= -En/Sqrt[(1-En^2) (r1-r3) (r2-r4)] (
4(r2-r3) (EllipticPi[hr,kr] qr/\[Pi]-EllipticPi[hr,\[Psi]r[qr],kr])
-4 (r2-r3)/(rp-rm) (
-(1/((-rm+r2) (-rm+r3)))(-2 a^2+rm (4-(a L)/En)) (EllipticPi[hm,kr] qr/\[Pi]-EllipticPi[hm,\[Psi]r[qr],kr] )
+1/((-rp+r2) (-rp+r3)) (-2 a^2+rp (4-(a L)/En)) (EllipticPi[hp,kr] qr/\[Pi]-EllipticPi[hp,\[Psi]r[qr],kr])
)
+(r2-r3) (r1+r2+r3+r4) (EllipticPi[hr,kr] qr/\[Pi]-EllipticPi[hr,\[Psi]r[qr],kr] )
+(r1-r3) (r2-r4) (EllipticE[kr] qr/\[Pi]-EllipticE[\[Psi]r[qr],kr]+hr((Sin[\[Psi]r[qr]]Cos[\[Psi]r[qr]] Sqrt[1-kr Sin[\[Psi]r[qr]]^2])/(1-hr Sin[\[Psi]r[qr]]^2))) );
tr[qr_]:= -En/Sqrt[(1-En^2) (r1-r3) (r2-r4)] (
4(r2-r3) (EllipticPi[hr,kr] qr/\[Pi]-EllipticPi[hr,\[Psi]r[qr],kr])
-4 (r2-r3)/(rp-rm) (
-(1/((-rm+r2) (-rm+r3)))(-2 a^2+rm (4-(a L)/En)) (EllipticPi[hm,kr] qr/\[Pi]-EllipticPi[hm,\[Psi]r[qr],kr] )
+1/((-rp+r2) (-rp+r3)) (-2 a^2+rp (4-(a L)/En)) (EllipticPi[hp,kr] qr/\[Pi]-EllipticPi[hp,\[Psi]r[qr],kr])
)
+(r2-r3) (r1+r2+r3+r4) (EllipticPi[hr,kr] qr/\[Pi]-EllipticPi[hr,\[Psi]r[qr],kr] )
+(r1-r3) (r2-r4) (EllipticE[kr] qr/\[Pi]-EllipticE[\[Psi]r[qr],kr]+hr((Sin[\[Psi]r[qr]]Cos[\[Psi]r[qr]] Sqrt[1-kr Sin[\[Psi]r[qr]]^2])/(1-hr Sin[\[Psi]r[qr]]^2))) );

\[Phi]r[qr_]:= (2 a En (-1/((-rm+r2) (-rm+r3))(2 rm-(a L)/En) (r2-r3) (EllipticPi[hm,kr] qr/\[Pi]-EllipticPi[hm,\[Psi]r[qr],kr])+1/((-rp+r2) (-rp+r3))(2 rp-(a L)/En) (r2-r3) (EllipticPi[hp,kr] qr/\[Pi]-EllipticPi[hp,\[Psi]r[qr],kr] )))/((-rm+rp) Sqrt[(1-En^2) (r1-r3) (r2-r4)]);
\[Phi]r[qr_]:= (2 a En (-1/((-rm+r2) (-rm+r3))(2 rm-(a L)/En) (r2-r3) (EllipticPi[hm,kr] qr/\[Pi]-EllipticPi[hm,\[Psi]r[qr],kr])+1/((-rp+r2) (-rp+r3))(2 rp-(a L)/En) (r2-r3) (EllipticPi[hp,kr] qr/\[Pi]-EllipticPi[hp,\[Psi]r[qr],kr] )))/((-rm+rp) Sqrt[(1-En^2) (r1-r3) (r2-r4)]);

\[Psi]z[qz_]:= \[Psi]z[qz] = JacobiAmplitude[EllipticK[k\[Theta]] 2/\[Pi] (qz+\[Pi]/2),k\[Theta]];
tz[qz_]:= 1/(1-En^2) En zp ( EllipticE[k\[Theta]]2((qz+\[Pi]/2)/\[Pi])-EllipticE[\[Psi]z[qz],k\[Theta]]);
\[Phi]z[qz_]:= -1/zp L ( EllipticPi[zm^2,k\[Theta]]2((qz+\[Pi]/2)/\[Pi])-EllipticPi[zm^2,\[Psi]z[qz],k\[Theta]]);
\[Psi]z[qz_]:= \[Psi]z[qz] = JacobiAmplitude[EllipticK[k\[Theta]] 2/\[Pi] (qz+\[Pi]/2),k\[Theta]];
tz[qz_]:= 1/(1-En^2) En zp ( EllipticE[k\[Theta]]2((qz+\[Pi]/2)/\[Pi])-EllipticE[\[Psi]z[qz],k\[Theta]]);
\[Phi]z[qz_]:= -1/zp L ( EllipticPi[zm^2,k\[Theta]]2((qz+\[Pi]/2)/\[Pi])-EllipticPi[zm^2,\[Psi]z[qz],k\[Theta]]);

{qt0, qr0, qz0, q\[Phi]0} = {initPhases[[1]], initPhases[[2]], initPhases[[3]], initPhases[[4]]};
(*Calculate normalization constants so that t=0 and \[Phi]=0 at \[Lambda]=0 when qt0=0 and q\[Phi]0=0 *)
Ct=tr[qr0]+tz[qz0]/.i_/;i==0:>0;
C\[Phi]=\[Phi]r[qr0]+\[Phi]z[qz0]/.i_/;i==0:>0;

t=Function[{Global`\[Lambda]}, Evaluate[ qt0 + \[CapitalUpsilon]t Global`\[Lambda] + tr[\[CapitalUpsilon]r Global`\[Lambda] + qr0] + tz[\[CapitalUpsilon]\[Theta] Global`\[Lambda] + qz0]-Ct ], Listable];
r=Function[{Global`\[Lambda]}, Evaluate[ rq[\[CapitalUpsilon]r Global`\[Lambda]+ qr0] ], Listable];
\[Theta]=Function[{Global`\[Lambda]}, Evaluate[ ArcCos[zq[\[CapitalUpsilon]\[Theta] Global`\[Lambda] + qz0]] ], Listable];
\[Phi]=Function[{Global`\[Lambda]}, Evaluate[ q\[Phi]0 + \[CapitalUpsilon]\[Phi] Global`\[Lambda] + \[Phi]r[\[CapitalUpsilon]r Global`\[Lambda]+ qr0] + \[Phi]z[\[CapitalUpsilon]\[Theta] Global`\[Lambda] + qz0]-C\[Phi] ], Listable];

velocity = Values[KerrGeoFourVelocity[a,p,e,x,{initPhases[[2]],initPhases[[3]]}]];
t=Function[{Global`qt,Global`qr,Global`q\[Theta]}, Evaluate[ Global`qt+ tr[Global`qr] + tz[Global`q\[Theta]]], Listable];
r=Function[{Global`qr}, Evaluate[ rq[Global`qr] ], Listable];
\[Theta]=Function[{Global`q\[Theta]}, Evaluate[ ArcCos[zq[Global`q\[Theta]]] ], Listable];
\[Phi]=Function[{Global`q\[Phi],Global`qr,Global`q\[Theta]}, Evaluate[ Global`q\[Phi] + \[Phi]r[Global`qr] + \[Phi]z[Global`q\[Theta]]], Listable];

\[CapitalDelta]tr=Function[{Global`qr}, Evaluate[ tr[Global`qr] ], Listable];
\[CapitalDelta]t\[Theta]=Function[{Global`q\[Theta]}, Evaluate[ tz[Global`q\[Theta]] ], Listable];
\[CapitalDelta]\[Phi]r=Function[{Global`qr}, Evaluate[ \[Phi]r[Global`qr] ], Listable];
\[CapitalDelta]\[Phi]\[Theta]=Function[{Global`q\[Theta]}, Evaluate[ \[Phi]z[Global`q\[Theta]] ], Listable];

assoc = Association[
"a" -> a,
"p" -> p,
"e" -> e,
"Inclination" -> x,
"Parametrization"->"Mino",
"Parametrization"->"Phases",
"Energy" -> En,
"AngularMomentum" -> L,
"CarterConstant" -> Q,
"ConstantsOfMotion" -> consts,
"RadialFrequency" -> \[CapitalUpsilon]r,
"PolarFrequency" -> \[CapitalUpsilon]\[Theta],
"AzimuthalFrequency" -> \[CapitalUpsilon]\[Phi],
Copy link
Collaborator

Choose a reason for hiding this comment

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

Should we add a "TimeFrequency" key as well?

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 have added it to the "Frequencies", but did not added a top level key, because I did not know what to best call it.

"Frequencies" -> <|"\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(r\)]\)" -> \[CapitalUpsilon]r, "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Theta]\)]\)" -> \[CapitalUpsilon]\[Theta], "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Phi]\)]\)" -> \[CapitalUpsilon]\[Phi] |> ,
"Frequencies" -> <|"\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(r\)]\)" -> \[CapitalUpsilon]r, "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Theta]\)]\)" -> \[CapitalUpsilon]\[Theta], "\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Phi]\)]\)" -> \[CapitalUpsilon]\[Phi] ,"\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(t\)]\)" -> \[CapitalUpsilon]t |> ,
"Trajectory" -> {t,r,\[Theta],\[Phi]},
"FourVelocity"-> velocity,
"TrajectoryDeltas" -> <|"\[CapitalDelta]tr"-> \[CapitalDelta]tr,"\[CapitalDelta]t\[Theta]"-> \[CapitalDelta]t\[Theta],"\[CapitalDelta]\[Phi]r"-> \[CapitalDelta]\[Phi]r,"\[CapitalDelta]\[Phi]\[Theta]"-> \[CapitalDelta]\[Phi]\[Theta]|>,
"RadialRoots" -> {r1,r2,r3,r4},
"a" -> a,
"p" -> p,
Expand All @@ -528,6 +526,39 @@
]


KerrGeoOrbitMino[a_,p_,e_,x_,initPhases:{_,_,_,_}:{0,0,0,0}]:=Module[{M=1,assoc,qt0,qr0,qz0,q\[Phi]0,tr,tz,t,r,\[Theta],\[Phi],\[Phi]z,\[Phi]r,Ct,C\[Phi],velocity,\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]t},
assoc= Last@KerrGeoOrbitPhases[a,p,e,x];

{qt0, qr0, qz0, q\[Phi]0} = {initPhases[[1]], initPhases[[2]], initPhases[[3]], initPhases[[4]]};

{\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta],\[CapitalUpsilon]\[Phi],\[CapitalUpsilon]t}=Values@assoc["Frequencies"];
{t,r,\[Theta],\[Phi]}=assoc["Trajectory"];

tr=assoc["TrajectoryDeltas"]["\[CapitalDelta]tr"];
tz=assoc["TrajectoryDeltas"]["\[CapitalDelta]t\[Theta]"];
\[Phi]r=assoc["TrajectoryDeltas"]["\[CapitalDelta]\[Phi]r"];
\[Phi]z=assoc["TrajectoryDeltas"]["\[CapitalDelta]\[Phi]\[Theta]"];

(*Calculate normalization constants so that t=0 and \[Phi]=0 at \[Lambda]=0 when qt0=0 and q\[Phi]0=0 *)
Ct=tr[qr0]+tz[qz0]/.i_/;i==0:>0;
C\[Phi]=\[Phi]r[qr0]+\[Phi]z[qz0]/.i_/;i==0:>0;

t=Function[{Global`\[Lambda]}, Evaluate[ t[qt0 + \[CapitalUpsilon]t Global`\[Lambda], \[CapitalUpsilon]r Global`\[Lambda] + qr0, \[CapitalUpsilon]\[Theta] Global`\[Lambda] + qz0]-Ct ], Listable];
r=Function[{Global`\[Lambda]}, Evaluate[ r[\[CapitalUpsilon]r Global`\[Lambda]+ qr0] ], Listable];
\[Theta]=Function[{Global`\[Lambda]}, Evaluate[ \[Theta][\[CapitalUpsilon]\[Theta] Global`\[Lambda] + qz0] ], Listable];
\[Phi]=Function[{Global`\[Lambda]}, Evaluate[ \[Phi][q\[Phi]0 + \[CapitalUpsilon]\[Phi] Global`\[Lambda] , \[CapitalUpsilon]r Global`\[Lambda]+ qr0, \[CapitalUpsilon]\[Theta] Global`\[Lambda] + qz0]-C\[Phi] ], Listable];

velocity = Values[KerrGeoFourVelocity[a,p,e,x,{qr0,qz0}]];

assoc["Parametrization"] = "Mino";
assoc["Trajectory"] = {t,r,\[Theta],\[Phi]};
assoc["FourVelocity"] = velocity;

KerrGeoOrbitFunction[a,p,e,x,assoc]

]


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

Expand Down Expand Up @@ -831,7 +862,7 @@
];


(* ::Subsubsection:: *)
(* ::Subsubsection::Closed:: *)
(*Subroutine that performs spectral integration on even functions*)


Expand Down Expand Up @@ -868,7 +899,7 @@

(* Create functions for discrete cosine series coefficient fn *)
samplePhase=(freq \[Lambda]+phaseList);
fn[n_]:=fn[n]=(sampledF.Cos[nn*samplePhase])/.nn->n;
fn[n_]:=fn[n]=(sampledF . Cos[nn*samplePhase])/.nn->n;

(* Calculate series coefficients until they equal 0 (with respect
to the precision being used) *)
Expand Down Expand Up @@ -897,7 +928,7 @@
];


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


Expand Down Expand Up @@ -1062,6 +1093,7 @@
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]]]
];
Expand Down
100 changes: 77 additions & 23 deletions Kernel/ParallelTransport.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,13 +156,13 @@


PrecessionPhase[orbit_KerrGeoOrbitFunction,\[CapitalUpsilon]\[Psi]_,qr0_,qz0_,q\[Psi]0_:0]:=
Function[{\[Lambda]},
Function[{Global`\[Lambda]},
Evaluate@With[
{
\[CapitalUpsilon]r=orbit["RadialFrequency"],
\[CapitalUpsilon]z=orbit["PolarFrequency"]
},
q\[Psi]0+\[CapitalUpsilon]\[Psi] \[Lambda] +PrecessionPhaser[orbit][\[CapitalUpsilon]r \[Lambda] + qr0]+PrecessionPhasez[orbit][\[CapitalUpsilon]z \[Lambda] + qz0]
q\[Psi]0+\[CapitalUpsilon]\[Psi] Global`\[Lambda] +PrecessionPhaser[orbit][\[CapitalUpsilon]r Global`\[Lambda] + qr0]+PrecessionPhasez[orbit][\[CapitalUpsilon]z Global`\[Lambda] + qz0]
]
]

Expand All @@ -174,21 +174,9 @@
(*Tetrads*)


MarckCarterFrame[orbit_KerrGeoOrbitFunction]:=
Module[
{
a,\[ScriptCapitalE],\[ScriptCapitalL],\[ScriptCapitalQ],\[ScriptCapitalK],rf,\[Theta]f
},
a=orbit["a"];
{\[ScriptCapitalE],\[ScriptCapitalL],\[ScriptCapitalQ]}=Values@orbit["ConstantsOfMotion"];
\[ScriptCapitalK]=\[ScriptCapitalQ]+(a \[ScriptCapitalE]-\[ScriptCapitalL])^2;
rf=Function[\[Lambda],orbit["Trajectory"][[2]][\[Lambda]]];
\[Theta]f=Function[\[Lambda],orbit["Trajectory"][[3]][\[Lambda]]];
Function[{\[Lambda]},
Evaluate@With[{
r=rf[\[Lambda]],\[Theta]=\[Theta]f[\[Lambda]],rd=rf'[\[Lambda]],\[Theta]d=\[Theta]f'[\[Lambda]]
},
{(*(Subscript[e, i])^(a)Subscript[\[Eta], (a)(b)]Subscript[(\[Omega]^(b)), \[Mu]]*)
MarckCarterFrame[{a_,\[ScriptCapitalE]_,\[ScriptCapitalL]_,\[ScriptCapitalK]_},{r_,rd_,\[Theta]_,\[Theta]d_}]:=
Module[{\[Alpha],\[Beta],\[CapitalSigma],\[CapitalDelta]},
{(*(Subscript[e, i])^(a)Subscript[\[Eta], (a)(b)]Subscript[(\[Omega]^(b)), \[Mu]]*)
{-\[ScriptCapitalE], rd/\[CapitalDelta][r],\[Theta]d, \[ScriptCapitalL]},
{(-r rd \[Alpha]-a^2 \[Beta] \[Theta]d Cos[\[Theta]] Sin[\[Theta]])/(Sqrt[\[ScriptCapitalK]] \[CapitalSigma][r,\[Theta]]),(r ((a^2+r^2) \[ScriptCapitalE]-a \[ScriptCapitalL]) \[Alpha])/(Sqrt[\[ScriptCapitalK]] \[CapitalDelta][r]),(a \[Beta] Cos[\[Theta]] (-\[ScriptCapitalL] Csc[\[Theta]]+a \[ScriptCapitalE] Sin[\[Theta]]))/Sqrt[\[ScriptCapitalK]],(a Sin[\[Theta]] ((a^2+r^2) \[Beta] \[Theta]d Cos[\[Theta]]+r rd \[Alpha] Sin[\[Theta]]))/(Sqrt[\[ScriptCapitalK]] \[CapitalSigma][r,\[Theta]])},
{(a^2 \[ScriptCapitalE] \[Alpha]+r^2 \[ScriptCapitalE] \[Alpha]+a \[ScriptCapitalL] (-\[Alpha]+\[Beta])-a^2 \[ScriptCapitalE] \[Beta] Sin[\[Theta]]^2)/\[CapitalSigma][r,\[Theta]],-((rd \[Alpha])/\[CapitalDelta][r]),-\[Beta] \[Theta]d,(-(a^2+r^2) \[ScriptCapitalL] \[Beta]+a (a \[ScriptCapitalL] \[Alpha]+a^2 \[ScriptCapitalE] (-\[Alpha]+\[Beta])+r^2 \[ScriptCapitalE] (-\[Alpha]+\[Beta])) Sin[\[Theta]]^2)/\[CapitalSigma][r,\[Theta]]},
Expand All @@ -199,8 +187,41 @@
\[Alpha]-> Sqrt[(\[ScriptCapitalK]-a^2 Cos[\[Theta]]^2)/(r^2+\[ScriptCapitalK])],
\[Beta]-> 1/Sqrt[((\[ScriptCapitalK]-a^2 Cos[\[Theta]]^2)/(r^2+\[ScriptCapitalK]))]
}
]


MarckCarterFrame[orbit_KerrGeoOrbitFunction]:=
Module[
{
a,\[ScriptCapitalE],\[ScriptCapitalL],\[ScriptCapitalQ],\[ScriptCapitalK],rf,\[Theta]f,mcf,\[CapitalUpsilon]r,\[CapitalUpsilon]\[Theta]
},
a=orbit["a"];
{\[ScriptCapitalE],\[ScriptCapitalL],\[ScriptCapitalQ]}=Values@orbit["ConstantsOfMotion"];
\[ScriptCapitalK]=\[ScriptCapitalQ]+(a \[ScriptCapitalE]-\[ScriptCapitalL])^2;
\[CapitalUpsilon]r=orbit["Frequencies"]["\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(r\)]\)"];
\[CapitalUpsilon]\[Theta]=orbit["Frequencies"]["\!\(\*SubscriptBox[\(\[CapitalUpsilon]\), \(\[Theta]\)]\)"];
rf=orbit["Trajectory"][[2]];
\[Theta]f=orbit["Trajectory"][[3]];
mcf=Switch[
orbit["Parametrization"],
"Mino",
Function[{Global`\[Lambda]},
Evaluate@With[{
r=rf[Global`\[Lambda]],\[Theta]=\[Theta]f[Global`\[Lambda]],rd=rf'[Global`\[Lambda]],\[Theta]d=\[Theta]f'[Global`\[Lambda]]
},
MarckCarterFrame[{a,\[ScriptCapitalE],\[ScriptCapitalL],\[ScriptCapitalK]},{r,rd,\[Theta],\[Theta]d}]
]
],
"Phases",
Function[{Global`qr,Global`q\[Theta]},
Evaluate@With[{
r=rf[Global`qr],\[Theta]=\[Theta]f[Global`q\[Theta]],rd=\[CapitalUpsilon]r rf'[Global`qr],\[Theta]d=\[CapitalUpsilon]\[Theta] \[Theta]f'[Global`q\[Theta]]
},
MarckCarterFrame[{a,\[ScriptCapitalE],\[ScriptCapitalL],\[ScriptCapitalK]},{r,rd,\[Theta],\[Theta]d}]
]
]
];
mcf
]


Expand All @@ -211,20 +232,51 @@
\[CapitalUpsilon]\[Psi] = MinoPrecessionFrequency[orbit];
pt\[Psi]=PrecessionPhase[orbit,\[CapitalUpsilon]\[Psi],initPhases[[2]],initPhases[[3]],initPhases[[5]]];

ptf=Function[{\[Lambda]},
ptf=Function[{Global`\[Lambda]},
Evaluate[
{
mcf[Global`\[Lambda]][[1]],
mcf[Global`\[Lambda]][[2]]Cos[pt\[Psi][Global`\[Lambda]]]+mcf[Global`\[Lambda]][[3]]Sin[pt\[Psi][Global`\[Lambda]]],
-mcf[Global`\[Lambda]][[2]]Sin[pt\[Psi][Global`\[Lambda]]]+mcf[Global`\[Lambda]][[3]]Cos[pt\[Psi][Global`\[Lambda]]],
mcf[Global`\[Lambda]][[4]]
}
]
];

assoc = Last@orbit;
assoc["PrecessionFrequency"]=\[CapitalUpsilon]\[Psi];
assoc["ParallelTransportedFrame"]=ptf;
assoc["MarckCarterFrame"]=mcf;
assoc["PrecessionPhase"]=pt\[Psi];

KerrParallelTransportFrameFunction[a,p,e,x,assoc]

]


KerrParallelTransportFramePhases[a_,p_,e_,x_]:=Module[
{orbit,mcf,\[CapitalUpsilon]\[Psi],pt\[Psi],ptf,assoc},
orbit=KerrGeoOrbit[a,p,e,x,"Method"->"Analytic","Parametrization"-> "Phases"];
mcf=Evaluate@MarckCarterFrame[orbit];
\[CapitalUpsilon]\[Psi] = MinoPrecessionFrequency[orbit];
pt\[Psi]=Function[{Global`qr,Global`q\[Theta],Global`q\[Psi]},Evaluate@PrecessionPhase[orbit,\[CapitalUpsilon]\[Psi],Global`qr,Global`q\[Theta],Global`q\[Psi]][0]];

ptf=Function[{Global`qr,Global`q\[Theta],Global`q\[Psi]},
Evaluate[
{
mcf[\[Lambda]][[1]],
mcf[\[Lambda]][[2]]Cos[pt\[Psi][\[Lambda]]]+mcf[\[Lambda]][[3]]Sin[pt\[Psi][\[Lambda]]],
-mcf[\[Lambda]][[2]]Sin[pt\[Psi][\[Lambda]]]+mcf[\[Lambda]][[3]]Cos[pt\[Psi][\[Lambda]]],
mcf[\[Lambda]][[4]]
mcf[Global`qr,Global`q\[Theta]][[1]],
mcf[Global`qr,Global`q\[Theta]][[2]]Cos[pt\[Psi][Global`qr,Global`q\[Theta],Global`q\[Psi]]]+mcf[Global`qr,Global`q\[Theta]][[3]]Sin[pt\[Psi][Global`qr,Global`q\[Theta],Global`q\[Psi]]],
-mcf[Global`qr,Global`q\[Theta]][[2]]Sin[pt\[Psi][Global`qr,Global`q\[Theta],Global`q\[Psi]]]+mcf[Global`qr,Global`q\[Theta]][[3]]Cos[pt\[Psi][Global`qr,Global`q\[Theta],Global`q\[Psi]]],
mcf[Global`qr,Global`q\[Theta]][[4]]
}
]
];

assoc = Last@orbit;
assoc["PrecessionFrequency"]=\[CapitalUpsilon]\[Psi];
assoc["ParallelTransportedFrame"]=ptf;
assoc["MarckCarterFrame"]=mcf;
assoc["PrecessionPhase"]=pt\[Psi];

KerrParallelTransportFrameFunction[a,p,e,x,assoc]

Expand All @@ -245,11 +297,12 @@
method = OptionValue["Method"];
param = OptionValue["Parametrization"];

If[param =!= "Mino", Print["Only Mino time parametrization has been implemented for parallel transport."];Return[];];
If[param =!= "Mino"&& param =!= "Phases", Print["Only Mino time parametrization has been implemented for parallel transport."];Return[];];
Copy link
Collaborator

Choose a reason for hiding this comment

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

We should change these Print statements to Messages

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'll take care of this (and other) changes to ParallelTransport.m

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Done

Copy link
Collaborator

Choose a reason for hiding this comment

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

Another instance of a Print statement that we probably want to change to a Message

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Done

If[method =!= "Analytic", Print["Only analytic method has been implemented for parallel transport."];Return[];];

If[method == "Analytic",
If[param == "Mino", Return[KerrParallelTransportFrameMino[a, p, e, x, initPhases]]];
If[param == "Phases", Return[KerrParallelTransportFramePhases[a, p, e, x]]];
];

Print["Unrecognized method: " <> method];
Expand Down Expand Up @@ -277,6 +330,7 @@


KerrParallelTransportFrameFunction[a_, p_, e_, x_, assoc_][\[Lambda]_/;StringQ[\[Lambda]] == False] := assoc["ParallelTransportedFrame"][\[Lambda]]
KerrParallelTransportFrameFunction[a_, p_, e_, x_, assoc_][\[Lambda]__] := assoc["ParallelTransportedFrame"][\[Lambda]]
KerrParallelTransportFrameFunction[a_, p_, e_, x_, assoc_][y_?StringQ] := assoc[y]
Copy link
Collaborator

Choose a reason for hiding this comment

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

We should add Keys support to KerrParallelTransportFrameFunction, just as we have for KerrGeoOrbitFunction

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Done



Expand Down