Seite 1 von 1

Kerr Newman De Sitter Metrik (KNdS)

Verfasst: Mo 21. Aug 2023, 00:42
von Yukterez
Bild
Bild Das ist die deutschsprachige Version.   Bild English versions will soon be available on en.yukterez.net and yukipedia. Bild
Bild
Verwandte Beiträge: Kerr Newman Metrik || Schwarzschild De Sitter Metrik || Einstein Maxwell Feldgleichungen || Relativistischer Raytracer Bild
Bild

Bild
Bild
Vorarbeit: Bild
Bild

Getrennte Darstellung der Horizonte & Ergosphären (in dem Fall Ergokuppeln bzw. Ergodome), wie oben für M=1, a=9/10, ℧=2/5, Λ=1/9, Θ=85°

Bild

Gallerie für verschiedene positive und negative Λ, wobei die Parameter des schwarzen Lochs M=1, a=9/10, ℧=2/5, Θ=85°gleich bleiben:

Bild

Der Simulator für die Bahnen von Photonen, geladenen und neutralen Partikeln in Boyer Lindquist Koordinaten funktioniert schon, wird in den nächsten Tagen aber noch auf all die Funktionen die auch der normale Kerr Newman Simulator hat upgegradet und auf die anderen Koordinaten erweitert. Diese Version hat nur die wichtigsten Funktionen und muss noch auf Herz und Nieren getestet werden, schönere Bilder mit interessanteren Bahnen als da oben folgen demnächst:

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Mathematica | knds.yukterez.net | Kerr Newman De Sitter Simulator Beta Release 1 || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
wp=MachinePrecision;
mt1=Automatic;
mt2={"EquationSimplification"-> "Residual"};
mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
mt4={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mta=mt2;                                                     (* mt1: Speed, mt3: Accuracy *)
dgl=1;
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
A=a;                                   (* pseudosphärisch [BL]: A=0, kartesisch [KS]: A=a *)
 
tmax=13.162954624632034;                                                     (* Eigenzeit *)
tMax=Min[tmax, plunge-1/100];                                         (* Integrationsende *)
 
r0 = 20/10;                                                                (* Startradius *)
r1 = 40/10;                                                                  (* Endradius *)
θ0 = π/2;                                                                  (* Breitengrad *)
φ0 = 0;                                                                     (* Längengrad *)
M  = 1;                                                                          (* Masse *)
a  = 9/10;                                                               (* Spinparameter *)
℧  = 2/5;                                       (* spezifische Ladung des schwarzen Lochs *)
Λ  = 1/9;                                                      (* kosmologische Konstante *)
q  = 0;                                             (* spezifische Ladung des Testkörpers *)
 
v0 = 64307651/100000000;                                        (* Anfangsgeschwindigkeit *)
α0 = 0;                                                      (* vertikaler Abschusswinkel *)
i0 = π/4;                                                       (* Bahninklinationswinkel *)
 
PR=r1;                                                                      (* Plot Range *)
VP={r0, r0, r0};                                                     (* Perspektive x,y,z *)
d1=1;                                                                     (* Schweiflänge *)
plp=40;                                                            (* Flächenplot Details *)
Plp=Automatic;                                                          (* Kurven Details *)
 
w1l=0; w2l=0; w1r=0; w2r=0;                                          (* Startperspektiven *)
Mrec=100; mrec=15;                                       (* Parametric Plot Subdivisionen *)
imgsize=380;                                                                 (* Bildgröße *)
 
s[text_]:=Style[text, FontFamily->"Consolas", FontSize->11];               (* Anzeigestil *)

μ=If[Abs[v0]==1, 0, If[Abs[v0]<1, -1, 1]];                   (* Baryon: μ=-1, Photon: μ=0 *)

vr0=v0 Sin[α0];                                     (* radiale Geschwindigkeitskomponente *)
vθ0=v0 Cos[α0] Sin[i0];                      (* longitudinale  Geschwindigkeitskomponente *)
vφ0=v0 Cos[α0] Cos[i0];                        (* latitudinale Geschwindigkeitskomponente *)

T[τ_]:=Evaluate[t[τ]/.sol][[1]];
R[τ_]:=Evaluate[r[τ]/.sol][[1]];
Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]];                           
Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];

Td[τ_]:=Evaluate[t'[τ]/.sol][[1]];
Rd[τ_]:=Evaluate[r'[τ]/.sol][[1]];
Φd[τ_]:=Evaluate[φ'[τ]/.sol][[1]];                           
Θd[τ_]:=Evaluate[θ'[τ]/.sol][[1]];

gTT[r_, θ_]:=-(((3+a^2 Λ)^2 (r^2+a^2 Cos[θ]^2) (3 (a^2+r^2)^2+a^2 (a^2+r^2)^2 Λ Cos[θ]^2+
a^2 (6 M r-3 r^2+r^4 Λ+a^2 (-3+r^2 Λ)-3 ℧^2) Sin[θ]^2))/(3 (6 M r-3 r^2+r^4 Λ+a^2 (-3+
r^2 Λ)-3 ℧^2) (3+a^2 Λ Cos[θ]^2) (a^2+r^2-a^2 Sin[θ]^2)^2));

gtt[r_, θ_]:=-((2 M r-(a^2+r^2) (1-(r^2 Λ)/3)-℧^2+a^2 (1+1/3 a^2 Λ Cos[θ]^2) Sin[θ]^2)/((1+
(a^2 Λ)/3)^2 (r^2+a^2 Cos[θ]^2)));

grr[r_, θ_]:=-((r^2+a^2 Cos[θ]^2)/(-2 M r+(a^2+r^2) (1-(r^2 Λ)/3)+℧^2));

gθθ[r_, θ_]:=-((r^2+a^2 Cos[θ]^2)/(1+1/3 a^2 Λ Cos[θ]^2));

gφφ[r_, θ_]:=-(((a^2+r^2)^2 (1+1/3 a^2 Λ Cos[θ]^2) Sin[θ]^2-a^2 (-2 M r+(a^2+r^2) (1-
(r^2 Λ)/3)+℧^2) Sin[θ]^4)/((1+(a^2 Λ)/3)^2 (r^2+a^2 Cos[θ]^2)));

gtφ[r_, θ_]:=-((a (-2 M r+(a^2+r^2) (1-(r^2 Λ)/3)+℧^2-(a^2+r^2) (1+
1/3 a^2 Λ Cos[θ]^2)) Sin[θ]^2)/((1+(a^2 Λ)/3)^2 (r^2+a^2 Cos[θ]^2)));

Aμ[r_, θ_]:={(3 r ℧)/((3+a^2 Λ) (r^2+a^2 Cos[θ]^2)),
0, 0, -((3 a r ℧ Sin[θ]^2)/((3+a^2 Λ) (r^2+a^2 Cos[θ]^2)))};

ς[τ_]:=Sqrt[-(((3+a^2 Λ)^2 (R[τ]^2+a^2 Cos[Θ[τ]]^2) (3 (a^2+R[τ]^2)^2+a^2 (a^2+
R[τ]^2)^2 Λ Cos[Θ[τ]]^2+a^2 (6 M R[τ]-3 R[τ]^2+R[τ]^4 Λ+a^2 (-3+R[τ]^2 Λ)-
3 ℧^2) Sin[Θ[τ]]^2))/(3 (6 M R[τ]-3 R[τ]^2+R[τ]^4 Λ+a^2 (-3+R[τ]^2 Λ)-3 ℧^2) (3+
a^2 Λ Cos[Θ[τ]]^2) (a^2+R[τ]^2-a^2 Sin[Θ[τ]]^2)^2))];

v[τ_]:= -(1/(2 Sqrt[3] μ))(\[Sqrt]((72 M (6+a^2 Λ+a^2 Λ Cos[2 Θ[τ]]) R[τ]^3 Td[τ]^2+
12 Λ (6+a^2 Λ+a^2 Λ Cos[2 Θ[τ]]) R[τ]^6 Td[τ]^2-6 a^2 M R[τ] (-8 (3+a^2 Λ)^2 Sin[Θ[τ]]^2-
12 Cos[Θ[τ]]^2 (6+a^2 Λ+a^2 Λ Cos[2 Θ[τ]]) Td[τ]^2)+a^2 (4 (3+a^2 Λ)^2 (3 a^2+a^4 Λ-3 ℧^2+
(3 a^2+a^4 Λ+3 ℧^2) Cos[2 Θ[τ]])-36 (a^2+℧^2) Cos[Θ[τ]]^2 (6+a^2 Λ+
a^2 Λ Cos[2 Θ[τ]]) Td[τ]^2)+R[τ]^2 (4 a^2 (3+a^2 Λ)^3 (3+Cos[2 Θ[τ]])+6 (6+a^2 Λ+
a^2 Λ Cos[2 Θ[τ]]) (-9 a^2+a^4 Λ-6 ℧^2+a^2 (-3+a^2 Λ) Cos[2 Θ[τ]]) Td[τ]^2)+R[τ]^4 (8 (3+
a^2 Λ)^3+3 (-72+24 a^2 Λ+7 a^4 Λ^2+8 a^4 Λ^2 Cos[2 Θ[τ]]+a^4 Λ^2 Cos[4 Θ[τ]]) Td[τ]^2))/((6+
a^2 Λ+a^2 Λ Cos[2 Θ[τ]]) (a^2 Cos[Θ[τ]]^2+R[τ]^2) (-3 (a^2+℧^2)+6 M R[τ]+(-3+a^2 Λ) R[τ]^2+
Λ R[τ]^4) Td[τ]^2)));

vrj[τ_]:=Sqrt[3-3 v[τ]^2 μ^2] Sqrt[-((R[τ]^2+a^2 Cos[Θ[τ]]^2)/(6 M R[τ]-3 R[τ]^2+R[τ]^4 Λ+
a^2 (-3+R[τ]^2 Λ)-3 ℧^2))] Rd[τ];

vθj[τ_]:=Sqrt[3-3 v[τ]^2 μ^2] Sqrt[(a^2 Cos[Θ[τ]]^2+R[τ]^2)/(3+a^2 Λ Cos[Θ[τ]]^2)] Θd[τ];

vφj[τ_]:=(Sqrt[3-3 v[τ]^2 μ^2] Sin[Θ[τ]]^2 (2 a q v[τ]^2 (3+a^2 Λ) μ^2 ℧ R[τ]-a (a^4 Λ-
6 ℧^2+12 M R[τ]+3 a^2 Λ R[τ]^2+2 Λ R[τ]^4+a^2 Λ Cos[2 Θ[τ]] (a^2+R[τ]^2)) Td[τ]+(a^6 Λ+
6 R[τ]^4+3 a^4 (1+Λ R[τ]^2)+a^2 (-3 ℧^2+6 M R[τ]+9 R[τ]^2+2 Λ R[τ]^4)+
a^2 Cos[2 Θ[τ]] (a^4 Λ+3 (℧^2-2 M R[τ]+R[τ]^2)+a^2 (3+Λ R[τ]^2))) Φd[τ]))/(2 (3+
a^2 Λ)^2 (a^2 Cos[Θ[τ]]^2+R[τ]^2) Sqrt[(Sin[Θ[τ]]^2 (3 (a^2+R[τ]^2)^2+
a^2 Λ Cos[Θ[τ]]^2 (a^2+R[τ]^2)^2+a^2 (-3 ℧^2+6 M R[τ]-3 R[τ]^2+Λ R[τ]^4+a^2 (-3+
Λ R[τ]^2)) Sin[Θ[τ]]^2))/((3+a^2 Λ)^2 (a^2 Cos[Θ[τ]]^2+R[τ]^2))]);

ω[τ_]:=-2 Sqrt[-((a (12 M R[τ]+a^4 Λ+3 a^2 R[τ]^2 Λ+2 R[τ]^4 Λ-6 ℧^2+
a^2 (a^2+R[τ]^2) Λ Cos[2 Θ[τ]]) Sin[Θ[τ]]^2)/(12 a^2-48 M R[τ]+24 R[τ]^2-a^4 Λ-
8 a^2 R[τ]^2 Λ-8 R[τ]^4 Λ+24 ℧^2+12 a^2 Cos[2 Θ[τ]]+a^4 Λ Cos[4 Θ[τ]]))];

ω0=-2 Sqrt[-((a (12 M r0+a^4 Λ+3 a^2 r0^2 Λ+2 r0^4 Λ-6 ℧^2+
a^2 (a^2+r0^2) Λ Cos[2 θ0]) Sin[θ0]^2)/(12 a^2-48 M r0+24 r0^2-a^4 Λ-
8 a^2 r0^2 Λ-8 r0^4 Λ+24 ℧^2+12 a^2 Cos[2 θ0]+a^4 Λ Cos[4 θ0]))];

ν[τ_]:=Sqrt[-((a^2 (6 M R[τ]+a^2 R[τ]^2 Λ+R[τ]^4 Λ-3 ℧^2+a^2 (a^2+
R[τ]^2) Λ Cos[Θ[τ]]^2)^2 Sin[Θ[τ]]^2)/((6 M R[τ]-3 R[τ]^2+R[τ]^4 Λ+a^2 (-3+R[τ]^2 Λ)-
3 ℧^2) (3+a^2 Λ Cos[Θ[τ]]^2) (a^2+R[τ]^2-a^2 Sin[Θ[τ]]^2)^2))];

dt=Sqrt[((3+a^2 Λ)^2 (r0^2+a^2 Cos[θ0]^2) (-(a^2+r0^2)^2 (3+a^2 Λ Cos[θ0]^2)-
a^2 (6 M r0+r0^4 Λ+a^2 (-3+r0^2 Λ)-3 (r0^2+℧^2)) Sin[θ0]^2))/((6 M r0+r0^4 Λ+a^2 (-3+
r0^2 Λ)-3 (r0^2+℧^2)) (3+a^2 Λ Cos[θ0]^2) (a^2+r0^2-
a^2 Sin[θ0]^2)^2)]/(Sqrt[3] Sqrt[1-v0^2 μ^2]);

initcon=Solve[
vr0==dr Sqrt[3-3 v0^2 μ^2] Sqrt[-((r0^2+a^2 Cos[θ0]^2)/(6 M r0-3 r0^2+
r0^4 Λ+a^2 (-3+r0^2 Λ)-3 ℧^2))]&&
vθ0==dθ Sqrt[3-3 v0^2 μ^2] Sqrt[(r0^2+a^2 Cos[θ0]^2)/(3+a^2 Λ Cos[θ0]^2)]&&
vφ0==(Sqrt[3-3 v0^2 μ^2] (2 a q r0 v0^2 (3+a^2 Λ) μ^2 ℧-a dt (12 M r0+a^4 Λ+
3 a^2 r0^2 Λ+2 r0^4 Λ-6 ℧^2+a^2 (a^2+r0^2) Λ Cos[2 θ0])+dφ (6 r0^4+a^6 Λ+3 a^4 (1+
r0^2 Λ)+a^2 (6 M r0+9 r0^2+2 r0^4 Λ-3 ℧^2)+a^2 (a^4 Λ+a^2 (3+r0^2 Λ)+3 (-2 M r0+r0^2+
℧^2)) Cos[2 θ0])) Sin[θ0]^2)/(2 (3+a^2 Λ)^2 (r0^2+
a^2 Cos[θ0]^2) Sqrt[(Sin[θ0]^2 (3 (a^2+r0^2)^2+a^2 (a^2+r0^2)^2 Λ Cos[θ0]^2+
a^2 (6 M r0-3 r0^2+r0^4 Λ+a^2 (-3+r0^2 Λ)-3 ℧^2) Sin[θ0]^2))/((3+a^2 Λ)^2 (r0^2+
a^2 Cos[θ0]^2))]),
{dr, dθ, dφ}];

dR=initcon[[1,1,2]];
dΘ=initcon[[1,2,2]];
dΦ=initcon[[1,3,2]];

e0=+gtt[r0, θ0] dt+gtφ[r0, θ0] dΦ+q Aμ[r0, θ0][[1]];
Lz=-gφφ[r0, θ0] dΦ-gtφ[r0, θ0] dt-q Aμ[r0, θ0][[4]];

pt[τ_]:=-gtt[R[τ], Θ[τ]] Td[τ]-gtφ[R[τ], Θ[τ]] Φd[τ]-q Aμ[R[τ], Θ[τ]][[1]];
pθ[τ_]:=-gθθ[R[τ], Θ[τ]] Θd[τ];
pr[τ_]:=-grr[R[τ], Θ[τ]] Rd[τ];
pφ[τ_]:=-gφφ[R[τ], Θ[τ]] Φd[τ]-gtφ[R[τ], Θ[τ]] Td[τ]-q Aμ[R[τ], Θ[τ]][[4]];
 
x0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0];
y0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
z0[A_]:=r0 Cos[θ0];

Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
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[ψ]};
 
rH1=1/(2 Sqrt[3]) (√(-2 a^2+6/Λ+(9+Λ (-42 a^2+a^4 Λ-36 ℧^2))/(Λ (486 M^2 Λ+(-3+
a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+
(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+℧^2))^2])^(1/3))+(486 M^2 Λ+(-3+
a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+
(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+℧^2))^2])^(1/3)/Λ)-√(-4 a^2+12/Λ+
(-9+Λ (42 a^2-a^4 Λ+36 ℧^2))/(Λ (486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+
1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+
a^2 Λ) (a^2+℧^2))^2])^(1/3))-(486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+
1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+
a^2 Λ) (a^2+℧^2))^2])^(1/3)/Λ-(36 Sqrt[3] M)/(Λ √(-2 a^2+6/Λ+(9+Λ (-42 a^2+a^4 Λ-
36 ℧^2))/(Λ (486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+
a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+
℧^2))^2])^(1/3))+(486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+
a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+
℧^2))^2])^(1/3)/Λ))));

RH1[A_, w1_, w2_]:=xyZ[{Sqrt[rH1^2+A^2] Sin[θ] Cos[φ],
Sqrt[rH1^2+A^2] Sin[θ] Sin[φ], rH1 Cos[θ]}, w1];

rH2=1/(2 Sqrt[3]) (√(-2 a^2+6/Λ+(9+Λ (-42 a^2+a^4 Λ-36 ℧^2))/(Λ (486 M^2 Λ+(-3+
a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+
(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+℧^2))^2])^(1/3))+(486 M^2 Λ+(-3+
a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+
(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+℧^2))^2])^(1/3)/Λ)+√(-4 a^2+12/Λ+
(-9+Λ (42 a^2-a^4 Λ+36 ℧^2))/(Λ (486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+
1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+
a^2 Λ) (a^2+℧^2))^2])^(1/3))-(486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+
1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+
a^2 Λ) (a^2+℧^2))^2])^(1/3)/Λ-(36 Sqrt[3] M)/(Λ √(-2 a^2+6/Λ+(9+Λ (-42 a^2+a^4 Λ-
36 ℧^2))/(Λ (486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+
a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+
℧^2))^2])^(1/3))+(486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+
a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+
℧^2))^2])^(1/3)/Λ))));

RH2[A_, w1_, w2_]:=xyZ[{Sqrt[rH2^2+A^2] Sin[θ] Cos[φ],
Sqrt[rH2^2+A^2] Sin[θ] Sin[φ], rH2 Cos[θ]}, w1];

rH3=1/(2 Sqrt[3]) (-√(-2 a^2+6/Λ+(9+Λ (-42 a^2+a^4 Λ-36 ℧^2))/(Λ (486 M^2 Λ+(-3+
a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+
(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+℧^2))^2])^(1/3))+(486 M^2 Λ+(-3+
a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+
(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+℧^2))^2])^(1/3)/Λ)+√(-4 a^2+12/Λ+
(-9+Λ (42 a^2-a^4 Λ+36 ℧^2))/(Λ (486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+
1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+
a^2 Λ) (a^2+℧^2))^2])^(1/3))-(486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+
1/2 Sqrt[-4 ((-3+a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+
a^2 Λ) (a^2+℧^2))^2])^(1/3)/Λ+(36 Sqrt[3] M)/(Λ √(-2 a^2+6/Λ+(9+Λ (-42 a^2+a^4 Λ-
36 ℧^2))/(Λ (486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+
a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+
℧^2))^2])^(1/3))+(486 M^2 Λ+(-3+a^2 Λ)^3+108 Λ (-3+a^2 Λ) (a^2+℧^2)+1/2 Sqrt[-4 ((-3+
a^2 Λ)^2-36 Λ (a^2+℧^2))^3+(972 M^2 Λ+2 (-3+a^2 Λ)^3+216 Λ (-3+a^2 Λ) (a^2+
℧^2))^2])^(1/3)/Λ))));

RH3[A_, w1_, w2_]:=xyZ[{Sqrt[rH3^2+A^2] Sin[θ] Cos[φ],
Sqrt[rH3^2+A^2] Sin[θ] Sin[φ], rH3 Cos[θ]}, w1];

rE1=1/2 √(-((2 a^2)/3)+2/Λ+(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+a^2 Λ)^3+288 Λ (-3+a^2 Λ) (-3 (a^2+
℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+2 Sqrt[2] √((-18+48 a^2 Λ-5 a^4 Λ^2+
72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-
36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-
(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+12 a^2 Sin[θ]^2+
a^4 Λ Sin[2 θ]^2)+Sqrt[(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+
12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2))^2]/(2 Sqrt[2]))^(1/3)))-1/2 √(-((4 a^2)/3)+
4/Λ+(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+
a^2 Λ Cos[θ]^2) Sin[θ]^2)+Sqrt[(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+
a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2]/(2 Sqrt[2]))^(1/3))-(-18+48 a^2 Λ-5 a^4 Λ^2+
72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+a^2 Λ)^3+
288 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+2 Sqrt[2] √((-18+
48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+
(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+
a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-(12 M)/(Λ √(-((2 a^2)/3)+2/Λ+(-18+48 a^2 Λ-
5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+
a^2 Λ)^3+288 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+
2 Sqrt[2] √((-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+
a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+
9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2)+
(1/(2 Sqrt[2]))(√((-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+
12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2))^2)))^(1/3)))));

RE1[A_, w1_, w2_]:=xyZ[{Sqrt[rE1^2+A^2] Sin[θ] Cos[φ],
Sqrt[rE1^2+A^2] Sin[θ] Sin[φ], rE1 Cos[θ]}, w1];

rE2=1/2 √(-((2 a^2)/3)+2/Λ+(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+a^2 Λ)^3+288 Λ (-3+a^2 Λ) (-3 (a^2+
℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+2 Sqrt[2] √((-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+
36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+
a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-
(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+12 a^2 Sin[θ]^2+
a^4 Λ Sin[2 θ]^2)+Sqrt[(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+
12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2))^2]/(2 Sqrt[2]))^(1/3)))+1/2 √(-((4 a^2)/3)+
4/Λ+(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+
a^2 Λ Cos[θ]^2) Sin[θ]^2)+Sqrt[(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+
a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2]/(2 Sqrt[2]))^(1/3))-(-18+48 a^2 Λ-5 a^4 Λ^2+
72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+a^2 Λ)^3+
288 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+2 Sqrt[2] √((-18+
48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+
(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+
a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-(12 M)/(Λ √(-((2 a^2)/3)+2/Λ+(-18+48 a^2 Λ-
5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+
a^2 Λ)^3+288 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+
2 Sqrt[2] √((-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+
a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+
9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2)+
(1/(2 Sqrt[2]))(√((-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+
12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2))^2)))^(1/3)))));

RE2[A_, w1_, w2_]:=xyZ[{Sqrt[rE2^2+A^2] Sin[θ] Cos[φ],
Sqrt[rE2^2+A^2] Sin[θ] Sin[φ], rE2 Cos[θ]}, w1];

rE3=-(1/2) √(-((2 a^2)/3)+2/Λ+(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+
36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+a^2 Λ)^3+288 Λ (-3+
a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+2 Sqrt[2] √((-18+
48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])^3+
8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+
a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+9 Λ (-3+
a^2 Λ) (-12 (a^2+℧^2)+12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2)+Sqrt[(-18+48 a^2 Λ-
5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+
a^2 Λ)^3-9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+12 a^2 Sin[θ]^2+
a^4 Λ Sin[2 θ]^2))^2]/(2 Sqrt[2]))^(1/3)))+1/2 √(-((4 a^2)/3)+4/Λ+
(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+
a^2 Λ Cos[θ]^2) Sin[θ]^2)+Sqrt[(-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+
a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2]/(2 Sqrt[2]))^(1/3))-(-18+48 a^2 Λ-5 a^4 Λ^2+
72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+a^2 Λ)^3+
288 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+2 Sqrt[2] √((-18+
48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+
(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+
a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))+(12 M)/(Λ √(-((2 a^2)/3)+2/Λ+(-18+48 a^2 Λ-
5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+3 a^4 Λ^2 Cos[4 θ])/(3 Λ (-3888 M^2 Λ-8 (-3+
a^2 Λ)^3+288 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2)+
2 Sqrt[2] √((-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-36 Λ (-3+a^2 Λ) (-3 (a^2+℧^2)+
a^2 (3+a^2 Λ Cos[θ]^2) Sin[θ]^2))^2))^(1/3))-(1/(3 Λ))((-486 M^2 Λ-(-3+a^2 Λ)^3+
9 Λ (-3+a^2 Λ) (-12 (a^2+℧^2)+12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2)+
(1/(2 Sqrt[2]))(√((-18+48 a^2 Λ-5 a^4 Λ^2+72 Λ ℧^2+36 a^2 Λ Cos[2 θ]+
3 a^4 Λ^2 Cos[4 θ])^3+8 (486 M^2 Λ+(-3+a^2 Λ)^3-9 Λ (-3+a^2 Λ) (-12 (a^2+
℧^2)+12 a^2 Sin[θ]^2+a^4 Λ Sin[2 θ]^2))^2)))^(1/3)))));

RE3[A_, w1_, w2_]:=xyZ[{Sqrt[rE3^2+A^2] Sin[θ] Cos[φ],
Sqrt[rE3^2+A^2] Sin[θ] Sin[φ], rE3 Cos[θ]}, w1];
 
rq[n_]:=Chop[n];
 
horizonsKNdS[A_, mesh_, w1_, w2_]:=Show[
ParametricPlot3D[rq[RE1[A, w1, w2]], {φ, 0, 2 π}, {θ, 0, π},
Mesh->None, PlotPoints->plp, PlotStyle->Directive[Red, Opacity[0.35]]],
ParametricPlot3D[rq[RH1[A, w1, w2]], {φ, 0, 2 π}, {θ, 0, π},
Mesh->None, PlotPoints->plp, PlotStyle->Directive[Cyan, Opacity[0.15]]],
ParametricPlot3D[rq[RH3[A, w1, w2]], {φ, 0, 2 π}, {θ, 0, π},
Mesh->None, PlotPoints->plp, PlotStyle->Directive[Red, Opacity[0.25]]],
ParametricPlot3D[rq[RH2[A, w1, w2]], {φ, 0, 2 π}, {θ, 0, π},
Mesh->None, PlotPoints->plp, PlotStyle->Directive[Gray, Opacity[0.05]]],
ParametricPlot3D[rq[RE3[A, w1, w2]], {φ, 0, 2 π}, {θ, 0, π},
Mesh->None, PlotPoints->plp, PlotStyle->Directive[Blue, Opacity[0.05]]],
ParametricPlot3D[rq[RE2[A, w1, w2]], {φ, 0, 2 π}, {θ, 0, π},
Mesh->None, PlotPoints->plp, PlotStyle->Directive[Blue, Opacity[0.05]]],
If[A==0, {}, ParametricPlot3D[Xyz[xyZ[{Sin[prm] A, Cos[prm] A, 0}, w1], w2],
{prm, 0, 2π}, PlotStyle->{Thickness[0.0025], Orange}]]];

rEΚΝ=1+Sqrt[1-a^2 Cos[θ]^2-℧^2]; 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]; 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]; 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]; RIΚΝ[A_, w1_, w2_]:=Xyz[xyZ[
{Sqrt[rIΚΝ^2+A^2] Sin[θ] Cos[φ], Sqrt[rIΚΝ^2+A^2] Sin[θ] Sin[φ], rIΚΝ Cos[θ]}, w1], w2];

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]]]];

rHRNdS1=1/2 (Sqrt[(-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))]-(1/Sqrt[2])(\[Sqrt]((1/Λ)(8+(-2+8 Λ ℧^2)/(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)-2 (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)-(24 M)/Sqrt[(-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))]))));

rHRNdS2=1/2 (Sqrt[(-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))]+(1/Sqrt[2])(\[Sqrt]((1/Λ)(8+(-2+8 Λ ℧^2)/(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)-2 (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)-(24 M)/Sqrt[(-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))]))));

rHRNdS3=1/2 (-Sqrt[((-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)))]-(1/Sqrt[2])(\[Sqrt]((1/Λ)(8+(-2+8 Λ ℧^2)/(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)-2 (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)+(24 M)/Sqrt[(-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))]))));

rHRNdS4=1/2 (-Sqrt[((-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)))]+(1/Sqrt[2])(\[Sqrt]((1/Λ)(8+(-2+8 Λ ℧^2)/(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)-2 (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3)+(24 M)/Sqrt[(-4 Λ ℧^2+(1+(-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))^2)/(Λ (-1+18 M^2 Λ-12 Λ ℧^2+2 Sqrt[Λ (81 M^4 Λ+℧^2 (3+4 Λ ℧^2)^2-9 M^2 (1+12 Λ ℧^2))])^(1/3))]))));

horizonsRNdS[A_, mesh_, w1_, w2_]:=Show[
Graphics3D[{Red, Opacity[0.2], Sphere[{0,0,0}, rHRNdS1]}],
Graphics3D[{Red, Opacity[0.2], Sphere[{0,0,0}, rHRNdS2]}],
Graphics3D[{Red, Opacity[0.2], Sphere[{0,0,0}, rHRNdS3]}],
Graphics3D[{Red, Opacity[0.2], Sphere[{0,0,0}, rHRNdS4]}]];

rHRSSdS1=(1+(-3 M Sqrt[Λ]+Sqrt[-1+9 M^2 Λ])^(2/3))/(-3 M Λ^2+Sqrt[Λ^3 (-1+9 M^2 Λ)])^(1/3);
rH->(I (I-Sqrt[3]+(I+Sqrt[3]) (-3 M Sqrt[Λ]+Sqrt[-1+9 M^2 Λ])^(2/3)))/(2 (-3 M Λ^2+Sqrt[Λ^3 (-1+9 M^2 Λ)])^(1/3));

rHRSSdS2=((I (-((I+Sqrt[3]) Λ)+(-I+Sqrt[3]) Λ (-3 M Sqrt[Λ]+Sqrt[-1+9 M^2 Λ])^(2/3)))/(2 Λ (-3 M Λ^2+Sqrt[Λ^3 (-1+9 M^2 Λ)])^(1/3)));

horizonsSSdS[A_, mesh_, w1_, w2_]:=Show[
Graphics3D[{Red, Opacity[0.2], Sphere[{0,0,0}, rHSSdS1]}],
Graphics3D[{Red, Opacity[0.2], Sphere[{0,0,0}, rHSSdS2]}]];

rHRN1=1+Sqrt[1-℧^2];
rHRN2=1-Sqrt[1-℧^2];

horizonsRN[A_, mesh_, w1_, w2_]:=Show[
Graphics3D[{Gray, Opacity[0.2], Sphere[{0,0,0}, rHRN1]}],
Graphics3D[{Red, Opacity[0.25], Sphere[{0,0,0}, rHRN2]}]];

horizonsSS[A_, mesh_, w1_, w2_]:=Graphics3D[{Gray, Opacity[0.25], Sphere[{0,0,0}, 2]}]
 
j[v_]:=Sqrt[1-μ^2 v^2];
т[τ_]:=Evaluate[t[τ]/.sol][[1]];
д[ξ_]:=Quiet[zt /.FindRoot[т[zt]-ξ, {zt, 0}]];
 
dp= Style[\!\(\*SuperscriptBox[\(Y\),\(Y\)]\), White]; n0[z_] := Chop[Re[N[Simplify[z]]]];
                                       
DG1={
 
t''[τ]==(2 q (3+a^2 Λ) ℧ (a^2+a^2 Cos[2 θ[τ]]-2 r[τ]^2) (a^2+r[τ]^2) r'[τ])/((a^2+
a^2 Cos[2 θ[τ]]+2 r[τ]^2)^2 (-3 ℧^2+6 M r[τ]-3 r[τ]^2+Λ r[τ]^4+a^2 (-3+Λ r[τ]^2)))-((a^2+
r[τ]^2) (3 a^4 Λ r[τ]+a^4 Λ Cos[4 θ[τ]] r[τ]+4 a^2 (3 M+2 Λ r[τ]^3)+8 r[τ] (3 ℧^2-3 M r[τ]+
Λ r[τ]^4)+4 a^2 Cos[2 θ[τ]] (3 M+Λ r[τ] (a^2+2 r[τ]^2))) r'[τ] t'[τ])/((a^2+a^2 Cos[2 θ[τ]]+
2 r[τ]^2)^2 (6 M r[τ]+Λ r[τ]^4-3 (℧^2+r[τ]^2)+a^2 (-3+Λ r[τ]^2)))-(8 a^2 q (3+
a^2 Λ) ℧ r[τ] Sin[2 θ[τ]] θ'[τ])/((6+a^2 Λ+a^2 Λ Cos[2 θ[τ]]) (a^2+a^2 Cos[2 θ[τ]]+
2 r[τ]^2)^2)+(a^2 (3 (a^4 Λ-8 ℧^2)+a^4 Λ Cos[4 θ[τ]]+4 a^2 Λ Cos[2 θ[τ]] (a^2+2 r[τ]^2)+
8 r[τ] (6 M+Λ r[τ] (a^2+r[τ]^2))) Sin[2 θ[τ]] t'[τ] θ'[τ])/((6+a^2 Λ+
a^2 Λ Cos[2 θ[τ]]) (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^2)+(12 a (a^4 M+4 ℧^2 r[τ]^3-6 M r[τ]^4+
3 a^2 r[τ] (℧^2-M r[τ])+a^2 Cos[2 θ[τ]] (a^2 M+r[τ] (℧^2-
M r[τ]))) Sin[θ[τ]]^2 r'[τ] φ'[τ])/((a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^2 (6 M r[τ]+Λ r[τ]^4-
3 (℧^2+r[τ]^2)+a^2 (-3+Λ r[τ]^2)))-(6 a^3 Cos[θ[τ]] (-℧^2+
2 M r[τ]) Sin[θ[τ]]^3 θ'[τ] φ'[τ])/((3+a^2 Λ Cos[θ[τ]]^2) (a^2+r[τ]^2-a^2 Sin[θ[τ]]^2)^2),

t'[0]==dt,
t[0]==0,
 
r''[τ]==-((r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)-(3 M+r[τ] (-3+a^2 Λ+2 Λ r[τ]^2))/(6 M r[τ]+
Λ r[τ]^4-3 (℧^2+r[τ]^2)+a^2 (-3+Λ r[τ]^2))) r'[τ]^2)-(8 a^2 q ℧ (a^2+a^2 Cos[2 θ[τ]]-
2 r[τ]^2) (-3 ℧^2+6 M r[τ]+a^2 Λ r[τ]^2+Λ r[τ]^4+a^2 Λ Cos[θ[τ]]^2 (a^2+
r[τ]^2)) Sin[θ[τ]]^2 t'[τ])/((3+a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^4)+(8 q ℧ (a^2+
a^2 Cos[2 θ[τ]]-2 r[τ]^2) (a^2+r[τ]^2) (-3 ℧^2+6 M r[τ]-3 r[τ]^2+Λ r[τ]^4+a^2 (-3+Λ r[τ]^2)+
a^2 (3+a^2 Λ Cos[θ[τ]]^2) Sin[θ[τ]]^2) t'[τ])/((3+a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^4)-
(1/(8 (3+a^2 Λ)^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(6 M r[τ]+Λ r[τ]^4-3 (℧^2+r[τ]^2)+a^2 (-3+
Λ r[τ]^2)) (3 a^4 Λ r[τ]+a^4 Λ Cos[4 θ[τ]] r[τ]+4 a^2 (3 M+2 Λ r[τ]^3)+8 r[τ] (3 ℧^2-
3 M r[τ]+Λ r[τ]^4)+4 a^2 Cos[2 θ[τ]] (3 M+Λ r[τ] (a^2+2 r[τ]^2))) t'[τ]^2+
(a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(3 r[τ] (℧^2-2 M r[τ]+(a^2+
r[τ]^2) (1-1/3 Λ r[τ]^2)) θ'[τ]^2)/((3+a^2 Λ Cos[θ[τ]]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2))-
(8 a q ℧ (a^2+a^2 Cos[2 θ[τ]]-2 r[τ]^2) (a^2+r[τ]^2) (-3 ℧^2+6 M r[τ]+a^2 Λ r[τ]^2+Λ r[τ]^4+
a^2 Λ Cos[θ[τ]]^2 (a^2+r[τ]^2)) Sin[θ[τ]]^2 φ'[τ])/((3+a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+
2 r[τ]^2)^4)+(1/((3+a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^4))4 a q ℧ (a^2+a^2 Cos[2 θ[τ]]-
2 r[τ]^2) (a^6 Λ+6 r[τ]^4+3 a^4 (1+Λ r[τ]^2)+a^2 (-3 ℧^2+6 M r[τ]+9 r[τ]^2+2 Λ r[τ]^4)+
a^2 Cos[2 θ[τ]] (a^4 Λ+3 (℧^2-2 M r[τ]+r[τ]^2)+a^2 (3+Λ r[τ]^2))) Sin[θ[τ]]^2 φ'[τ]+
(2 a (6 M r[τ]+Λ r[τ]^4-3 (℧^2+r[τ]^2)+a^2 (-3+Λ r[τ]^2)) (a^4 Λ Cos[θ[τ]]^4 r[τ]+
a^2 Cos[θ[τ]]^2 (3 M+2 Λ r[τ]^3)+r[τ] (3 ℧^2-3 M r[τ]+
Λ r[τ]^4)) Sin[θ[τ]]^2 t'[τ] φ'[τ])/((3+a^2 Λ)^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3)-
(1/(8 (3+a^2 Λ)^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(6 M r[τ]+
Λ r[τ]^4-3 (℧^2+r[τ]^2)+a^2 (-3+Λ r[τ]^2)) (3 a^6 Λ r[τ]+24 r[τ]^5+a^4 Cos[4 θ[τ]] (-3 M+(3+
a^2 Λ) r[τ])+a^4 (3 M+9 r[τ]+8 Λ r[τ]^3)+4 a^2 r[τ] (3 ℧^2-3 M r[τ]+6 r[τ]^2+2 Λ r[τ]^4)+
4 a^2 Cos[2 θ[τ]] r[τ] (a^4 Λ-3 ℧^2+3 r[τ] (M+2 r[τ])+a^2 (3+
2 Λ r[τ]^2))) Sin[θ[τ]]^2 φ'[τ]^2,
 
r'[0]==dR,
r[0]==r0,
 
θ''[τ]==-(1/(96 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))((32 a^2 Cos[θ[τ]] (3+
a^2 Λ Cos[θ[τ]]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] r'[τ]^2)/(℧^2-2 M r[τ]+(a^2+
r[τ]^2) (1-1/3 Λ r[τ]^2))+(24 a^2 q ℧ r[τ] (12 a^2-a^4 Λ+24 ℧^2+12 a^2 Cos[2 θ[τ]]+
a^4 Λ Cos[4 θ[τ]]-48 M r[τ]+24 r[τ]^2-8 a^2 Λ r[τ]^2-8 Λ r[τ]^4) Sin[2 θ[τ]] t'[τ])/((3+
a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2))+(96 a^2 q ℧ r[τ] (a^4 Λ-6 ℧^2+12 M r[τ]+
3 a^2 Λ r[τ]^2+2 Λ r[τ]^4+a^2 Λ Cos[2 θ[τ]] (a^2+r[τ]^2)) Sin[2 θ[τ]] t'[τ])/((3+
a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2))-(3 a^2 (6+a^2 Λ+a^2 Λ Cos[2 θ[τ]]) (3 (a^4 Λ-
8 ℧^2)+a^4 Λ Cos[4 θ[τ]]+4 a^2 Λ Cos[2 θ[τ]] (a^2+2 r[τ]^2)+8 r[τ] (6 M+
Λ r[τ] (a^2+r[τ]^2))) Sin[2 θ[τ]] t'[τ]^2)/(3+
a^2 Λ)^2+192 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+(96 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 (-3+Λ r[τ]^2) Sin[θ[τ]] θ'[τ]^2)/(3+a^2 Λ Cos[θ[τ]]^2)+
(192 a^3 q ℧ Cos[θ[τ]] r[τ] (a^4 Λ-6 ℧^2+12 M r[τ]+3 a^2 Λ r[τ]^2+2 Λ r[τ]^4+
a^2 Λ Cos[2 θ[τ]] (a^2+r[τ]^2)) Sin[θ[τ]]^3 φ'[τ])/((3+a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+
2 r[τ]^2))-(1/((3+a^2 Λ) (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)))96 a q ℧ r[τ] (a^6 Λ+6 r[τ]^4+
3 a^4 (1+Λ r[τ]^2)+a^2 (-3 ℧^2+6 M r[τ]+9 r[τ]^2+2 Λ r[τ]^4)+a^2 Cos[2 θ[τ]] (a^4 Λ+
3 (℧^2-2 M r[τ]+r[τ]^2)+a^2 (3+Λ r[τ]^2))) Sin[2 θ[τ]] φ'[τ]+(6 a (6+a^2 Λ+
a^2 Λ Cos[2 θ[τ]]) (a^2+r[τ]^2) (3 (a^4 Λ-8 ℧^2)+
a^4 Λ Cos[4 θ[τ]]+4 a^2 Λ Cos[2 θ[τ]] (a^2+2 r[τ]^2)+8 r[τ] (6 M+Λ r[τ] (a^2+
r[τ]^2))) Sin[2 θ[τ]] t'[τ] φ'[τ])/(3+a^2 Λ)^2-(1/((3+a^2 Λ)^2))3 (6+a^2 Λ+
a^2 Λ Cos[2 θ[τ]]) (3 a^8 Λ+24 r[τ]^6+a^6 (9+11 Λ r[τ]^2)+a^4 (-15 ℧^2+30 M r[τ]+
33 r[τ]^2+16 Λ r[τ]^4)+8 a^2 r[τ]^2 (-3 ℧^2+Λ r[τ]^4+6 r[τ] (M+r[τ]))+
a^2 (a^2 Cos[4 θ[τ]]+4 Cos[2 θ[τ]] (a^2+2 r[τ]^2)) (3 (℧^2-2 M r[τ]+r[τ]^2)+a^2 (3+
Λ (a^2+r[τ]^2)))) Sin[2 θ[τ]] φ'[τ]^2),
 
θ'[0]==dΘ,
θ[0]==θ0,
 
φ''[τ]==-(1/(6 (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^3))((8 a q (3+a^2 Λ) ℧ (a^2+a^2 Cos[2 θ[τ]]-
2 r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2) r'[τ])/(℧^2-2 M r[τ]+(a^2+r[τ]^2) (1-1/3 Λ r[τ]^2))+
(6 (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2) (3 a^5 Λ r[τ]+a^5 Λ Cos[4 θ[τ]] r[τ]+4 a^3 (3 M+
2 Λ r[τ]^3)+8 a r[τ] (3 ℧^2-3 M r[τ]+Λ r[τ]^4)+4 a^3 Cos[2 θ[τ]] (3 M+
Λ r[τ] (a^2+2 r[τ]^2))) r'[τ] t'[τ])/(6 M r[τ]+Λ r[τ]^4-3 (℧^2+r[τ]^2)+a^2 (-3+Λ r[τ]^2))+
(96 a q (3+a^2 Λ) ℧ Cot[θ[τ]] r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2) θ'[τ])/(3+a^2 Λ Cos[θ[τ]]^2)-
(12 a Cot[θ[τ]] (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2) (3 (a^4 Λ-8 ℧^2)+a^4 Λ Cos[4 θ[τ]]+
4 a^2 Λ Cos[2 θ[τ]] (a^2+2 r[τ]^2)+8 r[τ] (6 M+Λ r[τ] (a^2+r[τ]^2))) t'[τ] θ'[τ])/(6+a^2 Λ+
a^2 Λ Cos[2 θ[τ]])+(6 (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2) (-3 a^4 (M+3 r[τ]-Λ r[τ]^3)+
a^4 Cos[4 θ[τ]] (3 M-3 r[τ]+Λ r[τ]^3)+8 r[τ]^3 (-3 ℧^2+6 M r[τ]-3 r[τ]^2+Λ r[τ]^4)+
4 a^2 r[τ] (-6 ℧^2+9 M r[τ]-6 r[τ]^2+2 Λ r[τ]^4)+4 a^2 Cos[2 θ[τ]] r[τ] (3 (M-2 r[τ]) r[τ]+
2 Λ r[τ]^4+a^2 (-3+Λ r[τ]^2))) r'[τ] φ'[τ])/(6 M r[τ]+Λ r[τ]^4-3 (℧^2+r[τ]^2)+a^2 (-3+
Λ r[τ]^2))-(1/(6+a^2 Λ+a^2 Λ Cos[2 θ[τ]]))3 (a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2) (-32 (3+
a^2 Λ) Cot[θ[τ]] (a^2+r[τ]^2)^2+a^2 (29 a^4 Λ+16 Λ r[τ]^4+24 a^2 (3+2 Λ r[τ]^2)+48 (℧^2+
2 r[τ] (-M+r[τ]))) Sin[2 θ[τ]]+4 a^4 (3+2 Λ (a^2+r[τ]^2)) Sin[4 θ[τ]]+
a^6 Λ Sin[6 θ[τ]]) θ'[τ] φ'[τ]),
 
φ'[0]==dΦ,
φ[0]==φ0
 
};

DGL = If[dgl==1, DG1, If[dgl==2, DG2, If[dgl==3, DG3, If[dgl==4, DG4, DG5]]]];
 
sol=NDSolve[DGL, {t, r, θ, φ}, {τ, 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])}];
 
X[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
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]];
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];
 
Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
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[ψ]};

"Plots"

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

Plot[Mod[180/Pi Θ[tt], 360], {tt, 0, plunge},
Frame->True, PlotStyle->Cyan, AspectRatio->1/5, ImagePadding->{{60, 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->{{60, 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->{{60, 10}, {20, 10}},
ImageSize->600, PlotRange->{{0, plunge}, All}, GridLines->{{}, {0, 1/2, 1}},
PlotLabel -> "v(τ)"]
 
displayP[tp_]:=Grid[{
{If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[tp]], s["GM/c³"], s[dp]},
{s[" t coord"], " = ", s[n0[T[tp]]], s["GM/c²"], s[dp]},
{s[" ṫ coord"], " = ", s[n0[Td[tp]]], s["1"], s[dp]},
{s[" ς lapse"], " = ", s[n0[ς[tp]]], s["1"], s[dp]},
{s[" γ kinet"], " = ", s[n0[1/Sqrt[1-v[tp]^2]]], s["1"], 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[" 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[Rd[tp]]], s["c"], s[dp]},
{s[" d¹φ/dτ¹"], " = ", s[n0[Φd[tp]]], s["c\.b3/G/M"], s[dp]},
{s[" d¹θ/dτ¹"], " = ", s[n0[Θd[tp]]], s["c\.b3/G/M"], s[dp]},

{s[" d²r/dτ²"], " = ", s[n0[Rd'[tp]]], s["c⁴/G/M"], s[dp]},
{s[" d²φ/dτ²"], " = ", s[n0[Φd'[tp]]], s["c⁶/G²/M²"], s[dp]},
{s[" d²θ/dτ²"], " = ", s[n0[Θd'[tp]]], s["c⁶/G²/M²"], s[dp]},

{s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
{s[" Λ cosmo"], " = ", s[n0[Λ]], s["c⁴/G²/M²"], s[dp]},
{s[" H Hubbl"], " = ", s[n0[Sqrt[Λ/3]]], s["c⁶/G²/M²"], s[dp]},
{s[" ℧ cntrl"], " = ", s[n0[℧]], s["M✓(G/kε)"], s[dp]},
{s[" q prtcl"], " = ", s[n0[q]], s["m✓(G/kε)"], s[dp]},

{s[" E total"], " = ", s[n0[e0]], s["mc²"], s[dp]},
{s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},

{s[" ω fdrag"], " = ", s[n0[ω[tp]]], s["c³/G/M"], s[dp]},
{s[" ν fdrag"], " = ", s[n0[ν[tp]]], s["c"], s[dp]},

{s[" v r,loc"], " = ", s[n0[vrj[tp]]], s["c"], s[dp]},
{s[" v θ,loc"], " = ", s[n0[vθj[tp]]], s["c"], s[dp]},
{s[" v φ,loc"], " = ", s[n0[vφj[tp]]], s["c"], s[dp]},
{s[" v local"], " = ", s[n0[v[tp]]], s["c"], s[dp]},
{s[" "], s[" "], s["                       "], s["         "]}},
Alignment-> Left, Spacings-> {0, 0}];
 
hzKNdS=horizonsKNdS[A, None, 0, 0]; hzKN=horizonsΚΝ[A, None, 0, 0];
hzRNdS=horizonsRNdS[A, None, 0, 0]; hzSSdS=horizonsSSdS[A, None, 0, 0];
hzRN=horizonsRN[A, None, 0, 0]; hzSS=horizonsSS[A, None, 0, 0];
 
plot1b[{xx_, yy_, zz_, tk_, w1_, w2_}]:=
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],

If[Λ==0, If[a==0, If[℧==0, hsSS, hzRN],
If[w1==0, hzKN, horizonsΚΝ[A, None, w1, w2]]],
If[w1==0, hzKNdS, horizonsKNdS[A, None, w1, w2]]],

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}];

"Animation"

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/3}]