(* Simulator-Code für Photonen, geladene und neutrale Teilchen                            *)
(* in Boyer Lindquist Koordinaten                                                         *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||| Mathematica | kerr.newman.yukterez.net | 06.08.2017 - 13.06.2020, Version 25 |||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
wp=MachinePrecision;
mt1=Automatic;
mt2={"EquationSimplification"-> "Residual"};
mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
mt4={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mt5={"EventLocator", "Event"-> (r[τ]-1001/1000 rA)};
mta=mt1;                                                     (* mt1: Speed, mt3: Accuracy *)
dgl=1;         (* 1: Full d²x/dτ², 2: Mixed dx/dτ, 3: Testmodul, 4: Weak Field, 5: Newton *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
A=a;                                   (* pseudosphärisch [BL]: A=0, kartesisch [KS]: A=a *)
 
tmax=300;                                                                    (* Eigenzeit *)
Tmax=300;                                                              (* Koordinatenzeit *)
TMax=Min[Tmax, т[plunge-1/100]]; tMax=Min[tmax, plunge-1/100];        (* Integrationsende *)
 
r0 = Sqrt[7^2-a^2];                                                        (* Startradius *)
r1 = r0+2;                                                   (* Endradius wenn v0=vr0=vr1 *)
θ0 = π/2;                                                                  (* Breitengrad *)
φ0 = 0;                                                                     (* Längengrad *)
a  = 9/10;                                                               (* Spinparameter *)
℧  = 2/5;                                       (* spezifische Ladung des schwarzen Lochs *)
q  = 0;                                             (* spezifische Ladung des Testkörpers *)
 
v0 = 2/5;                                                       (* Anfangsgeschwindigkeit *)
α0 = 0;                                                      (* vertikaler Abschusswinkel *)
i0 = ArcTan[5/6];                                               (* Bahninklinationswinkel *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 2) GESCHWINDIGKEITS-, ENERGIE UND DREHIMPULSKOMPONENTEN ||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
vr0=v0 Sin[α0];                                     (* radiale Geschwindigkeitskomponente *)
vθ0=v0 Cos[α0] Sin[i0];                      (* longitudinale  Geschwindigkeitskomponente *)
vφ0=v0 Cos[α0] Cos[i0];                        (* latitudinale Geschwindigkeitskomponente *)

vrj[τ_]:=R'[τ]/Sqrt[Δi[τ]] Sqrt[Σi[τ] (1+μ v[τ]^2)];
vθj[τ_]:=Θ'[τ] Sqrt[Σi[τ] (1+μ v[τ]^2)];
vφj[τ_]:=Evaluate[(-(((a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]] Sqrt[1-
μ^2 v[τ]^2] (-φ'[τ]-(a q ℧ r[τ])/((a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2))+
(ε Csc[θ[τ]]^2 (a (-a^2-℧^2+2 r[τ]-r[τ]^2) Sin[θ[τ]]^2+a (a^2+
r[τ]^2) Sin[θ[τ]]^2))/((a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2))+(a q ℧ r[τ] (a^2+
℧^2-2 r[τ]+r[τ]^2-a^2 Sin[θ[τ]]^2))/((a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (a^2+℧^2-2 r[τ]+
r[τ]^2) (1-μ^2 v[τ]^2))))/((a^2+℧^2-2 r[τ]+r[τ]^2-a^2 Sin[θ[τ]]^2) Sqrt[((a^2+r[τ]^2)^2-
a^2 (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2)/(a^2 Cos[θ[τ]]^2+r[τ]^2)]))) /. sol][[1]]
vtj[τ_]:=Sqrt[vrj[τ]^2+vθj[τ]^2+vφj[τ]^2];
vr[τ_]:=vrj[τ]/vtj[τ]*v[τ];
vθ[τ_]:=vθj[τ]/vtj[τ]*v[τ];
vφ[τ_]:=vφj[τ]/vtj[τ]*v[τ];
VΦ[τ_]:=Sqrt[v[τ]^2-vθ[τ]^2-vr[τ]^2];
Vφ[τ_]:=If[q==0, Vφ[τ], VΦ[τ]];
 
x0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0];                             (* Anfangskoordinaten *)
y0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
z0[A_]:=r0 Cos[θ0];
 
ε0=Sqrt[δ Ξ/χ]/j[v0]+Lz ω0;
ε=ε0+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2));
εζ:=Sqrt[Δ Σ/Χ]/j[ν]+Lz ωζ+((q r[τ] ℧)/(r[τ]^2+a^2 Cos[θ[τ]]^2));
LZ=vφ0 Ы/j[v0];
Lz=LZ+((q a r0 ℧ Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)) j[v0]^2;
Lζ:=vφ0 я/j[ν]+0((q a r[τ] ℧ Sin[θ[τ]]^2)/(r[τ]^2+a^2 Cos[θ[τ]]^2));
pθ0=vθ0 Sqrt[Ξ]/j[v0]; pθζ:=θ'[τ] Σ;
pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
Qk=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2, θ1->θ0];       (* Carter Konstante *)
Q=Limit[pθ0^2+(Lz^2 Csc[θ1]^2-a^2 (ε^2+μ)) Cos[θ1]^2, θ1->θ0];
Qζ:=pθζ^2+(Lz^2 Csc[θ[τ]]^2-a^2 (εζ^2+μ)) Cos[θ[τ]]^2;
k=Q+Lz^2+a^2 (ε^2+μ); kζ:=Qζ+Lz^2+a^2 (εζ^2+μ);
                                                                                  (* ISCO *)
isco = rISCO/.Solve[0 == rISCO (6 rISCO-rISCO^2-9 ℧^2+3 a^2)+4 ℧^2 (℧^2-a^2)-
8 a (rISCO-℧^2)^(3/2) && rISCO>=rA, rISCO][[1]];
μ=If[Abs[v0]==1, 0, If[Abs[v0]<1, -1, 1]];                   (* Baryon: μ=-1, Photon: μ=0 *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 3) FLUCHTGESCHWINDIGKEIT UND BENÖTIGTER ABSCHUSSWINKEL |||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

vEsc=If[q==0, ж0, Abs[(\[Sqrt](r0^2 (r0^2 (δ Ξ-χ)+2 q r0 χ ℧-q^2 χ ℧^2)+
2 a^2 r0 (r0 δ Ξ-r0 χ+q χ ℧) Cos[θ0]^2+a^4 (δ Ξ-
χ) Cos[θ0]^4))/(Sqrt[χ] (r0 (r0-q ℧)+a^2 Cos[θ0]^2))]];
                                                  (* horizontaler Photonenkreiswinkel, i0 *)
iP[r0_, a_]:=Normal[iPh/.NSolve[1/(8 (r0^2+a^2 Cos[θ0]^2)^3) (a^2+(-2+r0) r0+
℧^2) (8 r0 (r0^2+a^2 Cos[θ0]^2) Sin[iPh]^2+1/((a^2-2 r0+r0^2+℧^2) (r0^2+
a^2 Cos[θ0]^2)) 8 a (Cos[iPh] Sin[θ0] (a^2-2 r0+r0^2+℧^2-a^2 Sin[θ0]^2) Sqrt[((a^2+
r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+(a (a^2+r0^2) Sin[θ0]^2+
a (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^2) (Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+
r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a (2 r0-℧^2) Cos[iPh] Sin[θ0] Sqrt[((a^2+
r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-
2 r0+r0^2+℧^2) Sin[θ0]^2))) (-(1/((a^2-2 r0+r0^2+℧^2) (r0^2+
a^2 Cos[θ0]^2)))2 a^2 Cot[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-2 r0+
r0^2+℧^2) Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-
℧^2) Sin[θ0]^4) (Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-
2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a (2 r0-℧^2) Cos[iPh] Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-
2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^2)))+1/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)) 2 r0 (r0-
℧^2) Csc[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^4) (Sqrt[((a^2-
2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+
(a (2 r0-℧^2) Cos[iPh] Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))))+1/((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θ0]^2)) 8 Csc[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-
2 r0+r0^2+℧^2) Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^4) (Sqrt[((a^2-
2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+
(a (2 r0-℧^2) Cos[iPh] Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))) (1/((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θ0]^2)) a^2 Cot[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+
a (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-
℧^2) Sin[θ0]^4) (Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-
a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a (2 r0-℧^2) Cos[iPh] Sin[θ0] Sqrt[((a^2+
r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-
2 r0+r0^2+℧^2) Sin[θ0]^2)))+1/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)) r0 (-r0+
℧^2) Csc[θ0]^2 (Cos[iPh] Sin[θ0] (-a (a^2+r0^2) Sin[θ0]^2+a (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+
((a^2+r0^2)^2 Sin[θ0]^2+a^2 (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^4) (Sqrt[((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a (2 r0-
℧^2) Cos[iPh] Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))))+1/((a^2-2 r0+r0^2+
℧^2)^2 (r0^2+a^2 Cos[θ0]^2)^2) Csc[θ0]^2 (Cos[iPh] Sin[θ0] (a^2-2 r0+r0^2+℧^2-
a^2 Sin[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+
(a (a^2+r0^2) Sin[θ0]^2+a (-a^2+2 r0-r0^2-℧^2) Sin[θ0]^2) (Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+
a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a (2 r0-
℧^2) Cos[iPh] Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+
a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)))^2 (r0 (a^2 (3 a^2+
4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ0]+a^2 Cos[4 θ0])+8 r0 (r0^3+2 a^2 r0 Cos[θ0]^2-
a^2 Sin[θ0]^2))+2 a^4 Sin[2 θ0]^2))==0,iPh,Reals]][[1]]/.C[1]->0
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 4) HORIZONTE UND ERGOSPHÄREN RADIEN ||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
rE=1+Sqrt[1-a^2 Cos[θ]^2-℧^2];                                       (* äußere Ergosphäre *)
RE[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rE^2+A^2] Sin[θ] Cos[φ], Sqrt[rE^2+A^2] Sin[θ] Sin[φ], rE Cos[θ]}, w1], w2];
rG=1-Sqrt[1-a^2 Cos[θ]^2-℧^2];                                       (* innere Ergosphäre *)
RG[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2];
rA=1+Sqrt[1-a^2-℧^2];                                                 (* äußerer Horizont *)
RA[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rA^2+A^2] Sin[θ] Cos[φ], Sqrt[rA^2+A^2] Sin[θ] Sin[φ], rA Cos[θ]}, w1], w2];
rI=1-Sqrt[1-a^2-℧^2];                                                 (* innerer Horizont *)
RI[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rI^2+A^2] Sin[θ] Cos[φ], Sqrt[rI^2+A^2] Sin[θ] Sin[φ], rI Cos[θ]}, w1], w2];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 5) HORIZONTE UND ERGOSPHÄREN PLOT ||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
horizons[A_, mesh_, w1_, w2_]:=Show[
ParametricPlot3D[RE[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> mesh, PlotPoints -> plp, PlotStyle -> Directive[Blue, Opacity[0.10]]],
ParametricPlot3D[RA[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
ParametricPlot3D[RI[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.25]]],
ParametricPlot3D[RG[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.35]]]];
BLKS:=Grid[{{horizons[a, 35, 0, 0], horizons[0, 35, 0, 0]}}];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 6) FUNKTIONEN ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
j[v_]:=Sqrt[1-μ^2 v^2];                                                  (* Lorentzfaktor *)
mirr=Sqrt[2-℧^2+2 Sqrt[1-a^2-℧^2]]/2;                                (* irreduzible Masse *)
я=Sqrt[Χ/Σ]Sin[θ[τ]];                                            (* axialer Umfangsradius *)
яi[τ_]:=Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
Ы=Sqrt[χ/Ξ]Sin[θ0];

Σ=r[τ]^2+a^2 Cos[θ[τ]]^2;                                    (* poloidialer Umfangsradius *)
Σi[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2;
Ξ=r0^2+a^2 Cos[θ0]^2;

Δ=r[τ]^2-2r[τ]+a^2+℧^2;
Δr[r_]:=r^2-2r+a^2+℧^2;
Δi[τ_]:=R[τ]^2-2R[τ]+a^2+℧^2;
δ=r0^2-2r0+a^2+℧^2;

Χ=(r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 Δ;
Χi[τ_]:=(R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 Δi[τ];
χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;

Xj=a Sin[θ0]^2;
xJ[τ_]:=a Sin[Θ[τ]]^2;
XJ=a Sin[θ[τ]]^2;

Pr[r_]:=ε(r^2+a^2)+℧ q r-a Lz;
Pt[τ_]:=ε(R[τ]^2+a^2)+℧ q R[τ]-a Lz;
Pτ=ε(r[τ]^2+a^2)+℧ q r[τ]-a Lz;
pτ=ε(r0^2+a^2)+℧ q r0-a Lz;

Vr[r_]:=Pr[r]^2-Δr[r](μ^2 r^2+(Lz-a ε)^2+Q);             (* effektives radiales Potential *)
Vτ=Pτ^2-Δ(μ^2 r[τ]^2+(Lz-a ε)^2+Q);
Vθ[θ_]:=Q-Cos[θ]^2(a^2 (μ^2-ε^2)+Lz^2 Sin[θ]^(-2)); (* effektives latitudinales Potential *)

т[τ_]:=Evaluate[t[τ]/.sol][[1]];                        (* Koordinatenzeit nach Eigenzeit *)
д[ξ_]:=Quiet[zt /.FindRoot[т[zt]-ξ, {zt, 0}]];          (* Eigenzeit nach Koordinatenzeit *)
T :=Quiet[д[tk]];                       
 
ю[τ_]:=Evaluate[t'[τ]/.sol][[1]];
γ[τ_]:=If[μ==0, "Infinity", ю[τ]];                                           (* totale ZD *)
R[τ_]:=Evaluate[r[τ]/.sol][[1]];                                (* Boyer-Lindquist Radius *)
Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]];                           
Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];
ß[τ_]:=Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ];

ς[τ_]:=Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0=Sqrt[χ/δ/Ξ];                         (* gravitative ZD *)
ω[τ_]:=(a(2R[τ]-℧^2))/Χi[τ]; ω0=(a(2r0-℧^2))/χ; ωζ=(a(2r[τ]-℧^2))/Χ;   (* F-Drag Winkelg *)
Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2];            (* Frame Dragging beobachtete Geschwindigkeit *)
й[τ_]:=ω[τ] яi[τ] ς[τ]; й0=ω0 Ы ς0;              (* Frame Dragging lokale Geschwindigkeit *)
 
ж[τ_]:=Sqrt[ς[τ]^2-1]/ς[τ]; ж0=Sqrt[ς0^2-1]/ς0;                  (* Fluchtgeschwindigkeit *)
V[τ_]:=If[μ==0, 1, Re[Sqrt[-ς[τ]^2+ю[τ]^2]/ю[τ]]];
                                                  (* Fluchtgeschwindigkeit von r0 nach r1 *)
vd1:=v1/.NSolve[Sqrt[δ Ξ/χ]/Sqrt[1-v1^2]+((q r0 ℧)/(r0^2+a^2 Cos[θ0]^2))==Sqrt[((a^2+
(-2+r1) r1+℧^2) (r1^2+a^2 Cos[θ0]^2))/((a^2+r1^2)^2-a^2 (a^2+(-2+r1) r1+℧^2) Sin[θ0]^2)]+
(a^2 q r1 ℧ (2 r1-℧^2) Sin[θ0]^2)/((r1^2+a^2 Cos[θ0]^2) ((a^2+r1^2)^2-a^2 (a^2+(-2+r1) r1+
℧^2) Sin[θ0]^2))+((q r1 ℧)/(r1^2+a^2 Cos[θ0]^2))&&v1>0,v1][[1]];
                                                          (* lokale Dreiergeschwindigkeit *)
vd[τ_]:=Abs[(Sqrt[Δ Σ^3 Χ-ε^2 Σ^2 Χ^2-2 a Lz ε Σ^2 Χ ℧^2-a^2 Lz^2 Σ^2 ℧^4+
4 a Lz ε Σ^2 Χ r[τ]+2 q ε Σ Χ^2 ℧ r[τ]+4 a^2 Lz^2 Σ^2 ℧^2 r[τ]+2 a Lz q Σ Χ ℧^3 r[τ]-
4 a^2 Lz^2 Σ^2 r[τ]^2-4 a Lz q Σ Χ ℧ r[τ]^2-q^2 Χ^2 ℧^2 r[τ]^2])/(ε Σ Χ+
a Lz Σ ℧^2-2 a Lz Σ r[τ]-q Χ ℧ r[τ])];

v[τ_]:=If[μ==0, 1, Evaluate[vlt'[τ]/.sol][[1]]];

vnt[τ_]:=Evaluate[Sqrt[(φ'[τ] r[τ]/Csc[θ[τ]])^2+(θ'[τ] r[τ])^2+r'[τ]^2]/.sol][[1]]
ν:=If[μ==0, 1, Re[Sqrt[(Δ Σ-Χ(εζ-Lζ ωζ)^2)/(μ Χ (εζ-Lζ ωζ)^2)]]];
vesc[τ_]:=Abs[(\[Sqrt](R[τ]^2 (R[τ]^2 (Δi[τ] Σi[τ]-Χi[τ])+2 q R[τ] Χi[τ] ℧-q^2 Χi[τ] ℧^2)+
2 a^2 R[τ] (R[τ] Δi[τ] Σi[τ]-R[τ] Χi[τ]+q Χi[τ] ℧) Cos[Θ[τ]]^2+a^4 (Δi[τ] Σi[τ]-
Χi[τ]) Cos[Θ[τ]]^4))/(Sqrt[Χi[τ]] (R[τ] (R[τ]-q ℧)+a^2 Cos[Θ[τ]]^2))];

dst[τ_]:=Evaluate[str[τ]/.sol][[1]];                                           (* Strecke *)
 
pΘ[τ_]:=Evaluate[Ξ θ'[τ] /. sol][[1]];
pR[τ_]:=Evaluate[r'[τ] Ξ/δ /. sol][[1]];

epot[τ_]:=ε+μ-ekin[τ];                                             (* potentielle Energie *)
ekin[τ_]:=If[μ==0, ς[τ], 1/Sqrt[1-v[τ]^2]-1];                       (* kinetische Energie *)

drwf=vr0/(Sqrt[(2+r0)/r0] Sqrt[1-v0^2 μ^2]);
duwf=vθ0/(Sqrt[r0^2] Sqrt[1-v0^2 μ^2]);
dfwf=(vφ0 Csc[θ0]^4 (r0^2 Sin[θ0]^2)^(3/2))/(r0^4 Sqrt[1-v0^2 μ^2]);
dtwf=1/(-2+r0) (-2 a Sin[θ0]^2 dfwf+r0 Sqrt[-μ+(2 μ)/r0+(1-4/r0^2) drwf^2+(-2+r0) r0 duwf^2-
2 r0 Sin[θ0]^2 dfwf^2+r0^2 Sin[θ0]^2 dfwf^2+(4 a^2 Sin[θ0]^4 dfwf^2)/r0^2]);
 
                                                               (* beobachtete Inklination *)
ink0:=б/. Solve[Z'[0]/ю[0] Cos[б]==-Y'[0]/ю[0] Sin[б]&&б>0&&б<2π&&б<δp[r0, a], б][[1]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 7) DIFFERENTIALGLEICHUNG |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
dp= Style[\!\(\*SuperscriptBox[\(Y\),\(Y\)]\), White]; n0[z_] := Chop[Re[N[Simplify[z]]]];

initcon = NSolve[
dr0 == pr0 δ/Ξ
&&
dθ0 == pθ0/Ξ
&&
dφ0 == 1/(δ Ξ Sin[θ0]^2) (ε (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+Lz (δ-a^2 Sin[θ0]^2)-
q ℧ r0 a Sin[θ0]^2)
&&
-μ  == -(((r0^2+a^2 Cos[θ0]^2) dr0^2)/(a^2-2 r0+r0^2+℧^2))+((a^2-2 r0+
r0^2+℧^2-a^2 Sin[θ0]^2) (dt0)^2)/(r0^2+a^2 Cos[θ0]^2)+(-r0^2-
a^2 Cos[θ0]^2) dθ0^2+(2 a (2 r0-℧^2) Sin[θ0]^2 dt0 dφ0)/(r0^2+
a^2 Cos[θ0]^2)+((-(a^2+r0^2)^2 Sin[θ0]^2+a^2 (a^2-2 r0+r0^2+
℧^2) Sin[θ0]^4) dφ0^2)/(r0^2+a^2 Cos[θ0]^2)
&&
dt0 > 0,
{dθ0, dr0, dt0, dφ0}, Reals];

initkon = NSolve[
dr0 == pr0 δ/Ξ
&&
dθ0 == pθ0/Ξ
&&
dφ0 == 1/(δ Ξ Sin[θ0]^2) (ε (-δ Xj+a Sin[θ0]^2 (r0^2+a^2))+Lz (δ-a^2 Sin[θ0]^2)-
q ℧ r0 a Sin[θ0]^2)
&&
dt0 == ς0/Sqrt[1-μ^2 v0^2]
&&
dt0 > 0,
{dθ0, dr0, dt0, dφ0}, Reals];
                                       
DG1={
 
t''[τ]==-(((r'[τ] ((a^2+r[τ]^2) (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+
2 (-℧^2+r[τ]) t'[τ]))+a (2 a^4 Cos[θ[τ]]^2+a^2 ℧^2 (3+Cos[2 θ[τ]]) r[τ]-
a^2 (3+Cos[2 θ[τ]]) r[τ]^2+4 ℧^2 r[τ]^3-6 r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))/(a^2+℧^2+(-2+
r[τ]) r[τ])+a^2 θ'[τ] (Sin[2 θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])-2 a Cos[θ[τ]] (℧^2-
2 r[τ]) Sin[θ[τ]]^3 φ'[τ]))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^2),

t'[0]==If[μ==0, dt0/.initcon[[1]], ς0/Sqrt[1-v0^2]],
t[0]==0,
 
r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
(a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
 
r'[0]==dr0/.initcon[[1]],
r[0]==r0,
 
θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
10 a^4 r[τ]+11 a^4 r[τ]^2-8 a^2 ℧^2 r[τ]^2+16 a^2 r[τ]^3+16 a^2 r[τ]^4+8 r[τ]^6+
a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
 
θ'[0]==dθ0/.initcon[[1]],
θ[0]==θ0,
 
φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),
 
φ'[0]==dφ0/.initcon[[1]],
φ[0]==φ0,
 
str'[τ]==If[μ==0, 1, vd[τ]/Abs[Sqrt[1-vd[τ]^2]]],
str[0]==0,
vlt'[τ]==If[μ==0, 1, vd[τ]],
vlt[0]==0
 
};

DG2={
 
t'[τ]==1/(Δ Σ Sin[θ[τ]]^2) (Lz (Δ XJ-a Sin[θ[τ]]^2 (r[τ]^2+a^2))+ε (-Δ XJ^2+
Sin[θ[τ]]^2 (r[τ]^2+a^2)^2)-q ℧ r[τ] Sin[θ[τ]]^2 (r[τ]^2+a^2)),
t[0]==0,
 
r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
(a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
 
r'[0]==(pr0 δ)/Ξ,
r[0]==r0,
 
θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
10 a^4 r[τ]+11 a^4 r[τ]^2-8 a^2 ℧^2 r[τ]^2+16 a^2 r[τ]^3+16 a^2 r[τ]^4+8 r[τ]^6+
a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
 
θ'[0]==pθ0/Ξ,
θ[0]==θ0,
 
φ'[τ]==1/(Δ Σ Sin[θ[τ]]^2) (ε (-Δ XJ+a Sin[θ[τ]]^2 (r[τ]^2+a^2))+Lz (Δ-a^2 Sin[θ[τ]]^2)-
q ℧ r[τ] a Sin[θ[τ]]^2),
φ[0]==φ0,
 
str'[τ]==If[μ==0, 1, vd[τ]/Abs[Sqrt[1-vd[τ]^2]]],
str[0]==0,
vlt'[τ]==If[μ==0, 1, vd[τ]],
vlt[0]==0
 
};

DG3={
 
t'[τ]==1/(Δ Σ Sin[θ[τ]]^2) (Lz (Δ XJ-a Sin[θ[τ]]^2 (r[τ]^2+a^2))+ε (-Δ XJ^2+
Sin[θ[τ]]^2 (r[τ]^2+a^2)^2)-q ℧ r[τ] Sin[θ[τ]]^2 (r[τ]^2+a^2)),
t[0]==0,
 
r'[τ]==Sign[vr0] Sqrt[Vτ]/Σ,
r[0]==r0,

θ'[τ]==Sign[vθ0] Sqrt[Q-Cos[θ[τ]]^2 (a^2 (μ^2-ε^2)+Lz^2/Sin[θ[τ]]^2)]/Σ,
θ[0]==θ0,

φ'[τ]==1/(Δ Σ Sin[θ[τ]]^2) (ε (-Δ XJ+a Sin[θ[τ]]^2 (r[τ]^2+a^2))+Lz (Δ-a^2 Sin[θ[τ]]^2)-
q ℧ r[τ] a Sin[θ[τ]]^2),
φ[0]==φ0,
 
str'[τ]==If[μ==0, 1, vd[τ]/Abs[Sqrt[1-vd[τ]^2]]],
str[0]==0,
vlt'[τ]==If[μ==0, 1, vd[τ]],
vlt[0]==0
 
};

DG4={
 
t''[τ]==((-4 a^2 r[τ] Sin[2 θ[τ]] t'[τ] θ'[τ]+r'[τ] (4 a^2 Sin[θ[τ]]^2 t'[τ]+
r[τ]^3 (q ℧-2 t'[τ]+6 a Sin[θ[τ]]^2 φ'[τ])))/((-2+r[τ]) r[τ]^4+4 a^2 r[τ] Sin[θ[τ]]^2)),
t'[0]==dtwf,
t[0]==0,

r''[τ]==((r'[τ]^2-t'[τ]^2+t'[τ] (q ℧+2 a Sin[θ[τ]]^2 φ'[τ])+r[τ]^3 (θ'[τ]^2+
Sin[θ[τ]]^2 φ'[τ]^2))/(r[τ] (2+r[τ]))),
r'[0]==drwf,
r[0]==r0,

θ''[τ]==((-2 r[τ]^2 r'[τ] θ'[τ]+Cos[θ[τ]] Sin[θ[τ]] φ'[τ] (-4 a t'[τ]+
r[τ]^3 φ'[τ]))/(r[τ]^3)),
θ'[0]==duwf,
θ[0]==θ0,

φ''[τ]==-((2 (r'[τ] (-a q ℧+a r[τ] t'[τ]+((-2+r[τ]) r[τ]^3-2 a^2 Sin[θ[τ]]^2) φ'[τ])+
r[τ] θ'[τ] (-2 a Cot[θ[τ]] (-2+r[τ]) t'[τ]+(Cot[θ[τ]] (-2+r[τ]) r[τ]^3+
2 a^2 Sin[2 θ[τ]]) φ'[τ])))/((-2+r[τ]) r[τ]^4+4 a^2 r[τ] Sin[θ[τ]]^2)),
φ'[0]==dfwf,
φ[0]==φ0,
 
str'[τ]==vd[τ],
str[0]==0,
vlt'[τ]==vd[τ],
vlt[0]==0
 
};

DG5={
 
t''[τ]==0,
 
t'[0]==1,
t[0]==0,
 
r''[τ]==-(1-℧ q)/r[τ]^2+(Sqrt[vφ0^2+vθ0^2] r0)^2/r[τ]^3,
 
r'[0]==vr0,
r[0]==r0,
 
θ''[τ]==-2 θ'[τ] r'[τ]/r[τ]+Sin[θ[τ]] Cos[θ[τ]] φ'[τ]^2,
 
θ'[0]==vθ0/r0,
θ[0]==θ0,
 
φ''[τ]==-2 φ'[τ] (r'[τ]+r[τ] θ'[τ] Cot[θ[τ]])/r[τ],
 
φ'[0]==vφ0/r0 Csc[θ0],
φ[0]==φ0,
 
vlt'[τ]==Sqrt[r'[τ]^2+φ'[τ]^2 r[τ]^2],
vlt[0]==0,
str'[τ]==Sqrt[r'[τ]^2+φ'[τ]^2 r[τ]^2],
str[0]==0
 
};

DGL = If[dgl==1, DG1, If[dgl==2, DG2, If[dgl==3, DG3, If[dgl==4, DG4, DG5]]]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 8) INTEGRATION |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
sol=NDSolve[DGL, {t, r, θ, φ, vlt, str}, {τ, 0, tmax+1/1000},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
Method-> mta,
InterpolationOrder-> All,
StepMonitor :> (laststep=plunge; plunge=τ;
stepsize=plunge-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<1*^-4, 0, 1])}];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 9) KOORDINATEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
X[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];            (* kartesisch *)
Y[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
Z[τ_]:=Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];
 
x[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];       (* Plotkoordinaten *)
y[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
z[τ_]:=Z[τ];
 
XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2];  (* kartesischer Radius *)
 
Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};      (* Rotationsmatrix *)
xYz[{x_, y_, z_}, β_]:={x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
xyZ[{x_, y_, z_}, ψ_]:={x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 10) PLOT EINSTELLUNGEN |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
PR=r1;                                                                      (* Plot Range *)
VP={r0, r0, r0};                                                     (* Perspektive x,y,z *)
d1=10;                                                                    (* Schweiflänge *)
plp=50;                                                            (* Flächenplot Details *)
Plp=Automatic;                                                          (* Kurven Details *)
 
w1l=0; w2l=0; w1r=0; w2r=0;                                          (* Startperspektiven *)
Mrec=100; mrec=10;                                       (* Parametric Plot Subdivisionen *)
imgsize=380;                                                                 (* Bildgröße *)
 
s[text_]:=Style[text, FontFamily->"Consolas", FontSize->11];               (* Anzeigestil *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 11) PLOT NACH EIGENZEIT ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Plot[R[tt], {tt, 0, plunge},
Frame->True, PlotStyle->Red, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, plunge}, All}, GridLines->{{}, {rA, rI}},
PlotLabel -> "r(τ)"]

Plot[Mod[180/Pi Θ[tt], 360], {tt, 0, plunge},
Frame->True, PlotStyle->Cyan, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, plunge}, {0, 360}}, GridLines->{{}, {90, 180, 270}},
PlotLabel -> "θ(τ)"]

Plot[Mod[180/Pi Φ[tt], 360], {tt, 0, plunge},
Frame->True, PlotStyle->Magenta, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, plunge}, {0, 360}}, GridLines->{{}, {90, 180, 270}},
PlotLabel -> "φ(τ)"]

Plot[v[tt], {tt, 0, plunge},
Frame->True, PlotStyle->Orange, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, plunge}, All}, GridLines->{{}, {0, 1}},
PlotLabel -> "v(τ)"]
 
displayP[T_]:=Grid[{
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[tp]], s["GM/c³"], s[dp]},
{s[" t coord"], " = ", s[n0[т[tp]]], s["GM/c³"], s[dp]},
{s[" ṫ total"], " = ", s[n0[ю[tp]]], s["dt/dτ"], s[dp]},
{s[" ς gravt"], " = ", s[n0[If[dgl==5, 1, ς[tp]]]], s["dt/dτ"], s[dp]},
{s[" γ kinet"], " = ", If[μ==0, s[n0[ς[tp]]], s[n0[If[dgl==5,
vnt[tp]^2/2, 1/Sqrt[1-v[tp]^2]]]]], s["dt/dτ"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[tp]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[X[tp]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[tp]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[tp]]], s["GM/c²"], s[dp]},
{s[" s dstnc"], " = ", s[n0[dst[tp]]], s["GM/c²"], s[dp]},

{s[" r coord"], " = ", s[n0[R[tp]]], s["GM/c²"], s[dp]},
{s[" φ longd"], " = ", s[n0[Φ[tp] 180/π]], s["deg"], s[dp]},
{s[" θ lattd"], " = ", s[n0[Θ[tp] 180/π]], s["deg"], s[dp]},
{s[" d¹r/dτ¹"], " = ", s[n0[R'[tp]]], s["c"], s[dp]},
{s[" d¹φ/dτ¹"], " = ", s[n0[Φ'[tp]]], s["c\.b3/G/M"], s[dp]},
{s[" d¹θ/dτ¹"], " = ", s[n0[Θ'[tp]]], s["c\.b3/G/M"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ", s[n0[R''[tp]]], s["c⁴/G/M"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Φ''[tp]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Θ''[tp]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},

{s[" ℧ cntrl"], " = ", s[n0[℧]], s["Q/M"], s[dp]},
{s[" q prtcl"], " = ", s[n0[q]], s["q/m"], s[dp]},
{s[" M irred"], " = ", s[N[mirr]], s["M"], s[dp]},
{s[" E kinet"], " = ", s[n0[If[dgl==5, vnt[tp]^2/2, ekin[tp]]]], s["mc²"], s[dp]},
{s[" E poten"], " = ", s[n0[If[dgl==5, -(1-℧ q)/R[tp], epot[tp]]]], s["mc²"], s[dp]},
{s[" E total"], " = ", s[n0[If[dgl==5, v0^2/2-(1-℧ q)/r0+1, ε]]], s["mc²"], s[dp]},
{s[" CarterQ"], " = ", s[n0[Qk]], s["GMm/c"], s[dp]},
{s[" L axial"], " = ", s[n0[If[dgl==5, vφ0 r0, Lz]]], s["GMm/c"], s[dp]},
{s[" L polar"], " = ", s[n0[If[dgl==5, vθ0 r0, pΘ[tp]]]], s["GMm/c"], s[dp]},
{s[" α dv/dτ"], " = ", s[n0[v'[tp]]], s["c⁴/G/M"], s[dp]},

{s[" ω fdrag"], " = ", s[n0[Abs[ω[tp]]]], s["c³/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[Abs[й[tp]]]], s["c"], s[dp]},
{s[" Ω fdrag"], " = ", s[n0[Abs[Ω[tp]]]], s["c"], s[dp]},
{s[" v propr"], " = ", s[n0[If[dgl==5, vnt[tp], v[tp]/Sqrt[1-v[tp]^2]]]], s["c"], s[dp]},
{s[" v escpe"], " = ", s[n0[vesc[tp]]], s["c"], s[dp]},
{s[" v obsvd"], " = ", s[n0[ß[tp]]], s["c"], s[dp]},
{s[" v r,loc"], " = ", s[n0[If[dgl==5, R'[tp], vr[tp]]]], s["c"], s[dp]},
{s[" v θ,loc"], " = ", s[n0[If[dgl==5, Θ'[tp] R[tp], vθ[tp]]]], s["c"], s[dp]},
{s[" v φ,loc"], " = ", s[n0[If[dgl==5, Φ'[tp] R[tp]/Csc[Θ[tp]], vφ[tp]]]], s["c"], s[dp]},
{s[" v local"], " = ", s[n0[If[dgl==5, vnt[tp], v[tp]]]], s["c"], s[dp]},
{s[" "], s[" "], s["                       "], s["         "]}},
Alignment-> Left, Spacings-> {0, 0}];
 
plot1b[{xx_, yy_, zz_, tk_, w1_, w2_}]:=                                     (* Animation *)
Show[

Graphics3D[{
{PointSize[0.011], Red, Point[
Xyz[xyZ[{x[tp], y[tp], z[tp]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> {
{-(2 Sign[Abs[xx]]+1) PR, +(2 Sign[Abs[xx]]+1) PR},
{-(2 Sign[Abs[yy]]+1) PR, +(2 Sign[Abs[yy]]+1) PR},
{-(2 Sign[Abs[zz]]+1) PR, +(2 Sign[Abs[zz]]+1) PR}
},
SphericalRegion->False,
ImagePadding-> 1],

horizons[A, None, w1, w2],

If[a==0, {}, ParametricPlot3D[
Xyz[xyZ[{
Sin[prm] a,
Cos[prm] a,
0}, w1], w2],
{prm, 0, 2π},
PlotStyle -> {Thickness[0.005], Orange}]],

If[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]],

If[tk==0, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2],
{tt, Max[0, д[т[tp]-1/2 π/ω0]], tp},
PlotStyle -> {Thickness[0.001], Dashed, Purple},
PlotPoints-> Automatic,
MaxRecursion-> 12]]],

If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp},
PlotStyle-> {Thickness[0.004]},
ColorFunction-> Function[{x, y, z, t},
Hue[0, 1, 0.5, If[tp<0,
Max[Min[(+tp+(-t+d1))/d1, 1], 0], Max[Min[(-tp+(t+d1))/d1, 1], 0]]]],
ColorFunctionScaling-> False,
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],

If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2],
{tt, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]},
PlotStyle-> {Thickness[0.004], Opacity[0.6], Darker[Gray]},
PlotPoints-> Plp,
MaxRecursion-> mrec]]],

ViewPoint-> {xx, yy, zz}];
 
Do[
Print[Rasterize[Grid[{{
plot1b[{0, -Infinity, 0, tp, w1l, w2l}],
plot1b[{0, 0, +Infinity, tp, w1r, w2r}],
displayP[tp]
}, {"  ", "  ", "                                             "}
}, Alignment->Left]]],
{tp, 0, tMax, tMax/1}]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 12) PLOT NACH KOORDINATENZEIT ||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Plot[R[д[tt]], {tt, 0, TMax},
Frame->True, PlotStyle->Red, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, TMax}, All}, GridLines->{{}, {rA, rI}},
PlotLabel -> "r(t)"]

Plot[Mod[180/Pi Θ[д[tt]], 360], {tt, 0, TMax},
Frame->True, PlotStyle->Cyan, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, TMax}, {0, 360}}, GridLines->{{}, {90, 180, 270}},
PlotLabel -> "θ(t)"]

Plot[Mod[180/Pi Φ[д[tt]], 360], {tt, 0, TMax},
Frame->True, PlotStyle->Magenta, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, TMax}, {0, 360}}, GridLines->{{}, {90, 180, 270}},
PlotLabel -> "φ(t)"]

Plot[v[д[tt]], {tt, 0, TMax},
Frame->True, PlotStyle->Orange, AspectRatio->1/5, ImagePadding->{{40, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, TMax}, All}, GridLines->{{}, {0, 1}},
PlotLabel -> "v(t)"]
 
displayC[T_]:=Grid[{
{s[" t coord"], " = ", s[n0[tk]], s["GM/c³"], s[dp]},
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c³"], s[dp]},
{s[" ṫ total"], " = ", s[n0[ю[T]]], s["dt/dτ"], s[dp]},
{s[" ς gravt"], " = ", s[n0[If[dgl==5, 1, ς[T]]]], s["dt/dτ"], s[dp]},
{s[" γ kinet"], " = ", If[μ==0, s[n0[ς[T]]], s[n0[If[dgl==5,
vnt[T]^2/2, 1/Sqrt[1-v[T]^2]]]]], s["dt/dτ"], s[dp]},
{s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c²"], s[dp]},
{s[" x carts"], " = ", s[n0[X[T]]], s["GM/c²"], s[dp]},
{s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c²"], s[dp]},
{s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c²"], s[dp]},
{s[" s dstnc"], " = ", s[n0[dst[T]]], s["GM/c²"], s[dp]},

{s[" r coord"], " = ", s[n0[R[T]]], s["GM/c²"], s[dp]},
{s[" φ longd"], " = ", s[n0[Φ[T] 180/π]], s["deg"], s[dp]},
{s[" θ lattd"], " = ", s[n0[Θ[T] 180/π]], s["deg"], s[dp]},
{s[" d¹r/dτ¹"], " = ", s[n0[R'[T]]], s["c"], s[dp]},
{s[" d¹φ/dτ¹"], " = ", s[n0[Φ'[T]]], s["c\.b3/G/M"], s[dp]},
{s[" d¹θ/dτ¹"], " = ", s[n0[Θ'[T]]], s["c\.b3/G/M"], s[dp]},
{s[" d\.b2r/dτ\.b2"], " = ", s[n0[R''[T]]], s["c⁴/G/M"], s[dp]},
{s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Φ''[T]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Θ''[T]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},

{s[" ℧ cntrl"], " = ", s[n0[℧]], s["Q/M"], s[dp]},
{s[" q prtcl"], " = ", s[n0[q]], s["q/m"], s[dp]},
{s[" M irred"], " = ", s[N[mirr]], s["M"], s[dp]},
{s[" E kinet"], " = ", s[n0[If[dgl==5, vnt[T]^2/2, ekin[T]]]], s["mc²"], s[dp]},
{s[" E poten"], " = ", s[n0[If[dgl==5, -(1-℧ q)/R[T], epot[T]]]], s["mc²"], s[dp]},
{s[" E total"], " = ", s[n0[If[dgl==5, v0^2/2-(1-℧ q)/r0+1, ε]]], s["mc²"], s[dp]},
{s[" CarterQ"], " = ", s[n0[Qk]], s["GMm/c"], s[dp]},
{s[" L axial"], " = ", s[n0[If[dgl==5, vφ0 r0, Lz]]], s["GMm/c"], s[dp]},
{s[" L polar"], " = ", s[n0[If[dgl==5, vθ0 r0, pΘ[T]]]], s["GMm/c"], s[dp]},
{s[" α dv/dτ"], " = ", s[n0[v'[T]]], s["c⁴/G/M"], s[dp]},

{s[" ω fdrag"], " = ", s[n0[Abs[ω[T]]]], s["c³/G/M"], s[dp]},
{s[" v fdrag"], " = ", s[n0[Abs[й[T]]]], s["c"], s[dp]},
{s[" Ω fdrag"], " = ", s[n0[Abs[Ω[T]]]], s["c"], s[dp]},
{s[" v propr"], " = ", s[n0[If[dgl==5, vnt[T], v[T]/Sqrt[1-v[T]^2]]]], s["c"], s[dp]},
{s[" v escpe"], " = ", s[n0[vesc[T]]], s["c"], s[dp]},
{s[" v obsvd"], " = ", s[n0[ß[T]]], s["c"], s[dp]},
{s[" v r,loc"], " = ", s[n0[If[dgl==5, R'[T], vr[T]]]], s["c"], s[dp]},
{s[" v θ,loc"], " = ", s[n0[If[dgl==5, Θ'[T] R[T], vθ[T]]]], s["c"], s[dp]},
{s[" v φ,loc"], " = ", s[n0[If[dgl==5, Φ'[T] R[T]/Csc[Θ[T]], vφ[T]]]], s["c"], s[dp]},
{s[" v local"], " = ", s[n0[If[dgl==5, vnt[T], v[T]]]], s["c"], s[dp]},
{s[" "], s[" "], s["                       "], s["         "]}},
Alignment-> Left, Spacings-> {0, 0}];
 
plot1a[{xx_, yy_, zz_, tk_, w1_, w2_}]:=                                     (* Animation *)
Show[

Graphics3D[{
{PointSize[0.011], Red, Point[
Xyz[xyZ[{x[T], y[T], z[T]}, w1], w2]]}},
ImageSize-> imgsize,
PlotRange-> {
{-(2 Sign[Abs[xx]]+1) PR, +(2 Sign[Abs[xx]]+1) PR},
{-(2 Sign[Abs[yy]]+1) PR, +(2 Sign[Abs[yy]]+1) PR},
{-(2 Sign[Abs[zz]]+1) PR, +(2 Sign[Abs[zz]]+1) PR}
},
SphericalRegion->False,
ImagePadding-> 1],

horizons[A, None, w1, w2],

If[a==0, {}, ParametricPlot3D[
Xyz[xyZ[{
Sin[prm] a,
Cos[prm] a,
0}, w1], w2],
{prm, 0, 2π},
PlotStyle -> {Thickness[0.005], Orange}]],

If[a==0, {},
Graphics3D[{{PointSize[0.009], Purple, Point[
Xyz[xyZ[{
Sin[-φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2]]}}]],

If[tk==0, {}, If[a==0, {},
ParametricPlot3D[
Xyz[xyZ[{
Sin[-φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
Cos[-φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
z0[A]}, w1], w2],
{tt, Max[0, tk-1/2 π/ω0], tk},
PlotStyle -> {Thickness[0.001], Dashed, Purple},
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],

Block[{$RecursionLimit = Mrec},
If[tk==0, {},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[TMax<0, Min[0, T+d1], Max[0, T-d1]], T},
PlotStyle-> {Thickness[0.004]},
ColorFunction-> Function[{x, y, z, t},
Hue[0, 1, 0.5, If[TMax<0, Max[Min[(+T+(-t+d1))/d1, 1], 0], Max[Min[(-T+(t+d1))/d1, 1], 0]]]],
ColorFunctionScaling-> False,
PlotPoints-> Automatic,
MaxRecursion-> mrec]]],

If[tk==0, {},
Block[{$RecursionLimit = Mrec},
ParametricPlot3D[
Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2],
{tt, 0, If[Tmax<0, Min[-1*^-16, T+d1/3], Max[1*^-16, T-d1/3]]},
PlotStyle-> {Thickness[0.004], Opacity[0.6], Darker[Gray]},
PlotPoints-> Plp,
MaxRecursion-> mrec]]],

ViewPoint-> {xx, yy, zz}];
 
Quiet[Do[
Print[Rasterize[Grid[{{
plot1a[{0, -Infinity, 0, tk, w1l, w2l}],
plot1a[{0, 0, Infinity, tk, w1r, w2r}],
displayC[Quiet[д[tk]]]
}, {"  ", "  ", "                                             "}
}, Alignment->Left]]],
{tk, 0, TMax, TMax/1}]]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 13) EXPORTOPTIONEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* Export als HTML Dokument *)
(* Export["Y:\\export\\dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"]  *)
(* Export direkt als Bildsequenz *)
(* Do[Export["Y:\\export\\dateiname" <> ToString[tk] <> ".png", Rasterize[...]            *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||| http://kerr.newman.yukerez.net ||||| Simon Tyran, Vienna |||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)