Kerr Black Hole orbits

Newbie first post, not sure how to pretty print the code. Graph the orbits of a Kerr Black Hole. Here's a Wolfram Language version in 22 lines:

Manipulate[(*Initialization*)If[!sE,{pT,a0,r0,L0,\[Theta]0,p\[Theta]0,f0,tL,zM}=pv;sE=True;];
(*Viewpoints setup*)vp=10;v\[Theta]=0.85 \[Pi]/2;v\[Phi]=0.35 \[Pi]/2;d=0.05 \[Pi]/2;
rvp=vp {Sin[v\[Theta]] Cos[v\[Phi]],Sin[v\[Theta]] Sin[v\[Phi]],Cos[v\[Theta]]};
lvp=vp {Sin[v\[Theta]] Cos[v\[Phi]-d],Sin[v\[Theta]] Sin[v\[Phi]-d],Cos[v\[Theta]]};
(*Energy calculation*)E0=\|01d486/. Solve[(-a0^2 p\[Theta]0^2+2 L0^2 r0+2 p\[Theta]0^2 r0-a0^2 r0^2-L0^2 r0^2-p\[Theta]0^2 r0^2+2 r0^3-r0^4-4 a0 L0 r0 \|01d486+2 a0^2 r0 \|01d486^2+a0^2 r0^2 \|01d486^2+r0^4 \|01d486^2+a0^2 (a0^2+(-2+r0) r0) (-1+\|01d486^2) Cos[\[Theta]0]^2-L0^2 (a0^2+(-2+r0) r0) Cot[\[Theta]0]^2)==0,\|01d486][[2]];
(*Carter constant*)Q=p\[Theta]0^2+Cos[\[Theta]0]^2 (a0^2 (1-E0^2)+L0^2 Csc[\[Theta]0]^2);
(*Equations of motion*)dynEq={r'[\[Tau]]==(pr[\[Tau]] (a^2-2 r[\[Tau]]+r[\[Tau]]^2))/(a^2 Cos[\[Theta][\[Tau]]]^2+r[\[Tau]]^2),pr'[\[Tau]]==(a^4 (-a E0+L)^2 Cos[\[Theta][\[Tau]]]^2+a^4 (L^2 Cos[\[Theta][\[Tau]]]^2 Cot[\[Theta][\[Tau]]]^2+p\[Theta][\[Tau]]^2) r[\[Tau]]+a^2 (-a^2 E0^2+2 a E0 L-L^2+2 a E0 (a E0+L) Cos[\[Theta][\[Tau]]]^2-4 L^2 Cot[\[Theta][\[Tau]]]^2-4 p\[Theta][\[Tau]]^2) r[\[Tau]]^2+(4 a^2 E0^2-8 a E0 L+4 L^2-4 a^2 E0^2 Cos[\[Theta][\[Tau]]]^2+4 L^2 Cot[\[Theta][\[Tau]]]^2+2 a^2 L^2 Cot[\[Theta][\[Tau]]]^2+2 (2+a^2) p\[Theta][\[Tau]]^2) r[\[Tau]]^3+(-2 a^2 E0^2+6 a E0 L-4 L^2+a^2 E0^2 Cos[\[Theta][\[Tau]]]^2-4 L^2 Cot[\[Theta][\[Tau]]]^2-4 p\[Theta][\[Tau]]^2) r[\[Tau]]^4+(L^2 Csc[\[Theta][\[Tau]]]^2+p\[Theta][\[Tau]]^2) r[\[Tau]]^5-E0^2 r[\[Tau]]^6+pr[\[Tau]]^2 (a^2-2 r[\[Tau]]+r[\[Tau]]^2)^2 (a^2 Cos[\[Theta][\[Tau]]]^2-r[\[Tau]]^2+a^2 r[\[Tau]] Sin[\[Theta][\[Tau]]]^2))/((a^2 Cos[\[Theta][\[Tau]]]^2+r[\[Tau]]^2)^2 (a^2-2 r[\[Tau]]+r[\[Tau]]^2)^2),\[Phi]'[\[Tau]]==(a^2 L Cot[\[Theta][\[Tau]]]^2+2 (a E0-L-L Cot[\[Theta][\[Tau]]]^2) r[\[Tau]]+L Csc[\[Theta][\[Tau]]]^2 r[\[Tau]]^2)/((a^2 Cos[\[Theta][\[Tau]]]^2+r[\[Tau]]^2) (a^2-2 r[\[Tau]]+r[\[Tau]]^2)),\[Theta]'[\[Tau]]==p\[Theta][\[Tau]]/(a^2 Cos[\[Theta][\[Tau]]]^2+r[\[Tau]]^2),p\[Theta]'[\[Tau]]==((2 a^2 Cos[\[Theta][\[Tau]]] ((Q+a^2 (-1+E0^2) Cos[\[Theta][\[Tau]]]^2-L^2 Cot[\[Theta][\[Tau]]]^2) (a^2-2 r[\[Tau]]+r[\[Tau]]^2)-(Q+(-a E0+L)^2+r[\[Tau]]^2) (a^2-2 r[\[Tau]]+r[\[Tau]]^2)+(a L-E0 (a^2+r[\[Tau]]^2))^2) Sin[\[Theta][\[Tau]]])/(a^2-2 r[\[Tau]]+r[\[Tau]]^2)-a^2 p\[Theta][\[Tau]]^2 Sin[2 \[Theta][\[Tau]]]-a^2 pr[\[Tau]]^2 (a^2-2 r[\[Tau]]+r[\[Tau]]^2) Sin[2 \[Theta][\[Tau]]]+(a^2 Cos[\[Theta][\[Tau]]]^2+r[\[Tau]]^2) (2 L^2 Cot[\[Theta][\[Tau]]]+2 L^2 Cot[\[Theta][\[Tau]]]^3-a^2 (-1+E0^2) Sin[2 \[Theta][\[Tau]]]))/(2 (a^2 Cos[\[Theta][\[Tau]]]^2+r[\[Tau]]^2)^2)};
(*Initial conditions*)iC={r[0]==r0,pr[0]==0,\[Theta][0]==\[Theta]0,p\[Theta][0]==p\[Theta]0,\[Phi][0]==0};
(*Solve equations of motion*)Quiet[HSolve=NDSolve[Flatten[{dynEq/. {a->a0,L->L0},iC}],{r,\[Phi],\[Theta],pr,p\[Theta]},{\[Tau],0,pT},Method->{EventLocator,"Event"->(r[\[Tau]]-1.02 hs)}]];
dom=(r/. HSolve[[1]])["Domain"][[1]];
e=dom[[2]];
pHP=Abs[(r[e]/. HSolve)[[1]]-hs]<=0.05 hs;
sP=If[e-tL<=0,0,e-tL];
If[!zM,oR=1.05 If[(r[e]/. HSolve)[[1]]>r0,(r[e]/. HSolve)[[1]],r0];f0=If[oR>f0,oR,f0];];
pp={r[e] Sin[\[Theta][e]] Cos[\[Phi][e]],r[e] Sin[\[Theta][e]] Sin[\[Phi][e]],r[e] Cos[\[Theta][e]]}/. HSolve;
oPlot=ParametricPlot3D[{r[\[Tau]] Sin[\[Theta][\[Tau]]] Cos[\[Phi][\[Tau]]],r[\[Tau]] Sin[\[Theta][\[Tau]]] Sin[\[Phi][\[Tau]]],r[\[Tau]] Cos[\[Theta][\[Tau]]]}/. HSolve,{\[Tau],sP,e},PlotRange->Table[{-f0,f0},{3}],PerformanceGoal->"Speed",PlotPoints->200,MaxRecursion->8,SphericalRegion->True,Mesh->4,Ticks->Automatic];
hs=1+Sqrt[1-a0^2];pS=0.02 f0;apS=If[pHP,0,pS];
(*Graphics elements*)graphics={Green,Sphere[pp,apS],Black,Sphere[{0,0,0},hs],Black,Opacity[0.2],Scale[Sphere[],{2,2,hs},{0,0,0}],Text[StringForm["E=``",E0],{1.5 f0,0,-1.1 f0}],Text[StringForm["Q=``",Chop[Q]],{1.5 f0,0,-1.3 f0}]};
(*Output:two plots in a column*)Column[{Show[oPlot,Graphics3D[graphics],ViewPoint->rvp,ImageSize->400],Show[oPlot,Graphics3D[graphics],ViewPoint->lvp,ImageSize->400]}],(*Controls*)Grid[{{Control[{{pT,150,"time"},150,50000,ImageSize->Tiny}],SpanFromLeft},{Control[{{a0,0.99,"spin"},0,0.99,0.01,Appearance->"Labeled",ImageSize->Tiny}],SpanFromLeft},{Control[{{r0,4,"r"},2.1,30,0.01,Appearance->"Labeled",ImageSize->Tiny}],Control[{{L0,2,"L"},-4.5,4.5,0.01,Appearance->"Labeled",ImageSize->Tiny}]},{Control[{{\[Theta]0,\[Pi]/3,"\[Theta]I"},-\[Pi]/1080,6\[Pi]/7,\[Pi]/1080,Appearance->"Labeled",ImageSize->Tiny}],Control[{{p\[Theta]0,0.76,"p\[Theta]I"},-3,3,0.01,Appearance->"Labeled",ImageSize->Tiny}]},{Control[{{tL,1200,"tail"},150,50000,ImageSize->Tiny}],SpanFromLeft},{Control[{{f0,4.5,"zoom"},-2.5,100,0.01,Appearance->"Labeled",ImageSize->Tiny,Enabled->zM}],Control[{{zM,False},{False->"auto",True->"manual"}}]},{Control[{{sE,True},{False->"orbit preset"},SetterBar,ImageSize->Tiny}],Control[{{pv,{300,0.9,4,2.148,1.037,0,4.2,350,False}},{{300,0.9,4,2.148,1.037,0,4.2,350,False}->"closed orbit",{1200,0.99,4,2,\[Pi]/3,0.767851,4.5,1200,False}->"constant radius orbit",{150,0.0,10,3.5,\[Pi]/2,0,4.5,350,False}->"spiral capture orbit",{150,0.0,4,3.99999,\[Pi]/2,0,4.5,350,False}->"unstable circular capture",{100,0.0,4,4.00001,\[Pi]/2,0,4.5,350,False}->"unstable circular escape",{330,0.99,25,2.427,\[Pi]/2,0,25,330,False}->"equatorial zoom",{150,0.9,4,-4.5,\[Pi]/2,0,4.2,350,False}->"reverse capture",{150,0.99,10,1.05769,\[Pi]/2,2.89,4,150,True}->"3D whirl"},PopupMenu,ImageSize->Small}]}},Alignment->Left],TrackedSymbols:>{pT,a0,r0,L0,\[Theta]0,p\[Theta]0,f0,tL,zM,pv,sE},SynchronousUpdating->False,AutorunSequencing->{1,2,3,4,6,7}]