def.

SetDirectory[NotebookDirectory[]];
$PlotTheme=None;(* plots in Mathematica10 looks like in previous versions *)

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Off[NDSolve::"ndsz"](*vypnuti hlasky pri dosazeni horizontu*)
Off[NDSolve::maxit]
Off[NDSolve::"noout"]

fonty=15;fontyW={FontSize->15,Background->White};xInfinity=1000;
G=6.67300×(10^(-11));c=299792458 ;sM=1.98892×(10^30) ;
ImSz=300;
ImRez=200;

tA[]:=(CasZ=SessionTime[];DateZ=DateString[{ "Start ","Time", " ","DayName", " ","Day", "/","Month"}];Date1=Date[];);
tB[]:=(CasK=SessionTime[];DateK=DateString[{ "Konec ","Time", " ","DayName", " ","Day", "/","Month"}];Date2=Date[];
Print[DateString[Date2-Date1,{ "TIME"," ","Time"," ","Day"}]<>"  #  "<>DateZ<>" / "<>DateK];);

(* 2D - polar coordinates *)
XX[r_,ϕ_]:=r Cos[ϕ];YY[r_,ϕ_]:=r Sin[ϕ];
mXX[r_,θ_]:=(Sqrt[(r^2)+(aa^2)]) Sin[θ];mZZ[r_,θ_]:=r Cos[θ];

(* 3D - spherical coordinates *)
XXX[r_,θ_,ϕ_]:=r Sin[θ] Cos[ϕ];
YYY[r_,θ_,ϕ_]:=r Sin[θ] Sin[ϕ];
ZZZ[r_,θ_,ϕ_]:=r Cos[θ];

mXXX[r_,θ_,ϕ_]:=(Sqrt[(r^2)+(aa^2)]) Sin[θ] Cos[ϕ-ArcTan[aa/r]];
mYYY[r_,θ_,ϕ_]:=(Sqrt[(r^2)+(aa^2)])Sin[θ] Sin[ϕ-ArcTan[aa/r]];
mZZZ[r_,θ_,ϕ_]:=r Cos[θ];

particle moving around black hole

M=1;

(*velka sigma*) R2[r_,θ_,a_]:=(r^2)+(a^2)((Cos[θ])^2);(*trojuhelnicek*) TR[r_,a_]:=(r^2)-2M r +(a^2);

(*metrika s dolnimi indexy*)
 gtt[r_,θ_,a_]:=-(1-((2M r)/(R2[r,θ,a])));gtϕ[r_,θ_,a_]:=-((2 r M a ((Sin[θ])^2))/( R2[r,θ,a]));gϕϕ[r_,θ_,a_]:=((r^2)+(a^2)+((2  M r (a^2) )/( R2[r,θ,a]))((Sin[θ])^2))((Sin[θ])^2);
grr[r_,θ_,a_]:=((R2[r,θ,a])/(TR[r,a]));gθθ[r_,θ_,a_]:= R2[r,θ,a];
(*metrika s hornimi indexy*) 
gUtt[r_,θ_,a_]:=-(1/(TR[r,a]))((r^2)+(a^2)+((2  M r (a^2) )/( R2[r,θ,a]))((Sin[θ])^2));
gUtϕ[r_,θ_,a_]:=-((2 M a r)/(TR[r,a] R2[r,θ,a]));gUϕϕ[r_,θ_,a_]:=((TR[r,a]- (a^2)((Sin[θ])^2))/(TR[r,a]R2[r,θ,a]((Sin[θ])^2)));gUrr[r_,θ_,a_]:=((TR[r,a])/(R2[r,θ,a]));gUθθ[r_,θ_,a_]:=(1/(R2[r,θ,a]));


(*effective potential*) 
amA[r_,θ_,a_,LL_]:=-gUtt[r,θ,a];
amB[r_,θ_,a_,LL_]:=2gUtϕ[r,θ,a]LL;
amC[r_,θ_,a_,LL_]:=-gUϕϕ[r,θ,a](((LL))^2)-1;
Veff[r_,θ_,a_,LL_]:=((-amB[r,θ,a,LL]+(Sqrt[((amB[r,θ,a,LL])^2)-4amA[r,θ,a,LL]amC[r,θ,a,LL]]))/(2amA[r,θ,a,LL]));

RfromXZ[x_,z_,a_]:=(Sqrt[(1/2) (-(a^2)+(x^2)+(z^2)+(Sqrt[4 (a^2) (z^2)+((((a^2)-(x^2)-(z^2)))^2)]))]);TfromXZ[x_,z_,a_]:=ArcTan[(Sqrt[((-(a^2)+(x^2)-(z^2)+(Sqrt[4 (a^2) (z^2)+(((-(a^2)+(x^2)+(z^2)))^2)]))/(2 (z^2)))])];VeffXZ[x_,z_,a_,L_]:=Veff[RfromXZ[x,z,a],TfromXZ[x,z,a],a,L];


GiveMeTrajectory[R0_,θ0_,ϕ0_,aa_,LL_,xEND_]:=(
Clear[r,pr,θ, pθ,ϕ,EE];

RH1=Part[r/.NSolve[(r^2)-2r+(aa^2)==0,r],1];RH2=Part[r/.NSolve[(r^2)-2r+(aa^2)==0,r],2];If[Element[RH2,Reals],RH=RH2,RH=-1];
Rend =1.001RH;

EE= Veff[R0,θ0,aa,LL];(*particle energy*)(* works only for zero momenta - STARTUJEME Z KLIDU - *)

(*---- Hamiltonian a equations of motion ----*)HoHo[r_,θ_,pr_,pθ_]:=(1/2)gUrr[r,θ,aa](pr^2)+(1/2)gUθθ[r,θ,aa]((pθ)^2)+(1/2)+(1/2)( gUtt[r,θ,aa](((-EE))^2)+2gUtϕ[r,θ,aa](-EE)(LL)+gUϕϕ[r,θ,aa](((LL))^2));

(* pr = (Subscript[p, r]) - dolni indexy *) (* r = (x^r), gUrr = (g^rr) - horni indexy *)
DrHoHo[r_,θ_,pr_,pθ_]:=(∂_r)HoHo[r,θ,pr,pθ];DθHoHo[r_,θ_,pr_,pθ_]:=(∂_θ)HoHo[r,θ,pr,pθ];HamI=HoHo[r[τ],θ[τ],pr[τ],pθ[τ]];

Equations ={
r'[τ]   ==    gUrr[r[τ] ,θ[τ],aa ]pr[τ],pr'[τ]  == -DrHoHo[r[τ],θ[τ],pr[τ],pθ[τ]],
θ'[τ]   ==   gUθθ[r[τ] ,θ[τ] ,aa]pθ[τ],pθ'[τ]  ==-DθHoHo[r[τ],θ[τ],pr[τ],pθ[τ]],
ϕ'[τ]   ==   -gUtϕ[r[τ] ,θ[τ],aa]EE+ gUϕϕ[r[τ] ,θ[τ],aa]LL
(* ϕ' = (p^ϕ) = (g^(ϕt))(Subscript[p, t]) + (g^(ϕϕ))(Subscript[p, ϕ]) = -(g^(ϕt))E + (g^(ϕϕ))L *) 
(* t' = we do not need to calculate, becouse the model is nedependent of t and ϕ (stacionarity, axill symmetry) and t coordinate is not used for plots (as ϕ does) *)
};
IntCond ={r[0]    == R0, pr[0]   ==  0,θ[0]   == θ0, pθ[0]  == 0,ϕ[0]   == ϕ0};

(*---- Hamiltonian a equations of motion ----*)

(*---- trajectory calculation ----*)
Metoda={"StiffnessSwitching", Method->{ {"Projection", Method->{"ExplicitRungeKutta","DifferenceOrder"->8},"Invariants"->HamI},Automatic}};
(* put your favourite numerical method here *)

draha=NDSolve[{Equations,IntCond},{r,pr,θ,pθ,ϕ},{τ,0,xEND},{Method->{"EventLocator","Event"->{r[τ]-Rend,r[τ]-xInfinity},"EventAction":>{Throw[ Null,"StopIntegration"],Throw[Null, "StopIntegration"]},"Method"->Metoda},{StartingStepSize->(10^(-1)),MaxSteps->∞}}];

ipub=r/.draha[[1,1]];domain=ipub["Domain"];{begin,end} = domain[[1]];
(*---- trajectory calculation ----*)

{r[τ],pr[τ],θ[τ],pθ[τ] ,ϕ[τ],EE}/.draha /.τ ->end

);

ShowMe[time_]:=(
GiveMeTrajectory[R0,θ0,0,aa,LL,time];

pic1=
ParametricPlot[ { mXXX[r[τ],θ[τ],ϕ[τ]], mYYY[r[τ],θ[τ],ϕ[τ]]}/.draha,{τ,begin,end},PlotStyle->{Blue},
Frame->True,BaseStyle->15,PlotRange->oblast2D,FrameLabel->{"x","z"},RotateLabel->False,Axes->False,ImageSize->ImSz,ImagePadding->{{43,1},{43,1}},Epilog->{LightGray,Disk[{0,0},RH],Black,PointSize[0.03],Point[{mXXX[R0,θ0,0], mYYY[R0,θ0,0]}]}
];

pic2=Show[
ParametricPlot[ {mXX[r[τ],θ[τ]],mZZ[r[τ],θ[τ]] }/.draha,{τ,begin,end},PlotStyle->{Blue}],
Quiet[ContourPlot[ VeffXZ[xx,zz,aa,LL]==EE,{xx,0,oblast1D[[2]]},{zz,-oblast1D[[2]]/2,oblast1D[[2]]/2},ContourStyle->{{Red,Thick,Dashed}},RegionFunction->Function[{x,z},RfromXZ[x,z,aa]>RH&&x>0.1]]],
blackhole2D,
Frame->True,BaseStyle->15,PlotRange->{{0,oblast1D[[2]]},{-oblast1D[[2]]/2,oblast1D[[2]]/2}},FrameLabel->{"x","z"},RotateLabel->False,Axes->False,ImageSize->ImSz,ImagePadding->{{43,1},{43,1}},Epilog->{PointSize[0.03],Point[{mXX[R0,θ0],mZZ[R0,θ0]}]}
];

pic3=Show[
ParametricPlot3D[{ mXXX[r[τ],θ[τ],ϕ[τ]], mYYY[r[τ],θ[τ],ϕ[τ]],  mZZZ[r[τ],θ[τ],ϕ[τ]] }/.draha,{τ,begin,end},PlotRange->oblast3D,BaseStyle->{15,Background->White},PlotStyle->{Blue},BoxStyle->{Gray},ViewPoint->{0,-3,1},Ticks->None,Lighting->"Neutral",ImageSize->ImSz],
blackhole3D,
Graphics3D[{sl=oblast3D[[1,2]];
{Text[Style["x",Background->White],{0.8sl,sl,-sl}],Text[Style["y",Background->White],{sl,0.8sl,-sl}],Text[Style["z",Background->White],{sl,sl,-0.8sl}]},
{Text[Framed[Style[Row[{"(Subscript[r, 0])=",R0," | (Subscript[θ, 0])=",θ0," | L≐",Round[LL,0.01]," | E≐",Round[EE,0.01]}],fontyW]],Scaled[{0.52,0.2,0}]],Black,PointSize[0.02],Point[{mXXX[R0,θ0,0],mYYY[R0,θ0,0],mZZZ[R0,θ0,0]}]}}]
];
{pic1,pic2,pic3}
);

Schw. a = 0

aa=0;RH1=Part[r/.NSolve[(r^2)-2r+(aa^2)==0,r],1];RH2=Part[r/.NSolve[(r^2)-2r+(aa^2)==0,r],2];RH=RH2;
oblast1D={-12.5,12.5};oblast2D={oblast1D,oblast1D};oblast3D={oblast1D,oblast1D,oblast1D};
blackhole2D=RegionPlot[RfromXZ[x,z,aa]<RH,{x,oblast1D[[1]],oblast1D[[2]]},{z,oblast1D[[1]],oblast1D[[2]]},PlotStyle->LightGray];
blackhole3D=RegionPlot3D[RfromXZ[(Sqrt[(x^2)+(y^2)]),z,aa]<RH,{x,oblast1D[[1]],oblast1D[[2]]},{y,oblast1D[[1]],oblast1D[[2]]},{z,oblast1D[[1]],oblast1D[[2]]},PlotStyle->Lighter[Gray,0.5],PlotPoints->100,Mesh->None,Lighting->"Neutral"];
{blackhole2D,blackhole3D}

R0=7;θ0=Pi/2-0.2;LL=3.35;ShowMe[300]

R0=7;θ0=Pi/2-0.2;LL=3.6;ShowMe[500]

R0=7;θ0=Pi/2-0.2;LL=4;ShowMe[500]

Kerr  a = 0.9

aa=0.9;RH1=Part[r/.NSolve[(r^2)-2r+(aa^2)==0,r],1];RH2=Part[r/.NSolve[(r^2)-2r+(aa^2)==0,r],2];RH=RH2;
oblast1D={-6.5,6.5};oblast2D={oblast1D,oblast1D};oblast3D={oblast1D,oblast1D,oblast1D};
blackhole2D=RegionPlot[RfromXZ[x,z,aa]<RH,{x,oblast1D[[1]],oblast1D[[2]]},{z,oblast1D[[1]],oblast1D[[2]]},PlotStyle->LightGray];
blackhole3D=RegionPlot3D[RfromXZ[(Sqrt[(x^2)+(y^2)]),z,aa]<RH,{x,oblast1D[[1]],oblast1D[[2]]},{y,oblast1D[[1]],oblast1D[[2]]},{z,oblast1D[[1]],oblast1D[[2]]},PlotStyle->Lighter[Gray,0.5],PlotPoints->100,Mesh->None,Lighting->"Neutral"];
{blackhole2D,blackhole3D}

R0=4;θ0=Pi/2-0.4;LL=2.4;ShowMe[500]

R0=4;θ0=Pi/2-0.4;LL=2;ShowMe[24]