Das Prinzip des Seins

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 | stable Beta, Version 4C | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

ClearAll["Global`*"]
rA = 1+Sqrt[1-a^2-℧^2];
wp = MachinePrecision;

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

r0   = 50;                                                                 (* Startradius *)
θ0   = π/2;                                                                (* Breitengrad *)
φ0   = 0;                                                                   (* Längengrad *)
a    = 1;                                                                (* Spinparameter *)
℧    = 0;                                       (* spezifische Ladung des schwarzen Lochs *)
v0   = 1;                                                       (* Anfangsgeschwindigkeit *)
μ    = 0;                                                                       (* Photon *)
R0   = 1000;                                               (* Ebene des verzerrten Bildes *)
tmax =-2R0;                                             (* zeitlicher Integrationsbereich *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 2) Metrische Koeffizienten ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Σ = r0^2+a^2 Cos[θ0]^2;
Δ = r0^2-2 r0+a^2+℧^2;
Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;

gtt = (2r0-℧^2)/Σ-1;
grr = Σ/Δ;
gθθ = Σ;
gφφ = Χ/Σ Sin[θ0]^2;
gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 3) Rotationsmatrix für die auf der Sichtebene eintreffenden Strahlen ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

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

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 4) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

raytrace[{yy_,zz_}] :=

Quiet[Module[{vw, tMax, vr0i, vθ0i, vφ0i, vr0a, vθ0a, vφ0a, vt0a, DGL, sol, ε, Lz, pθ0, pr0, Q, k, t10, r10, Θ10, Φ10, т0, R1, t, r, θ, φ, τ, plunge, X, Y, Z, rt, θt, φt, т, ξ, stepsize, laststep, mta},

vw=xyZ[Xyz[{0,1,0},zz],yy+Pi/2];

vr0a = vw[[3]]/Sqrt[grr];
vφ0a = vw[[2]]/Sqrt[gφφ] r0;
vθ0a = vw[[1]]/Sqrt[gθθ] r0;

vt0a = Sqrt[vr0a^2+vφ0a^2+vθ0a^2];

vr0i = vr0a/vt0a;
vφ0i = vφ0a/vt0a;
vθ0i = vθ0a/vt0a;

mta={"EventLocator","Event"->(r[τ]-101/100 rA)};

DGL = {

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

t'[0]==-((a (2 r0-℧^2) Sin[θ0]^2 (vφ0i (-2 r0+r0^2+a^2 Cos[θ0]^2) Csc[θ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) (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 vφ0i (2 r0-℧^2) 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))))/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θ0]^2)))+\[Sqrt](((vr0i^2 (a^2-2 r0+r0^2+℧^2)+vθ0i^2 (a^2-2 r0+r0^2+℧^2)) (r0^2+a^2 Cos[θ0]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θ0]^2)+(a^2 (-2 r0+℧^2)^2 Sin[θ0]^4 (vφ0i (-2 r0+r0^2+a^2 Cos[θ0]^2) Csc[θ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) (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 vφ0i (2 r0-℧^2) 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)/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)^2)-((2 r0-r0^2-℧^2-a^2 Cos[θ0]^2) Sin[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2) (vφ0i (-2 r0+r0^2+a^2 Cos[θ0]^2) Csc[θ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) (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 vφ0i (2 r0-℧^2) 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)/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)^2))/((a^2-2 r0+r0^2+℧^2) (-2 r0+r0^2+℧^2+a^2 Cos[θ0]^2)^2)),
t[0]==0,

r''[τ]==(1/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((8 (a^2 Cos[θ[τ]]^2+(a^2+℧^2-a^2 Cos[θ[τ]]^2) r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) (t'[τ])^2+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] r'[τ] θ'[τ]+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (a^2+℧^2-2 r[τ]+r[τ]^2) (θ'[τ])^2-16 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+(a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 (a^2 (3 a^2+4 ℧^2+4 (a^2-℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]+16 a^2 Cos[θ[τ]]^2 r[τ]^3+8 r[τ]^5-8 a^2 r[τ]^2 Sin[θ[τ]]^2+2 a^4 Sin[2 θ[τ]]^2) (φ'[τ])^2),

r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θ0]^2)/(a^2+(-2+r0) r0+℧^2)],
r[0]==r0,

θ''[τ]==(1/(16 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))-8 a^2 (℧^2-2 r[τ]) Sin[2 θ[τ]] (t'[τ])^2-32 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (θ'[τ])^2+16 a (℧^2-2 r[τ]) (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ] φ'[τ]+(a^4 (3 a^2-5 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+(a^2+℧^2) Cos[4 θ[τ]])+a^2 (11 a^2-8 ℧^2+4 (3 a^2+2 ℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]^2+8 a^2 (2+Cos[2 θ[τ]]) r[τ]^4+8 r[τ]^6+8 a^4 (3+Cos[2 θ[τ]]) r[τ] Sin[θ[τ]]^2+32 a^2 r[τ]^3 Sin[θ[τ]]^2) Sin[2 θ[τ]] (φ'[τ])^2),

θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θ0]^2],
θ[0]==θ0,

φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))(-((8 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 a Cot[θ[τ]] (℧^2-2 r[τ]) t'[τ] θ'[τ]+((a^2 (3 a^2+8 ℧^2+4 a^2 Cos[2 θ[τ]]+a^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) r'[τ] φ'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+Cot[θ[τ]] (a^2 (3 a^2-4 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+16 a^2 Cos[θ[τ]]^2 r[τ]^2+8 r[τ]^4+16 a^2 r[τ] Sin[θ[τ]]^2) θ'[τ] φ'[τ]),

φ'[0]==(vφ0i ((-2+r0) r0+a^2 Cos[θ0]^2) Csc[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0i (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2)),
φ[0]==φ0

};

sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
Method-> mta,
InterpolationOrder-> All,
StepMonitor :> (laststep=plunge; plunge=τ;
stepsize=plunge-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<2*^-2, 0, 1])}];

tMax = Max[tmax, plunge+1/10];

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

rt[τ_] := Evaluate[r[τ]/.sol][[1]];
θt[τ_] := Evaluate[θ[τ]/.sol][[1]]+Pi/2;
φt[τ_] :=-Evaluate[φ[τ]/.sol][[1]]-Pi;

т[coord_,dist_] := Quiet[ξ/.FindRoot[coord[ξ]-dist, {ξ,tMax 9/10,tmax,-1}]];
т0 = т[rt,R0];

R1 = Evaluate[r[т0]/.sol][[1]];
Quiet[If[т0>0, {Pi,Pi/2}, If[Round[R1] == Round[R0],{φt[т0],θt[т0]},If[R1>3,{-Pi,-Pi/2},{Pi,Pi/2}]]]]]]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 5) Testbild laden und transformieren ||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

pic=Import["http://666kb.com/i/dsfr4u76175v8q277.png"]
fpt[{x_,y_}] := {If[y<0,x+1,x],If[y<0,-y,y]};
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi},{-Pi/2,3Pi/2}},PlotRange->{{-Pi,Pi},{-Pi/2,Pi/2}},Padding->"Periodic"]