Schwarzschild Geodesics

Schwarzschild Radius

(Subscript[r, s])=1

Schwarzschild Metric in Isotropic coordinates

g[{t_,x_,y_,z_}] := Module[{ R,wm,wp},
R=(Sqrt[(x^2)+(y^2)+(z^2)]);
wm=1-((Subscript[r, s])/(4R));
wp=1+((Subscript[r, s])/(4R));
{{(wm^2)/(wp^2),0,0,0},{0,-(wp^4),0,0},{0,0,-(wp^4),0},{0,0,0,-(wp^4)}}
]//FullSimplify

g[{t,x,y,z}] //MatrixForm

Inverse of the Metric

ginv[{t_,x_,y_,z_}] := g[{t,x,y,z}]//Inverse//FullSimplify

ginv[{t,x,y,z}] //MatrixForm

Derivatives of the Metric

gderivs[{t_,x_,y_,z_}]=D[g[{t,x,y,z}],{{t,x,y,z}}]//FullSimplify;

gderivs[{t,x,y,z}]//MatrixForm

Christoffel Symbols of the Metric

Γ[{t_,x_,y_,z_}]=With[{gd = gderivs[{t,x,y,z}]},
Table[Sum[(1/2)ginv[{t,x,y,z}][[a,d]](gd[[c,d,b]]+gd[[b,d,c]]-gd[[b,c,d]]),{d,4}],{a,4},{b,4},{c,4}]]//FullSimplify;

MatrixForm[Γ[{t,x,y,z}],TableDirections->{Row,Row,Column}]

xs={t,x,y,z};

Geodesic Equation ODEs

odes=Table[
0==xs[[a]]''[s]+Sum[Γ[{t[s],x[s],y[s],z[s]}][[a,b,c]]xs[[b]]'[s]xs[[c]]'[s],{b,4},{c,4}],
{a,4}
]//Simplify;
odes//MatrixForm

Geodesic

F=15800;

plot[initpos_,initvel_,S_]:=Module[
{solve,initial,xyz,plotrange,ticks,faces,facegrids,colorfunction},
init={
{t[0],x[0],y[0],z[0]}==initpos,
{t'[0],x'[0],y'[0],z'[0]}==initvel
};
solve=NDSolve[Append[odes,init],{t,x,y,z},{s,0,S}]//Quiet;
xyz[s_]:= Evaluate[{x[s],y[s],z[s]}/.solve]//First;
plotrange={{-50,50},{-50,50},{-50,50}};
ticks={-50,-25,0,25,50};
faces={{0,1,0},{-1,0,0},{0,0,-1}};
facegrids={#,{ticks,ticks}}&/@faces;
colorfunction = Function[{x,y,z,s},{Opacity[0.1+0.9LogisticSigmoid[0.001*(s*S-(S-1000))]],Hue[s*S/F]}];
Show[
 ParametricPlot3D[
xyz[s],
{s,0,S},
PlotRange->plotrange,
SphericalRegion->True,
AxesLabel->{"x","y","z"},
ColorFunction-> colorfunction ,
FaceGrids->facegrids,
Ticks->{ticks},
Background->None
],
 Graphics3D[{Black,Sphere[{0,0,0},(Subscript[r, s])/4]}]
]
];

Manipulate[
plot[{0,40,0,0},{1,vx,vy,vz},S],
{{vx,0},-0.5,0.5},{{vy,0.08},-0.5,0.5},{{vz,0},-0.5,0.5},{{S,F/5},0,F}]

fps=30;  (* frames per seconds *)
duration = 10; (* in seconds *)
animationplots = Table[plot[{0,40,0,0},{1,0,0.08,0},S],{S,1,F,F/(fps*duration)}];

Export["animate.gif", animationplots]