Code C3: Schwarzschild geodesics
Last edited by Jens Boos (jens.boos@kit.edu), Oct 6, 2024.
(For definitions and conventions refer to the lecture notes.)
Website: http://www.spintwo.net/Courses/EMRI-101/

(* Overall units: GM = 1 *)

(* particle energy and angular momentum *)
EE=0.97;
LL=2Sqrt[3]+0.3;

(* bound orbit: this quantity has to be negative *)
check=EE^2-2(4-LL(LL-Sqrt[LL^2-12]))^2/(LL(LL-Sqrt[LL^2-12])^3)

(* basic metric function f = f(r) and its derivative fp = f'(r) *)
f:=1-2/r[τ];
fp:=2/r[τ]^2;

(* Christoffel symbols, index position: Γttr = (Subscript[((OverTilde[Γ])^t), tr]) etc *)
Γttr:=fp/(2f);
Γrtt:=f fp/2;
Γrrr:=-fp/(2f);
Γrθθ:=-f r[τ];
Γrϕϕ:=-f r[τ]Sin[θ[τ]]^2;
Γθrθ:=1/r[τ];
Γϕrϕ:=1/r[τ];
Γθϕϕ:=-Sin[θ[τ]]Cos[θ[τ]];
Γϕθϕ:=Cos[θ[τ]]/Sin[θ[τ]];

(* simulation duration *)
τmax = 10000;

(* initial radial position (in units of GM=1) *)
r0=10;

(* convert these into initial 4-momentum *)
tdot0=EE/(1-2/r0);
rdot0=Sqrt[EE^2-(1+LL^2/r0^2)(1-2/r0)];
ϕdot0=LL/r0^2;

(* SANITY CHECK: verify that initial conditions describe timelike particle *)
test=(-1)(1-2/r0)tdot0^2+rdot0^2/(1-2/r0)+r0^2ϕdot0^2//N

(* integrate EOMs *)
sol=NDSolve[
{t''[τ]+2Γttr t'[τ]r'[τ]==0,
r''[τ]+Γrtt t'[τ]^2+Γrrr r'[τ]^2+Γrθθ θ'[τ]^2+Γrϕϕ ϕ'[τ]^2==0,
θ''[τ]+2Γθrθ r'[τ]θ'[τ]+Γθϕϕ ϕ'[τ]^2==0,
ϕ''[τ]+2Γϕrϕ r'[τ]ϕ'[τ]+2Γϕθϕ θ'[τ]ϕ'[τ]==0,
t[0]==0,r[0]==r0,θ[0]==π/2,ϕ[0]==0,
t'[0]==tdot0,r'[0]==rdot0,θ'[0]==0,ϕ'[0]==ϕdot0
},{t,r,θ,ϕ},{τ,0,τmax},PrecisionGoal->25
][[1]];
rmax=3r0;
curve[ττ_]:={r[τ]Cos[ϕ[τ]],r[τ]Sin[ϕ[τ]]}/.sol/.τ->ττ;
p[ττmax_]:=ParametricPlot[curve[τ],{τ,0,ττmax},AspectRatio->1,PlotStyle->{Opacity[0.5]},PlotRange->{{-rmax,rmax},{-rmax,rmax}},Frame->True,GridLines->Automatic,FrameLabel->{"x/GM","y/GM"},PlotPoints->1000];
movie[τ0_]:=Show[p[τ0],Graphics[{Disk[{0,0},1]}],Graphics[{Red,Disk[curve[τ0],0.75 (rmax/60)]}]];
Animate[movie[τ],{τ,0.01,τmax},AnimationRate->1000]