(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 - 09.09.2020 | Version 20 | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Pause[1]                                                   (* Boyer Lindquist Koordinaten *)
 
wp  = MachinePrecision;
mt0 = Automatic;
mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mta = mt0;                                             (* mt0 for speed, mt1 for accuracy *)
 
kernels = 6;                                                          (* Parallelisierung *)
grain   = 5;                            (* Subparallelisierung auf kernels*grain Streifen *)
rsp     = "Nearest";                                                        (* Resampling *)
 
breite  = 60;                                                (* Zielabmessungen in Pixeln *)
hoehe   = 30;           (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom    = 1;                                  (* doppelter Zoom ergibt halben Sichtwinkel *)
 
LaunchKernels[kernels]
wp = MachinePrecision;                                                     (* Genauigkeit *)
st = 0.1;                                                     (* Auflösung des Gradienten *)
 
pic1 = Import["http://www.yukterez.net/mw/2/flip70.png"];    (* Hintergrundpanorama laden *)
pic2 = Import["http://www.yukterez.net/mw/2/worldmap.png"];    (* Sphärenoberfläche laden *)
pic3 = Import["http://www.yukterez.net/mw/akkretionsscheibe.jpg"];(* Scheibentextur laden *)
pic4 = Import["http://www.yukterez.net/mw/disk.png"];          (* Scheibengeometrie laden *)
pic5 = Import["http://www.yukterez.net/mw/gradient1.png"];   (* Helligkeitsgradient laden *)
pic6 = Import["http://www.yukterez.net/mw/gradient2.png"];          (* Farbgradient laden *)
pic7 = Import["http://www.yukterez.net/mw/2/bw.png"];               (* Checkerboard laden *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 10;                                              (* Radialkoordinate des Beobachters *)
R0 = 500;                                (* Radius des umspannenden Kugelschalenpanoramas *)
R1 = Max[1, Re[1.0001 rA]];                                                     (* Sphäre *)
rT = 10;                                                 (* Lichtlaufzeit Synchronisation *)
t0 = 0;                                                      (* Eigenzeit des Beobachters *)
 
si = isco+0.1;                                           (* Akkretionsscheibe Innenradius *)
sr = 7;                                                  (* Akkretionsscheibe Außenradius *)
 
θ0 = 70 π/180;                                                             (* Breitengrad *)
φ0 = 0;                                                                     (* Längengrad *)
 
tmax = -100 R0;                                         (* zeitlicher Integrationsbereich *)
 
a  = 0.7;                                                                (* Spinparameter *)
℧  = 0.7;                                       (* spezifische Ladung des schwarzen Lochs *)
v0 = 1;                                                   (* Geschwindigkeit der Photonen *)
q  = 0;                                                            (* Ladung der Photonen *)
 
vr = -vя;                                      (* Radiale Geschwindigkeit des Beobachters *)
vϑ = 0;                                         (* Polare Geschwindigkeit des Beobachters *)
vφ = 0;       (* Azimutale Geschwindigkeit des Beobachters: 0 für ZAMO, -й0 für stationär *)

vL = Sqrt[vr^2+vϑ^2+vφ^2];
 
hvs  = 0; (* ArcSin[vφ/vL] *)                        (* horizontaler Versatz in Radianten *)
vvs  = 0; (* ArcSin[vϑ/vL] *)                          (* vertikaler Versatz in Radianten *)

fmax = 2; fmin = 0;                                            (* Doppler Frequenzbereich *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 2) Bildreflektion |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]}
 
pcr1 = ParallelTable[
ImageTransformation[pic1, fpt, DataRange->{{-1, 1}, {0, 1}},
PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Resampling->rsp, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct1 = ImageAssemble[{Table[{pcr1[[x]]}, {x, 2 kernels, 1, -1}]}];
 
pcr2 = ParallelTable[
ImageTransformation[pic2, fpt, DataRange->{{-1, 1}, {0, 1}},
PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Resampling->rsp, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct2 = ImageAssemble[{Table[{pcr2[[x]]}, {x, 2 kernels, 1, -1}]}];
 
pcr3 = ParallelTable[
ImageTransformation[pic7, fpt, DataRange->{{-1, 1}, {0, 1}},
PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Resampling->rsp, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct3 = ImageAssemble[{Table[{pcr3[[x]]}, {x, 2 kernels, 1, -1}]}];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 3) Funktionen |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
gtt = (2r0-℧^2)/Σ-1;
grr = Σ/Δ;
gθθ = Σ;
gφφ = Χ/Σ Sin[θ0]^2;
gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;

rA = 1+Sqrt[1-a^2-℧^2];
rE = 1+Sqrt[1-℧^2];
 
Σ = r0^2+a^2 Cos[θ0]^2;
Δ = r0^2-2 r0+a^2+℧^2;
Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
 
Σs[rs_] := rs^2;
Δs[rs_] := rs^2-2 rs+a^2+℧^2;
Χs[rs_] := (rs^2+a^2)^2-a^2 Δs[rs];
κs[rs_] := a;
 
Σj[rt_, θt_] := rt^2+a^2 Cos[θt]^2;
Δj[rt_, θt_] := rt^2-2 rt+a^2+℧^2;
Χj[rt_, θt_] := (rt^2+a^2)^2-a^2 Sin[θt]^2 Δj[rt, θt];
ωj[rt_, θt_] := (a(2 rt-℧^2))/Χj[rt, θt];
ωR1 = ωj[R1, π/2];
 
ς = Abs@Sqrt[Χ/Δ/Σ]; ςr[rs_] := Abs@Sqrt[Χs[rs]/Δs[rs]/Σs[rs]];
μ = If[Abs[N@v0]==1.0, 0, -1];
 
j[v_] := Sqrt[1-v^2];
Ы[rs_]  := Sqrt[Χs[rs]/Σs[rs]];
ωs[rs_] := (a (2 rs - ℧^2))/Χs[rs];
 
ε[rs_, vt]  := Sqrt[Δs[rs] Σs[rs]/Χs[rs]]/j[vt]+Lz[rs, vt] ωs[rs];
Lz[rs_, vt] := vt Ы[rs]/j[vt];
 
nq[x_] := If[NumericQ[x], If[Element[x, Reals], x, 0], 0];
rl[x_] := If[Abs[Im[x]]>Abs[Re[x]], 1, x]

dΘF = Quiet[If[vL==0, 0, ArcCos[-vϑ/Sqrt[vr^2+vϑ^2+vφ^2]]]];
dθF = Quiet[If[vL==0, 0, If[vr>0, -dΘF, +dΘF]]];
dΦF = Quiet[If[vL==0, 0, ArcTan[Abs[vφ]/vr]]];
dφF = Quiet[If[vL==0, 0, If[vr!=0, -1, 1] If[vφ<0, +1, -1] If[NumericQ[dΦF], dΦF, π/2]]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 4) Geschwindigkeitskomponenten ||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
εj[rt_, δt_, δr_, δθ_, δφ_] := δt (1-(2 rt-℧^2)/rt^2)+(a δφ (2 rt-℧^2))/rt^2;
vrj[rt_, δt_, δr_, δθ_, δφ_] := δr/Sqrt[Δj[rt, θt]] Sqrt[Σj[rt, π/2]];
vθj[rt_, δt_, δr_, δθ_, δφ_] := δθ Sqrt[Σj[rt, π/2]];
vφj[rt_, δt_, δr_, δθ_, δφ_] := (-(((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+
rt^2) Sin[(π/2)] (-δφ-
(a q ℧ rt)/((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+rt^2))+
(εj[rt, δt, δr, δθ, δφ] Csc[(π/2)]^2 (a (-a^2-℧^2+2 rt-rt^2) Sin[(π/2)]^2+a (a^2+
rt^2) Sin[(π/2)]^2))/((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+rt^2))+(a q ℧ rt (a^2+
℧^2-2 rt+rt^2-a^2 Sin[(π/2)]^2))/((a^2 Cos[(π/2)]^2+rt^2)^2 (a^2+℧^2-2 rt+
rt^2))))/((a^2+℧^2-2 rt+rt^2-a^2 Sin[(π/2)]^2) Sqrt[((a^2+rt^2)^2-
a^2 (a^2+℧^2-2 rt+rt^2) Sin[(π/2)]^2)/(a^2 Cos[(π/2)]^2+rt^2)])));
 
vtj[rt_, δt_, δr_, δθ_, δφ_] := Sqrt[vrj[rt, δt, δr, δθ, δφ]^2+
vθj[rt, δt, δr, δθ, δφ]^2+vφj[rt, δt, δr, δθ, δφ]^2];
vφt[rt_, δt_, δr_, δθ_, δφ_] := vφj[rt, δt, δr, δθ, δφ]/vtj[rt, δt, δr, δθ, δφ];
shf[rt_, δt_, δr_, δθ_, δφ_] := Abs[ς/ςr[rt] Sqrt[1-vs[rt]^2]/(1-vs[rt] vφt[rt, δt, δr, δθ, δφ])];
 
dφv[rt_, tt_] := tt φs[rt]/ts[rt];
 
vθ  =-vϑ;
θs  = π/2;
θi  =-θ0+π;
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 5) Photonensphäre und ISCO ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
rp = ц/.Solve[4 a^2 (ц-℧^2)-(ц^2-3 ц+2 ℧^2)^2==0 && ц>=If[Element[rA, Reals], rA, 0], ц];
rP = 1.01 Min[rp]; Rp = 1.01 Max[rp];
isco = Quiet[Min[RI/.NSolve[0 == RI (6 RI-RI^2-9 ℧^2+3 a^2)+4 ℧^2 (℧^2-a^2)-8 a (RI-
℧^2)^(3/2) && RI>=If[Element[rA, Reals], rA, 0], RI]]];
{"r horizon" -> N@rA, "r ergosphere" -> N@rE, "r isco" -> N@isco,
"r photon pro" -> N@Min[rp], "r photon ret" -> N@Max[rp], "r disk" -> N@sr,
"r observer" -> N@r0, "θ observer" -> N@θ0 180/π}
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 6) Geschwindigkeit und Zeitdilatation auf der Scheibe |||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
red[rs_] := Quiet[Reduce[
dt == (Lz[rs, vt] (-a (a^2+rs^2)+Δs[rs] κs[rs])+ε[rs, vt] ((a^2+rs^2)^2-
Δs[rs] κs[rs]^2))/(Δs[rs] Σs[rs])
&&
0 == ((a^2+(-2+rs) rs+℧^2) (16 a dt dΦ rs (rs-℧^2)+8 dt^2 rs (-rs+℧^2)+
dΦ^2 rs (8 rs (-a^2+rs^3)+a^2 (4 a^2+4 ℧^2-4 (a-℧) (a+℧)))))/(8 rs^6)
&&
dΦ == (Lz[rs, vt] (-a^2+Δs[rs])+ε[rs, vt] (a (a^2+rs^2)-Δs[rs] κs[rs]))/(Δs[rs] Σs[rs])
&&
vt > 0,
{vt, dΦ, dt},
Reals]];
 
vs = Interpolation[ParallelTable[{rr, If[Quiet@NumericQ[red[rr][[1, 2]]],
red[rr][[1, 2]], 0]}, {rr, 0, sr+st, st}]];
φs = Interpolation[ParallelTable[{rr, If[Quiet@NumericQ[red[rr][[2, 2]]],
red[rr][[2, 2]], 0]}, {rr, 0, sr+st, st}]];
ts = Interpolation[ParallelTable[{rr, If[Quiet@NumericQ[red[rr][[3, 2]]],
red[rr][[3, 2]], 0]}, {rr, 0, sr+st, st}]];
 
"Accretion disk:"
 
plot[func_, label_] := Plot[func, {rr, si, sr},
GridLines -> {{nq[Min[rp]], nq[Max[rp]], nq[rA], nq[si], nq[isco], nq[rE], nq[sr]}, {}},
Frame -> True, ImagePadding -> {{40, 12}, {12, 12}}, ImageSize -> 340,
PlotLabel -> label, PlotRange->{{0, sr}, All}]
 
Grid[{{
plot[Sqrt[Χs[rr]/Δs[rr]/Σs[rr]],  "Gravitational time dilation: y=dt/dт, x=r"],
plot[φs[rr]/ts[rr], "Shapirodelayed angular velocity: y=dφ/dt, x=r"]},{
plot[ts[rr],  "Total time dilation: y=dt/dτ, x=r"],
plot[φs[rr],  "Coordinate speed: y=dφ/dτ, x=r"]}, {
plot[(a (2 rr-℧^2))/((a^2+rr^2)^2-a^2 (a^2-2 rr+rr^2+℧^2)), "Frame Dragging: y=dφ/dт, x=r"],
plot[vs[rr],  "Local velocity: y=v=dl/dτ, x=r"]}}]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 7) Frame Dragging und Gammafaktor |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
й0 = (a (2 r0-℧^2) Sin[θ0] Sqrt[((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))])/((r0^2+a^2 Cos[θ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)]);
vя = Sqrt[((a^2+r0^2)(℧^2-2r0))/(a^2 Sin[θ0]^2(a^2+(r0-2)r0+℧^2)-(a^2+r0^2)^2)];
vR = Sqrt[(2 r0-℧^2)/(a^2+r0^2)];
 
U = {+vr, +vθ, +vφ};
γ = 1/Sqrt[1-Norm[U]^2];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 8) 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[ψ]};
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 9) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
raytracer[{Ф_, ϑ_}] :=
 
Quiet[Module[{DGL, sol, εj, pθi, pr0, Q, k, V, W, vw,
vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
t10, r10, Θ10, Φ10, t, r, θ, φ, τ,
т, т0, т1, т2, т3, т4, т5,
plunge, plunge0, plunge1, plunge2, plunge3,
plunge4, plunge5, plunge6, plunge7,
dθ0, dφ0, δφ0, δθ0, δr0, δt0, tt0, rt0, θt0, φt0,
dθ1, dφ1, δφ1, δθ1, δr1, δt1, tt1, rt1, θt1, φt1,
dθ2, dφ2, δφ2, δθ2, δr2, δt2, tt2, rt2, θt2, φt2,
dθ3, dφ3, δφ3, δθ3, δr3, δt3, tt3, rt3, θt3, φt3,
dθ4, dφ4, δφ4, δθ4, δr4, δt4, tt4, rt4, θt4, φt4,
dθ5, dφ5, δφ5, δθ5, δr5, δt5, tt5, rt5, θt5, φt5,
dθ6, dφ6, δφ6, δθ6, δr6, δt6, tt6, rt6, θt6, φt6,
X, Y, Z, ξ, stepsize, laststep, mtl, ft, fτ, varb,
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0v, ft1v, ft2v, ft3v, ft4v, ft5v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b, ft5f, ft5x, ft5y,
dt0y, dr0y, dθ0y, dφ0y,
initcon, shF},
 
vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
                                 (* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
vr0a = vw[[3]];
vφ0a = vw[[2]];
vθia = vw[[1]];
                                                                                (* Betrag *)
vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2];
                                                                            (* Normierung *)
vr0n = vr0a/vt0a;
vφ0n = vφ0a/vt0a;
vθ0n = vθia/vt0a;
                                              (* Relativistische Geschwindigkeitsaddition *)
V={vr0n, vθ0n, vφ0n};
W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V);
                                                                            (* Aberration *)
vr0i = W[[1]];
vθ0i = W[[2]];
vφ0i = W[[3]];

shF = Abs[If[vL==0, 1, Sqrt[1-vr^2-vϑ^2-vφ^2]/(1+vL (-Cos[If[vr>0, -dΘF, +dΘF]] Sin[ϑ]-
Sin[If[vr>0, -dΘF, +dΘF]] (Cos[-Ф-π/2+If[vr>0, -1, +1] ArcSin[vφ/vL]] Cos[ϑ] Cos[1/2 π]-
Cos[ϑ] Sin[-Ф-π/2+If[vr>0, -1, +1] ArcSin[vφ/vL]] Sin[1/2 π])))]];
 
initcon = TimeConstrained[NSolve[
dt0y == ς
&&
vr0i ==
Sqrt[(a^2 Cos[θi]^2+r0^2)/(a^2+℧^2-2 r0+r0^2)] Sqrt[1-μ^2 v0^2] dr0y
&&
vθ0i ==
Sqrt[a^2 Cos[θi]^2+r0^2] Sqrt[1-μ^2 v0^2] dθ0y
&&
vφ0i ==
(Sin[θi]^2 Sqrt[1-μ^2 v0^2] (a q μ^2 ℧ r0 v0^2+a (℧^2-
2 r0) dt0y+((a^2+r0^2)^2-a^2 (a^2+℧^2-
2 r0+r0^2) Sin[θi]^2) dφ0y))/((a^2 Cos[θi]^2+r0^2) Sqrt[(Sin[θi]^2 ((a^2+
r0^2)^2-a^2 (a^2+℧^2-2 r0+r0^2) Sin[θi]^2))/(a^2 Cos[θi]^2+r0^2)]),
{dt0y, dr0y, dθ0y, dφ0y}, Reals], 3];
                     
DGL = If[NumericQ[dt0y /. initcon[[1]]] == False, {}, {           (* Bewegungsgleichungen *)
 
t''[τ] == Re[-(1/((a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((q ℧ (a^2 Cos[θ[τ]]^2-
r[τ]^2) (a^2+r[τ]^2) r'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)-(2 (a^2 Cos[θ[τ]]^2+
℧^2 r[τ]-r[τ]^2) (a^2+r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+
a^2 q ℧ r[τ] Sin[2 θ[τ]] θ'[τ]+a^2 (℧^2-2 r[τ]) Sin[2 θ[τ]] 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 r'[τ] φ'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)-
2 a^3 Cos[θ[τ]] (℧^2-2 r[τ]) Sin[θ[τ]]^3 θ'[τ] φ'[τ])],
 
t'[0]  == dt0y/.initcon[[1]],
t[0]   == 0,
 
r''[τ] == Re[-(1/(2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^4))((2 (-a^2 Cos[θ[τ]]^2 (-1+r[τ])+(a^2+
℧^2-r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+r[τ]^2)^3 r'[τ]^2)/(a^2+℧^2-2 r[τ]+r[τ]^2)+
q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2) (a^2+r[τ]^2) (a^2+2 ℧^2+a^2 Cos[2 θ[τ]]-4 r[τ]+
2 r[τ]^2) t'[τ]-2 a^2 q ℧ (℧^2-2 r[τ]) (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2 t'[τ]-
2 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+
r[τ]^2) t'[τ]^2-4 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^3 Sin[θ[τ]] r'[τ] θ'[τ]-
2 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^3 (a^2+℧^2-2 r[τ]+r[τ]^2) θ'[τ]^2-2 a q ℧ (℧^2-
2 r[τ]) (a^2 Cos[θ[τ]]^2-r[τ]^2) (a^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ]+
2 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2) (-(a^2+r[τ]^2)^2 Sin[θ[τ]]^2+a^2 (a^2+
℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^4) φ'[τ]+4 a (a^2 Cos[θ[τ]]^2+(℧^2-
r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+
2 (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 (-r[τ] (-a^4+r[τ]^4+
a^2 (a^2+℧^2-r[τ]) Sin[θ[τ]]^2)+Cos[θ[τ]]^2 (-2 a^2 r[τ] (a^2+r[τ]^2)+a^4 (-1+
r[τ]) Sin[θ[τ]]^2)) φ'[τ]^2)],
 
r'[0]  == dr0y/.initcon[[1]],
r[0]   == r0,
 
θ''[τ] == Re[-(1/(2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))((a^2 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 Sin[2 θ[τ]] r'[τ]^2)/(a^2+℧^2-2 r[τ]+r[τ]^2)+4 r[τ] (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 r'[τ] θ'[τ]-1/8 Sin[2 θ[τ]] (8 r[τ]^6 φ'[τ]^2+16 a r[τ]^3 φ'[τ] (q ℧-
2 t'[τ]+2 a Sin[θ[τ]]^2 φ'[τ])+8 a^2 r[τ]^4 (θ'[τ]^2+(2+Cos[2 θ[τ]]) φ'[τ]^2)+
a^2 (-8 ℧^2 t'[τ]^2+8 a^4 Cos[θ[τ]]^4 θ'[τ]^2+16 a ℧^2 t'[τ] φ'[τ]+a^2 (3 a^2-
5 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+(a^2+℧^2) Cos[4 θ[τ]]) φ'[τ]^2)+
r[τ]^2 (16 a^4 Cos[θ[τ]]^2 θ'[τ]^2+a φ'[τ] (16 ℧^2 t'[τ]+a (11 a^2-
8 ℧^2+4 (3 a^2+2 ℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) φ'[τ]))+
8 a^2 r[τ] (2 t'[τ]^2-2 t'[τ] (q ℧+2 a φ'[τ])+a φ'[τ] (2 q ℧+
a (3+Cos[2 θ[τ]]) Sin[θ[τ]]^2 φ'[τ]))))],
 
θ'[0]  == dθ0y/.initcon[[1]],
θ[0]   == θi,
 
φ''[τ] == Re[-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((4 a q ℧ (a^2 Cos[θ[τ]]^2-
r[τ]^2) r'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)-(8 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-
r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+8 a q ℧ Cot[θ[τ]] r[τ] θ'[τ]+
8 a Cot[θ[τ]] (℧^2-2 r[τ]) t'[τ] θ'[τ]+(1/(a^2+℧^2-2 r[τ]+r[τ]^2))(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'[τ] φ'[τ]+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]  == dφ0y/.initcon[[1]],
φ[0]   == φ0,
 
WhenEvent[Abs[r[τ]] == N[R0]||Abs[r[τ]]>R0 &&
NumericQ[plunge6] == False,
(plunge6=τ) &&
(tt6=t[τ]) && (rt6=r[τ]) && (θt6=θ[τ]) && (φt6=φ[τ]);
"StopIntegration"],
 
WhenEvent[Abs[r[τ]] == N[R1] &&
NumericQ[plunge5] == False,
(plunge5=τ) && (tt5=t[τ]) && (rt5=r[τ]) && (θt5=θ[τ]) && (φt5=φ[τ]);
If[r0>R1, "StopIntegration"]],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && τ < -0.05 &&
NumericQ[plunge0] == False,
(plunge0=τ) &&
(tt0=t[τ]) && (rt0=r[τ]) && (θt0=θ[τ]) && (φt0=φ[τ]) &&
(δt0=t'[τ]) && (δr0=r'[τ]) && (δθ0=θ'[τ]) && (δφ0=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && Re[θ'[τ]]>0 && τ < -0.05 &&
NumericQ[plunge1] == False,
(plunge1=τ) &&
(tt1=t[τ]) && (rt1=r[τ]) && (θt1=θ[τ]) && (φt1=φ[τ]) &&
(δt1=t'[τ]) && (δr1=r'[τ]) && (δθ1=θ'[τ]) && (δφ1=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 &&
r[τ]>N[si] && r[τ]<N[sr] && Re[θ'[τ]]<0 && τ < -0.05 &&
NumericQ[plunge2] == False,
(plunge2=τ) &&
(tt2=t[τ]) && (rt2=r[τ]) && (θt2=θ[τ]) && (φt2=φ[τ]) &&
(δt2=t'[τ]) && (δr2=r'[τ]) && (δθ2=θ'[τ]) && (δφ2=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && Re[θ'[τ]]>0 &&
τ < plunge1-0.05 &&
NumericQ[plunge3] == False && NumericQ[plunge1] == True,
(plunge3=τ) &&
(tt3=t[τ]) && (rt3=r[τ]) && (θt3=θ[τ]) && (φt3=φ[τ]) &&
(δt3=t'[τ]) && (δr3=r'[τ]) && (δθ3=θ'[τ]) && (δφ3=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && Re[θ'[τ]]<0 &&
τ < plunge2-0.05 &&
NumericQ[plunge4] == False && NumericQ[plunge2] == True,
(plunge4=τ) &&
(tt4=t[τ]) && (rt4=r[τ]) && (θt4=θ[τ]) && (φt4=φ[τ]) &&
(δt4=t'[τ]) && (δr4=r'[τ]) && (δθ4=θ'[τ]) && (δφ4=φ'[τ])]
 
}];
                                                                            (* Integrator *)
sol = TimeConstrained[If[NumericQ[dt0y /. initcon[[1]]] == False, {},
NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
Method-> mta,
MaxSteps-> Infinity,
StepMonitor :> (laststep=plunge7; plunge7=τ;
stepsize=plunge7-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<1*^-4, 0, 1])}
]], 10];
 
ft5h=If[NumericQ[plunge5], {φt5-π, θt5+π/2}, {0, -π/2}];
ft5b=If[NumericQ[plunge6], If[rt6<If[Element[rA, Reals], rA, 0], {0, -π/2},
If[rt6>4 R0, {0, -π/2}, {φt6-π, θt6+π/2}]], {0, -π/2}];
ft5f={0, Min[fmax, ς shF]};
ft5x=If[NumericQ[plunge6], {0, 1}, {0, 0}];
ft5y=If[NumericQ[plunge5], {0, 1}, {0, 0}];
 
ft0s=If[NumericQ[plunge0], {φt0, rt0}, {π, sr}];
ft1s=If[NumericQ[plunge1], {φt1, rt1}, {π, sr}];
ft2s=If[NumericQ[plunge2], {φt2, rt2}, {π, sr}];
ft3s=If[NumericQ[plunge3], {φt3, rt3}, {π, sr}];
ft4s=If[NumericQ[plunge4], {φt4, rt4}, {π, sr}];
 
ft0v=If[NumericQ[plunge0], {-dφv[rt0, tt0+t0+rT], 0}, {0, 0}];
ft1v=If[NumericQ[plunge1], {-dφv[rt1, tt1+t0+rT], 0}, {0, 0}];
ft2v=If[NumericQ[plunge2], {-dφv[rt2, tt2+t0+rT], 0}, {0, 0}];
ft3v=If[NumericQ[plunge3], {-dφv[rt3, tt3+t0+rT], 0}, {0, 0}];
ft4v=If[NumericQ[plunge4], {-dφv[rt4, tt4+t0+rT], 0}, {0, 0}];
ft5v=If[NumericQ[plunge5], {-(tt5+t0+rT) ωR1, 0}, {0, 0}];
 
ft0f=If[NumericQ[plunge0], {0, Min[fmax, shF shf[rt0, δt0, δr0, δθ0, δφ0]]}, {0, 0}];
ft1f=If[NumericQ[plunge1], {0, Min[fmax, shF shf[rt1, δt1, δr1, δθ1, δφ1]]}, {0, 0}];
ft2f=If[NumericQ[plunge2], {0, Min[fmax, shF shf[rt2, δt2, δr2, δθ2, δφ2]]}, {0, 0}];
ft3f=If[NumericQ[plunge3], {0, Min[fmax, shF shf[rt3, δt3, δr3, δθ3, δφ3]]}, {0, 0}];
ft4f=If[NumericQ[plunge4], {0, Min[fmax, shF shf[rt4, δt4, δr4, δθ4, δφ4]]}, {0, 0}];
{
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0s+ft0v, ft1s+ft1v, ft2s+ft2v, ft3s+ft3v, ft4s+ft4v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h+ft5v, ft5b, ft5f, ft5x, ft5y
}]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 10) Memoryfunktion ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
mem : raytrace[{Ф_, ϑ_}] := mem = Quiet[Re[raytracer[{Ф, ϑ}]]];
 
ray01[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[01]];
ray02[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[02]];
ray03[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[03]];
ray04[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[04]];
ray05[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[05]];
 
ray06[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[06]];
ray07[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[07]];
ray08[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[08]];
ray09[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[09]];
ray10[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[10]];
 
ray11[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[11]];
ray12[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[12]];
ray13[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[13]];
ray14[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[14]];
ray15[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[15]];
 
ray16[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[16]];
ray17[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[17]];
ray18[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[18]];
ray19[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[19]];
ray20[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[20]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 11) Proportionen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
width1  = ImageDimensions[pic1][[1]]; height1 = ImageDimensions[pic1][[2]];
width2  = ImageDimensions[pic2][[1]]; height2 = ImageDimensions[pic2][[2]];
width3  = ImageDimensions[pic3][[1]]; height3 = ImageDimensions[pic3][[2]];
width4  = ImageDimensions[pic4][[1]]; height4 = ImageDimensions[pic4][[2]];
width5  = ImageDimensions[pic5][[1]]; height5 = ImageDimensions[pic5][[2]];
width6  = ImageDimensions[pic7][[1]]; height6 = ImageDimensions[pic7][[2]];
hzoom   = If[breite>2 hoehe, 1/zoom, 1/zoom/2/hoehe*breite];
vzoom   = If[breite>2 hoehe, 1/zoom*2 hoehe/breite, 1/zoom];
 
"approximate time remaining" -> 1.2 (AbsoluteTiming[Do[raytracer[{
RandomReal[{-π, π}]/zoom, RandomReal[{-π/2, π/2}]/zoom
}], {ü, 1, 50}]][[1]])/50*hoehe*breite/kernels/60 "minutes"
 
FOV->{360.0 hzoom "degree", 180.0 vzoom "degree"}
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 12) Output ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
img = ParallelTable[{
 
(* 1 Hintergrundpanorama *)
Quiet@ImageTransformation[pct1, ray17, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-2π/width1},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+vvs/vzoom+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 2 Sphäre *)
Quiet@ImageTransformation[pct2, ray16, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-3π/2, π/2-2π/width2},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+π/kernels/grain+vvs/vzoom} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 3 Scheibe Textur komplett *)
Quiet@ImageTransformation[pic3, ray01, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 4 Scheibe Textur Front *)
Quiet@ImageTransformation[pic3, ray02, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 5 Scheibe Textur Echo 1 *)
Quiet@ImageTransformation[pic3, ray03, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 6 Scheibe Textur Echo 2 *)
Quiet@ImageTransformation[pic3, ray04, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 7 Scheibe Textur Echo 3 *)
Quiet@ImageTransformation[pic3, ray05, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 8 Scheibe Geometrie still komplett *)
Quiet@ImageTransformation[pic4, ray01, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 9 Scheibe Geometrie still Front *)
Quiet@ImageTransformation[pic4, ray02, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 10 Scheibe Geometrie still Echo 1 *)
Quiet@ImageTransformation[pic4, ray03, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 11 Scheibe Geometrie still Echo 2 *)
Quiet@ImageTransformation[pic4, ray04, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 12 Scheibe Geometrie still Echo 3 *)
Quiet@ImageTransformation[pic4, ray05, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 13 Scheibe Geometrie rotierend komplett *)
Quiet@ImageTransformation[pic4, ray06, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 14 Scheibe Geometrie rotierend Front *)
Quiet@ImageTransformation[pic4, ray07, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 15 Scheibe Geometrie rotierend Echo 1 *)
Quiet@ImageTransformation[pic4, ray08, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 16 Scheibe Geometrie rotierend Echo 2 *)
Quiet@ImageTransformation[pic4, ray09, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 17 Scheibe Geometrie rotierend Echo 3 *)
Quiet@ImageTransformation[pic4, ray10, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 18 Frequenzverschiebung komplett SW *)
Quiet@ImageTransformation[pic5, ray11, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-1, 1},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 19 Frequenzverschiebung Front SW *)
Quiet@ImageTransformation[pic5, ray12, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 20 Frequenzverschiebung Echo 1 SW *)
Quiet@ImageTransformation[pic5, ray13, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 21 Frequenzverschiebung Echo 2 SW *)
Quiet@ImageTransformation[pic5, ray14, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 22 Frequenzverschiebung Echo 3 SW *)
Quiet@ImageTransformation[pic5, ray15, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 23 Frequenzverschiebung komplett Farbe *)
Quiet@ImageTransformation[pic6, ray11, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 24 Frequenzverschiebung Front Farbe *)
Quiet@ImageTransformation[pic6, ray12, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 25 Frequenzverschiebung Echo 1 Farbe *)
Quiet@ImageTransformation[pic6, ray13, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 26 Frequenzverschiebung Echo 2 Farbe *)
Quiet@ImageTransformation[pic6, ray14, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 27 Frequenzverschiebung Echo 3 Farbe *)
Quiet@ImageTransformation[pic6, ray15, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 28 Hintergrundpanorama *)
ImageTransformation[pct3, ray17, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-2π/width6},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+vvs/vzoom+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 29 Sphäre *)
ImageTransformation[pct3, ray16, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-2π/width6},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+π/kernels/grain+vvs/vzoom} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 30 Scheibe Textur komplett *)
ImageTransformation[pic7, ray01, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width6},
{si, sr+(sr-si)/height6}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],

(* 31 Frequenzverschiebung Hintergrund Farbe *)
Quiet@ImageTransformation[pic6, ray18, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],

(* 32 Frequenzverschiebung Hintergrund SW *)
Quiet@ImageTransformation[pic5, ray18, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],

(* 33 Kontur BH *)
Quiet@ImageTransformation[pic5, ray19, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 1}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],

(* 34 Kontur IH *)
Quiet@ImageTransformation[pic5, ray20, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 1}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"]
 
}, {x, 0, π-π/kernels/grain, π/kernels/grain}];
 
image[num_] := ImageAssemble[{Table[{img[[x, num]]}, {x, kernels grain, 1, -1}]}];
fi[x_] := ColorNegate[x];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 13) Composite |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ParallelTable[image[n], {n, 1, 34, 1}]
 
ds1 = image[8]
ds2 = ImageMultiply[fi@image[09], fi@image[10]];
ds3 = ImageMultiply[fi@image[11], fi@image[12]];
ds4 = fi@ImageMultiply[ds2, ds3]
 
ds5 = image[13]
ds6 = ImageMultiply[fi@image[14], fi@image[15]];
ds7 = ImageMultiply[fi@image[16], fi@image[17]];
ds8 = fi@ImageMultiply[ds6, ds7]

bf1 = image[31]
bf2 = image[32]
bf3 = image[33]
bf4 = image[34]
 
bl1 = image[23]
bl2 = ImageMultiply[fi@image[24], fi@image[25]];
bl3 = ImageMultiply[fi@image[26], fi@image[27]];
bl4 = fi@ImageMultiply[bl2, bl3]
 
gr1 = image[18]
gr2 = ImageMultiply[fi@image[19], fi@image[20]];
gr3 = ImageMultiply[fi@image[21], fi@image[22]];
gr4 = fi@ImageMultiply[gr2, gr3]
 
sea = ImageMultiply[image[04], image[19]];
seb = ImageMultiply[image[05], image[20]];
sec = ImageMultiply[image[06], image[21]];
sed = ImageMultiply[image[07], image[22]];
 
se1 = image[03]
se2 = ImageMultiply[fi@image[04], fi@image[05]];
se3 = ImageMultiply[fi@image[06], fi@image[07]];
se4 = fi@ImageMultiply[se2, se3]
se5 = ImageMultiply[fi@sea, fi@seb];
se6 = ImageMultiply[fi@sec, fi@sed];
se7 = fi@ImageMultiply[se5, se6]
se8 = fi@ImageMultiply[fi@se7, fi@se7]
se9 = fi@ImageMultiply[fi@se8, fi@se7]
 
bkg = image[1]
hrz = image[2]
 
cb1 = image[28]
cb2 = image[29]
cb3 = image[30]
 
Quit[]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 14) Alternative initcon dt/dτ falls Partikel statt Photonen verwendet werden ||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
(* initcon = NSolve[
-μ  == -(((r0^2+a^2 Cos[θ0]^2) dr0y^2)/(a^2-2 r0+r0^2+℧^2))+((a^2-2 r0+
r0^2+℧^2-a^2 Sin[θ0]^2) (dt0y)^2)/(r0^2+a^2 Cos[θ0]^2)+(-r0^2-
a^2 Cos[θ0]^2) dθ0y^2+(2 a (2 r0-℧^2) Sin[θ0]^2 dt0y dφ0y)/(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φ0y^2)/(r0^2+a^2 Cos[θ0]^2)
&&
dt0y > 0
&&
vr0i ==
Sqrt[(a^2 Cos[θi]^2+r0^2)/(a^2+℧^2-2 r0+r0^2)] Sqrt[1-μ^2 v0^2] dr0y
&&
vθ0i ==
Sqrt[a^2 Cos[θi]^2+r0^2] Sqrt[1-μ^2 v0^2] dθ0y
&&
vφ0i ==
(Sin[θi]^2 Sqrt[1-μ^2 v0^2] (a q μ^2 ℧ r0 v0^2+a (℧^2-
2 r0) dt0y+((a^2+r0^2)^2-a^2 (a^2+℧^2-
2 r0+r0^2) Sin[θi]^2) dφ0y))/((a^2 Cos[θi]^2+r0^2) Sqrt[(Sin[θi]^2 ((a^2+
r0^2)^2-a^2 (a^2+℧^2-2 r0+r0^2) Sin[θi]^2))/(a^2 Cos[θi]^2+r0^2)]),
{dt0y, dr0y, dθ0y, dφ0y}, Reals];                                                         *)

























(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 - 09.09.2020 | Version 20 | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
Pause[1]                                                    (* Raindrop Doran Koordinaten *)
 
wp  = MachinePrecision;
mt0 = Automatic;
mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mta = mt0;                                             (* mt0 for speed, mt1 for accuracy *)
 
kernels = 6;                                                          (* Parallelisierung *)
grain   = 5;                            (* Subparallelisierung auf kernels*grain Streifen *)
rsp     = "Nearest";                                                        (* Resampling *)
 
breite  = 60;                                                (* Zielabmessungen in Pixeln *)
hoehe   = 30;           (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom    = 1;                                  (* doppelter Zoom ergibt halben Sichtwinkel *)
 
LaunchKernels[kernels]
wp = MachinePrecision;                                                     (* Genauigkeit *)
st = 0.1;                                                     (* Auflösung des Gradienten *)
 
pic1 = Import["http://www.yukterez.net/mw/2/flip70.png"];    (* Hintergrundpanorama laden *)
pic2 = Import["http://www.yukterez.net/mw/2/worldmap.png"];    (* Sphärenoberfläche laden *)
pic3 = Import["http://www.yukterez.net/mw/akkretionsscheibe.jpg"];(* Scheibentextur laden *)
pic4 = Import["http://www.yukterez.net/mw/disk.png"];          (* Scheibengeometrie laden *)
pic5 = Import["http://www.yukterez.net/mw/gradient1.png"];   (* Helligkeitsgradient laden *)
pic6 = Import["http://www.yukterez.net/mw/gradient2.png"];          (* Farbgradient laden *)
pic7 = Import["http://www.yukterez.net/mw/2/bw.png"];               (* Checkerboard laden *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0 = 10;                                              (* Radialkoordinate des Beobachters *)
R0 = 500;                                (* Radius des umspannenden Kugelschalenpanoramas *)
R1 = Max[1, Re[1.0001 rA]];                                                     (* Sphäre *)
rT = 10;                                                 (* Lichtlaufzeit Synchronisation *)
t0 = 0;                                                      (* Eigenzeit des Beobachters *)

si = isco+0.1;                                           (* Akkretionsscheibe Innenradius *)
sr = 7;                                                  (* Akkretionsscheibe Außenradius *)

θ0 = 70 π/180;                                                             (* Breitengrad *)
φ0 = 0;                                                                     (* Längengrad *)
 
tmax = -100 R0;                                         (* zeitlicher Integrationsbereich *)
 
a  = 0.7;                                                                (* Spinparameter *)
℧  = 0.7;                                       (* spezifische Ladung des schwarzen Lochs *)
v0 = 1;                                                   (* Geschwindigkeit der Photonen *)
 
vr = 0;                                   (* Radiale Geschwindigkeit relativ zum Raindrop *)
vϑ = 0;                                         (* Polare Geschwindigkeit des Beobachters *)
vφ = 0;       (* Azimutale Geschwindigkeit des Beobachters: 0 für ZAMO, -й0 für stationär *)

vL = Sqrt[vr^2+vϑ^2+vφ^2];
 
hvs  = 0; (* ArcSin[vφ/vL] *)                        (* horizontaler Versatz in Radianten *)
vvs  = 0; (* ArcSin[vϑ/vL] *)                          (* vertikaler Versatz in Radianten *)

fmax = 2; fmin = 0;                                            (* Doppler Frequenzbereich *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 2) Bildreflektion |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]}
 
pcr1 = ParallelTable[
Quiet@ImageTransformation[pic1, fpt, DataRange->{{-1, 1}, {0, 1}},
PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}},
Resampling->rsp, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct1 = ImageAssemble[{Table[{pcr1[[x]]}, {x, 2 kernels, 1, -1}]}];
 
pcr2 = ParallelTable[
Quiet@ImageTransformation[pic2, fpt, DataRange->{{-1, 1}, {0, 1}},
PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}},
Resampling->rsp, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct2 = ImageAssemble[{Table[{pcr2[[x]]}, {x, 2 kernels, 1, -1}]}];
 
pcr3 = ParallelTable[
ImageTransformation[pic7, fpt, DataRange->{{-1, 1}, {0, 1}},
PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Resampling->rsp, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct3 = ImageAssemble[{Table[{pcr3[[x]]}, {x, 2 kernels, 1, -1}]}];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 3) Funktionen |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
g11 = gtt = 1+(-4 r0+2 ℧^2)/(a^2+2 r0^2+a^2 Cos[2 θ0]);
g22 = grr = -((a^2+2 r0^2+a^2 Cos[2 θ0])/(2 (a^2+r0^2)));
g33 = gθθ = -r0^2-a^2 Cos[θ0]^2;
g44 = gφφ = (-(a^2+r0^2)^2 Sin[θ0]^2+a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^4)/(r0^2+a^2 Cos[θ0]^2);
g12 = gtr = -(Sqrt[2 r0-℧^2]/Sqrt[a^2+r0^2]);
g14 = gtφ = (a (2 r0-℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2);
g24 = grφ = (a Sqrt[2 r0-℧^2] Sin[θ0]^2)/Sqrt[a^2+r0^2];
g13 = g23 = g34 = 0;

rA  = 1+Sqrt[1-a^2-℧^2];
rI  = 1-Sqrt[1-a^2-℧^2];
rE  = 1+Sqrt[1-℧^2-a^2 Cos[θ0]^2];
rG  = 1-Sqrt[1-℧^2-a^2 Cos[θ0]^2];
 
q = 0;
 
Σ = r0^2+a^2 Cos[θ0]^2;
Δ = r0^2-2 r0+a^2+℧^2;
Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
 
Σs[rs_] := rs^2;
Δs[rs_] := rs^2-2 rs+a^2+℧^2;
Χs[rs_] := (rs^2+a^2)^2-a^2 Δs[rs];
κs[rs_] := a;
 
Σj[rt_, θt_] := rt^2+a^2 Cos[θt]^2;
Δj[rt_, θt_] := rt^2-2 rt+a^2+℧^2;
Χj[rt_, θt_] := (rt^2+a^2)^2-a^2 Sin[θt]^2 Δj[rt, θt];
ωj[rt_, θt_] := (a(2 rt-℧^2))/Χj[rt, θt];
ωR1 = ωj[R1, π/2];
 
ς = Abs@Sqrt[Χ/Δ/Σ]; ςr[rs_] := Abs@Sqrt[Χs[rs]/Δs[rs]/Σs[rs]];
μ = If[Abs[N@v0]==1.0, 0, -1];
 
j[v_] := Sqrt[1-v^2];
Ы[rs_]  := Sqrt[Χs[rs]/Σs[rs]];
ωs[rs_] := (a (2 rs - ℧^2))/Χs[rs];
 
ε[rs_, vt_]  := Sqrt[Δs[rs] Σs[rs]/Χs[rs]]/j[vt]+Lz[rs, vt] ωs[rs];
Lz[rs_, vt_] := vt Ы[rs]/j[vt];
 
nq[x_] := If[NumericQ[x], If[Element[x, Reals], x, 0], 0];

dΘF = Quiet[If[vLF==0, 0, ArcCos[-vθF/Sqrt[vrF^2+vθF^2+vφF^2]]]];
dθF = Quiet[If[vLF==0, 0, If[vrF>0, -dΘF, +dΘF]]];
dΦF = Quiet[If[vLF==0, 0, ArcTan[Abs[vφF]/vrF]]];
dφF = Quiet[If[vLF==0, 0, If[vrF!=0, -1, 1] If[vφF<0, +1, -1] If[NumericQ[dΦF], dΦF, π/2]]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 4) Geschwindigkeitskomponenten ||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
εj[rt_, δt_, δr_, δθ_, δφ_] := (δt-δr Sqrt[(2rt-℧^2)/(a^2+rt^2)]/(1-(2rt-℧^2)/(a^2+
rt^2))) (1-(2 rt-℧^2)/rt^2)+(a (δφ-δr a Sqrt[(2rt-℧^2)/(a^2+rt^2)]/(1-(2rt-℧^2)/(a^2+
rt^2))/(a^2+rt^2)) (2 rt-℧^2))/rt^2;

vrj[rt_, δt_, δr_, δθ_, δφ_] := δr/Sqrt[Δj[rt, θt]] Sqrt[Σj[rt, π/2]];
vθj[rt_, δt_, δr_, δθ_, δφ_] := δθ Sqrt[Σj[rt, π/2]];
vφj[rt_, δt_, δr_, δθ_, δφ_] := (-(((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+
rt^2) Sin[(π/2)] (-(δφ-δr a Sqrt[(2rt-℧^2)/(a^2+rt^2)]/(1-(2rt-℧^2)/(a^2+rt^2))/(a^2+rt^2))-
(a q ℧ rt)/((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+rt^2))+
(εj[rt, δt, δr, δθ, δφ] Csc[(π/2)]^2 (a (-a^2-℧^2+2 rt-rt^2) Sin[(π/2)]^2+a (a^2+
rt^2) Sin[(π/2)]^2))/((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+rt^2))+(a q ℧ rt (a^2+
℧^2-2 rt+rt^2-a^2 Sin[(π/2)]^2))/((a^2 Cos[(π/2)]^2+rt^2)^2 (a^2+℧^2-2 rt+
rt^2))))/((a^2+℧^2-2 rt+rt^2-a^2 Sin[(π/2)]^2) Sqrt[((a^2+rt^2)^2-
a^2 (a^2+℧^2-2 rt+rt^2) Sin[(π/2)]^2)/(a^2 Cos[(π/2)]^2+rt^2)])));

vtj[rt_, δt_, δr_, δθ_, δφ_] := Sqrt[vrj[rt, δt, δr, δθ, δφ]^2+
vθj[rt, δt, δr, δθ, δφ]^2+vφj[rt, δt, δr, δθ, δφ]^2];
vφt[rt_, δt_, δr_, δθ_, δφ_] := vφj[rt, δt, δr, δθ, δφ]/vtj[rt, δt, δr, δθ, δφ];
shf[rt_, δt_, δr_, δθ_, δφ_] := Abs[ς/ςr[rt] Sqrt[1-vs[rt]^2]/(1-vs[rt] vφt[rt, δt, δr, δθ, δφ])];

dφv[rt_, tt_] := tt φs[rt]/ts[rt];

vsolve = Block[{r=r0, θ=θ0, μ=-1, q=0},

ini = NSolve[
vr ==
(Sqrt[2] Sqrt[1-μ^2 vL^2] ((q μ^2 ℧ r Sqrt[-℧^2+2 r] vL^2)/((a^2+℧^2+
(-2+r) r) Sqrt[a^2+r^2])+((a^2+a^2 Cos[2 θ]+2 r^2) dR)/(2 (a^2+r^2))+(Sqrt[-℧^2+2 r] (dT-
a Sin[θ]^2 dΦ))/Sqrt[a^2+r^2]))/Sqrt[(a^2+a^2 Cos[2 θ]+2 r^2)/(a^2+r^2)]
&&
vϑ ==
Sqrt[a^2 Cos[θ]^2+r^2] Sqrt[1-μ^2 vL^2] dΘ
&&
vφ ==
(Sin[θ]^2 Sqrt[1-μ^2 vL^2] (a q μ^2 ℧ r Sqrt[a^2+r^2] vL^2-1/2 a Sqrt[-℧^2+2 r] (a^2+
a^2 Cos[2 θ]+2 r^2) dR+Sqrt[a^2+r^2] (a (℧^2-2 r) dT+((a^2+r^2)^2-a^2 (a^2+℧^2-2 r+
r^2) Sin[θ]^2) dΦ)))/(Sqrt[a^2+r^2] (a^2 Cos[θ]^2+r^2) Sqrt[(Sin[θ]^2 ((a^2+r^2)^2-
a^2 (a^2+℧^2-2 r+r^2) Sin[θ]^2))/(a^2 Cos[θ]^2+r^2)])
&&
-μ ==
-(((a^2+2 r^2+
a^2 Cos[2 θ]) (dR)^2)/(2 (a^2+r^2)))-(2 Sqrt[2 r-℧^2] dR dT)/Sqrt[a^2+r^2]+(1+(-4 r+
2 ℧^2)/(a^2+2 r^2+a^2 Cos[2 θ])) (dT)^2+(-r^2-a^2 Cos[θ]^2) dΘ^2+(2 a Sqrt[2 r-
℧^2] Sin[θ]^2 dR dΦ)/Sqrt[a^2+r^2]+(2 a (2 r-℧^2) Sin[θ]^2 dT dΦ)/(r^2+a^2 Cos[θ]^2)+
((-(a^2+r^2)^2 Sin[θ]^2+a^2 (a^2+(-2+r) r+℧^2) Sin[θ]^4) (dΦ)^2)/(r^2+a^2 Cos[θ]^2)
&&
dT > 0,
{dT, dR, dΘ, dΦ},
Reals];

β = -Sqrt[(2r-℧^2)/(a^2+r^2)];

drf = (dR/.ini[[1]]);
dθf = (dΘ/.ini[[1]]);
dφf = (dΦ/.ini[[1]])+If[r0<rE&&r0>rG, +1, -1] drf a β/(1-β^2)/(a^2+r^2);

NSolve[
vrf ==
Sqrt[(a^2 Cos[θ]^2+r^2)/(a^2+℧^2-2 r+r^2)] Sqrt[1-μ^2 v0f^2] drf
&&
vϑf == Sqrt[a^2 Cos[θ]^2+r^2] Sqrt[1-μ^2 v0f^2] dθf
&&
vφf == (Sin[θ]^2 Sqrt[1-μ^2 v0f^2] (a q μ^2 ℧ r v0f^2+a (℧^2-2 r) dtf+((a^2+r^2)^2-
a^2 (a^2+℧^2-2 r+r^2) Sin[θ]^2) dφf))/((a^2 Cos[θ]^2+r^2) Sqrt[(Sin[θ]^2 ((a^2+
r^2)^2-a^2 (a^2+℧^2-2 r+r^2) Sin[θ]^2))/(a^2 Cos[θ]^2+r^2)])
&&
v0f == Sqrt[vrf^2+vϑf^2+vφf^2]
&&
-μ == -(((r^2+a^2 Cos[θ]^2) drf^2)/(a^2-2 r+r^2+℧^2))+((a^2-2 r+r^2+℧^2-
a^2 Sin[θ]^2) (dtf)^2)/(r^2+a^2 Cos[θ]^2)+(-r^2-a^2 Cos[θ]^2) dθf^2+(2 a (2 r-
℧^2) Sin[θ]^2 dtf dφf)/(r^2+a^2 Cos[θ]^2)+((-(a^2+r^2)^2 Sin[θ]^2+a^2 (a^2-2 r+
r^2+℧^2) Sin[θ]^4) dφf^2)/(r^2+a^2 Cos[θ]^2),
{dtf, vrf, vϑf, vφf, v0f}]];

vrFx = vrf/.vsolve[[1]];
vθFx = vϑf/.vsolve[[1]];
vφFx = vφf/.vsolve[[1]];

vθF = Re[vθFx]+Im[vθFx];
vφF = Re[vφFx]+Im[vφFx];
vLF  = Sqrt[vrFx^2+vθF^2+vφF^2];
vrF = If[vLF<1,vrFx,-vrFx];

vθ  =-vϑ;
θs  = π/2;
θi  =-θ0+π;
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 5) Photonensphäre und ISCO ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
rp = ц/.Solve[4 a^2 (ц-℧^2)-(ц^2-3 ц+2 ℧^2)^2==0 && ц>=If[Element[rA, Reals], rA, 0], ц];
rP = 1.01 Min[rp]; Rp = 1.01 Max[rp];
isco = Quiet[Min[RI/.NSolve[0 == RI (6 RI-RI^2-9 ℧^2+3 a^2)+4 ℧^2 (℧^2-a^2)-8 a (RI-
℧^2)^(3/2) && RI>=If[Element[rA, Reals], rA, 0], RI]]];
{"r horizon" -> N@rA, "r ergosphere" -> N@rE, "r isco" -> N@isco,
"r photon pro" -> N@Min[rp], "r photon ret" -> N@Max[rp], "r disk" -> N@sr,
"r observer" -> N@r0, "θ observer" -> N@θ0 180/π}
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 6) Geschwindigkeit und Zeitdilatation auf der Scheibe |||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
red[rs_] := Quiet[Reduce[
dt == (Lz[rs, vt] (-a (a^2+rs^2)+Δs[rs] κs[rs])+ε[rs, vt] ((a^2+rs^2)^2-
Δs[rs] κs[rs]^2))/(Δs[rs] Σs[rs])
&&
0 == ((a^2+(-2+rs) rs+℧^2) (16 a dt dΦ rs (rs-℧^2)+8 dt^2 rs (-rs+℧^2)+
dΦ^2 rs (8 rs (-a^2+rs^3)+a^2 (4 a^2+4 ℧^2-4 (a-℧) (a+℧)))))/(8 rs^6)
&&
dΦ == (Lz[rs, vt] (-a^2+Δs[rs])+ε[rs, vt] (a (a^2+rs^2)-Δs[rs] κs[rs]))/(Δs[rs] Σs[rs])
&&
vt > 0,
{vt, dΦ, dt},
Reals]];
 
vs = Interpolation[ParallelTable[{rr, If[Quiet@NumericQ[red[rr][[1, 2]]],
red[rr][[1, 2]], 0]}, {rr, 0, sr+st, st}]];
φs = Interpolation[ParallelTable[{rr, If[Quiet@NumericQ[red[rr][[2, 2]]],
red[rr][[2, 2]], 0]}, {rr, 0, sr+st, st}]];
ts = Interpolation[ParallelTable[{rr, If[Quiet@NumericQ[red[rr][[3, 2]]],
red[rr][[3, 2]], 0]}, {rr, 0, sr+st, st}]];
 
"Accretion disk:"
 
plot[func_, label_] := Plot[func, {rr, si, sr},
GridLines -> {{nq[Min[rp]], nq[Max[rp]], nq[rA], nq[si], nq[isco], nq[rE], nq[sr]}, {}},
Frame -> True, ImagePadding -> {{40, 12}, {12, 12}}, ImageSize -> 340,
PlotLabel -> label, PlotRange->{{0, sr}, All}]
 
Grid[{{
plot[Sqrt[Χs[rr]/Δs[rr]/Σs[rr]],  "Gravitational time dilation: y=dt/dт, x=r"],
plot[φs[rr]/ts[rr], "Shapirodelayed angular velocity: y=dφ/dt, x=r"]},{
plot[ts[rr],  "Total time dilation: y=dt/dτ, x=r"],
plot[φs[rr],  "Coordinate speed: y=dφ/dτ, x=r"]}, {
plot[(a (2 rr-℧^2))/((a^2+rr^2)^2-a^2 (a^2-2 rr+rr^2+℧^2)), "Frame Dragging: y=dφ/dт, x=r"],
plot[vs[rr],  "Local velocity: y=v=dl/dτ, x=r"]}}]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 7) Frame Dragging und Gammafaktor |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
й0 = (a (2 r0-℧^2) Sin[θ0] Sqrt[((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))])/((r0^2+a^2 Cos[θ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)]);
vя = Sqrt[((a^2+r0^2)(℧^2-2r0))/(a^2 Sin[θ0]^2(a^2+(r0-2)r0+℧^2)-(a^2+r0^2)^2)];
vR = Sqrt[(2 r0-℧^2)/(a^2+r0^2)];
 
U = {+vr, +vθ, +vφ};
γ = 1/Sqrt[1-Norm[U]^2];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 8) 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[ψ]};
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 9) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
raytracer[{Ф_, ϑ_}] :=
 
Quiet[Module[{DGL, sol, εj, pθi, pr0, Q, k, V, W, vw,
vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
t10, r10, Θ10, Φ10, t, r, θ, φ, τ,
т, т0, т1, т2, т3, т4, т5,
plunge, plunge0, plunge1, plunge2, plunge3,
plunge4, plunge5, plunge6, plunge7,
dθ0, dφ0, δφ0, δθ0, δr0, δt0, tt0, rt0, θt0, φt0,
dθ1, dφ1, δφ1, δθ1, δr1, δt1, tt1, rt1, θt1, φt1,
dθ2, dφ2, δφ2, δθ2, δr2, δt2, tt2, rt2, θt2, φt2,
dθ3, dφ3, δφ3, δθ3, δr3, δt3, tt3, rt3, θt3, φt3,
dθ4, dφ4, δφ4, δθ4, δr4, δt4, tt4, rt4, θt4, φt4,
dθ5, dφ5, δφ5, δθ5, δr5, δt5, tt5, rt5, θt5, φt5,
dθ6, dφ6, δφ6, δθ6, δr6, δt6, tt6, rt6, θt6, φt6,
X, Y, Z, ξ, stepsize, laststep, mtl, ft, fτ, varb,
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0v, ft1v, ft2v, ft3v, ft4v, ft5v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b, ft5f, ft5x, ft5y,
dt0y, dr0y, dθ0y, dφ0y,
initcon, shF},
 
vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
                                 (* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
vr0a = vw[[3]];
vφ0a = vw[[2]];
vθia = vw[[1]];
                                                                                (* Betrag *)
vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2];
                                                                            (* Normierung *)
vr0n = vr0a/vt0a;
vφ0n = vφ0a/vt0a;
vθ0n = vθia/vt0a;
                                              (* Relativistische Geschwindigkeitsaddition *)
V={vr0n, vθ0n, vφ0n};
W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V);
                                                                            (* Aberration *)
vr0i = W[[1]];
vθ0i = W[[2]];
vφ0i = W[[3]];

shF = Abs[If[vLF==0, 1, Sqrt[1-vLF^2]/(1+vLF (-Cos[If[vrF>0, -dΘF, +dΘF]] Sin[ϑ]-
Sin[If[vrF>0, -dΘF, +dΘF]] (Cos[-Ф-π/2+If[vrF>0, -1,
1] ArcSin[vφF/vLF]] Cos[ϑ] Cos[1/2 π If[vφF<0,
1, -1] If[vrF!=0,-1,1]]-Cos[ϑ] Sin[-Ф-π/2+
If[vrF>0, -1, +1]ArcSin[vφF/vLF]] Sin[1/2 π])))]];
 
initcon = TimeConstrained[NSolve[
vr0i ==
N[(Sqrt[2] Sqrt[1-μ^2 v0^2] ((q μ^2 ℧ r0 Sqrt[-℧^2+2 r0] v0^2)/((a^2+℧^2+(-2+
r0) r0) Sqrt[a^2+r0^2])+((a^2+a^2 Cos[2 θi]+2 r0^2) dr0y)/(2 (a^2+r0^2))+(Sqrt[-℧^2+
2 r0] (dt0y-a Sin[θi]^2 dφ0y))/Sqrt[a^2+r0^2]))/Sqrt[(a^2+a^2 Cos[2 θi]+2 r0^2)/(a^2+r0^2)]]
&&
vθ0i ==
N[Sqrt[a^2 Cos[θi]^2+r0^2] Sqrt[1-μ^2 v0^2] dθ0y]
&&
vφ0i ==
N[(Sin[θi]^2 Sqrt[1-μ^2 v0^2] (a q μ^2 ℧ r0 Sqrt[a^2+r0^2] v0^2-1/2 a Sqrt[-℧^2+2 r0] (a^2+
a^2 Cos[2 θi]+2 r0^2) dr0y+Sqrt[a^2+r0^2] (a (℧^2-2 r0) dt0y+((a^2+r0^2)^2-a^2 (a^2+℧^2-
2 r0+r0^2) Sin[θi]^2) dφ0y)))/(Sqrt[a^2+r0^2] (a^2 Cos[θi]^2+r0^2) Sqrt[(Sin[θi]^2 ((a^2+
r0^2)^2-a^2 (a^2+℧^2-2 r0+r0^2) Sin[θi]^2))/(a^2 Cos[θi]^2+r0^2)])]
&&
-μ ==
N[-(((a^2+2 r0^2+a^2 Cos[2 θi]) (dr0y)^2)/(2 (a^2+r0^2)))-(2 Sqrt[2 r0-
℧^2] dr0y dt0y)/Sqrt[a^2+r0^2]+(1+(-4 r0+2 ℧^2)/(a^2+2 r0^2+a^2 Cos[2 θi])) (dt0y)^2+
(-r0^2-a^2 Cos[θi]^2) dθ0y^2+(2 a Sqrt[2 r0-℧^2] Sin[θi]^2 dr0y dφ0y)/Sqrt[a^2+
r0^2]+(2 a (2 r0-℧^2) Sin[θi]^2 dt0y dφ0y)/(r0^2+a^2 Cos[θi]^2)+((-(a^2+r0^2)^2 Sin[θi]^2+
a^2 (a^2+(-2+r0) r0+℧^2) Sin[θi]^4) (dφ0y)^2)/(r0^2+a^2 Cos[θi]^2)]
&&
dt0y > 0,
{dθ0y, dr0y, dt0y, dφ0y}, Reals], 3];
                     
DGL = If[NumericQ[dt0y /. initcon[[1]]] == False, {}, {           (* Bewegungsgleichungen *)
 
t''[τ] == Re[1/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3) (8 q ℧ (-a^4 Cos[θ[τ]]^4+r[τ]^4) r'[τ]+
(8 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ]^2)/(Sqrt[-℧^2+
2 r[τ]] Sqrt[a^2+r[τ]^2])+8 q ℧ Sqrt[-℧^2+2 r[τ]] Sqrt[a^2+r[τ]^2] (-a^2 Cos[θ[τ]]^2+
r[τ]^2) t'[τ]+16 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2) r'[τ] t'[τ]+
8 Sqrt[-℧^2+2 r[τ]] (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) Sqrt[a^2+r[τ]^2] t'[τ]^2-
8 a^2 q ℧ r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[2 θ[τ]] θ'[τ]+(8 a^2 Sqrt[-℧^2+
2 r[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/Sqrt[a^2+r[τ]^2]-8 a^2 (℧^2-
2 r[τ]) (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[2 θ[τ]] t'[τ] θ'[τ]+8 r[τ] Sqrt[-℧^2+
2 r[τ]] Sqrt[a^2+r[τ]^2] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-8 a q ℧ Sqrt[-℧^2+
2 r[τ]] Sqrt[a^2+r[τ]^2] (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ]-
16 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 r'[τ] φ'[τ]-
16 a Sqrt[-℧^2+2 r[τ]] (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) Sqrt[a^2+
r[τ]^2] Sin[θ[τ]]^2 t'[τ] φ'[τ]+16 a^3 Cos[θ[τ]] (℧^2-2 r[τ]) (a^2 Cos[θ[τ]]^2+
r[τ]^2) Sin[θ[τ]]^3 θ'[τ] φ'[τ]+Sqrt[-℧^2+2 r[τ]] Sqrt[a^2+r[τ]^2] (a^4+
a^4 Cos[4 θ[τ]] (-1+r[τ])+3 a^4 r[τ]+4 a^2 ℧^2 r[τ]-4 a^2 r[τ]^2+8 a^2 r[τ]^3+
8 r[τ]^5+4 a^2 Cos[2 θ[τ]] r[τ] (a^2-℧^2+r[τ]+2 r[τ]^2)) Sin[θ[τ]]^2 φ'[τ]^2)],
 
t'[0]  == dt0y/.initcon[[1]],
t[0]   == 0,
 
r''[τ] == Re[1/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3) (-8 q ℧ Sqrt[-℧^2+2 r[τ]] Sqrt[a^2+
r[τ]^2] (-a^2 Cos[θ[τ]]^2+r[τ]^2) r'[τ]+(8 a^2 q ℧ Sqrt[-℧^2+2 r[τ]] (-a^2 Cos[θ[τ]]^2+
r[τ]^2) Sin[θ[τ]]^2 r'[τ])/Sqrt[a^2+r[τ]^2]+(4 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 (a^2 Cos[2 θ[τ]] (-1+r[τ])-a^2 (1+r[τ])+2 r[τ] (-℧^2+r[τ])) r'[τ]^2)/(a^2+
r[τ]^2)-4 q ℧ (2 a^2 Cos[θ[τ]]^2-2 r[τ]^2) (a^2+r[τ]^2) (1+(℧^2-
2 r[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)) t'[τ]+(8 a^2 q ℧ (℧^2-2 r[τ]) (a^2 Cos[θ[τ]]^2-
r[τ]^2) Sin[θ[τ]]^2 t'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)-(16 Sqrt[-℧^2+
2 r[τ]] (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+
r[τ]^2) r'[τ] t'[τ])/Sqrt[a^2+r[τ]^2]+8 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+
℧^2-2 r[τ]+r[τ]^2) t'[τ]^2+8 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] r'[τ] θ'[τ]+
8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (a^2+℧^2-2 r[τ]+r[τ]^2) θ'[τ]^2+(8 a q ℧ (℧^2-
2 r[τ]) (a^2 Cos[θ[τ]]^2-r[τ]^2) (a^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ])/(a^2 Cos[θ[τ]]^2+
r[τ]^2)-(4 a q ℧ (2 a^2 Cos[θ[τ]]^2-2 r[τ]^2) (-(a^2+r[τ]^2)^2 Sin[θ[τ]]^2+a^2 (a^2+
℧^2+(-2+r[τ]) r[τ]) Sin[θ[τ]]^4) φ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)-(8 a Sqrt[-℧^2+
2 r[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2 (-1+r[τ])+a^2 Cos[2 θ[τ]] (-1+r[τ])+
2 r[τ] (-℧^2+r[τ]+r[τ]^2)) Sin[θ[τ]]^2 r'[τ] φ'[τ])/Sqrt[a^2+r[τ]^2]-
16 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) (a^2+℧^2-2 r[τ]+
r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+(a^2+℧^2-2 r[τ]+r[τ]^2) (a^4+a^4 Cos[4 θ[τ]] (-1+
r[τ])+3 a^4 r[τ]+4 a^2 ℧^2 r[τ]-4 a^2 r[τ]^2+8 a^2 r[τ]^3+8 r[τ]^5+
4 a^2 Cos[2 θ[τ]] r[τ] (a^2-℧^2+r[τ]+2 r[τ]^2)) Sin[θ[τ]]^2 φ'[τ]^2)],
 
r'[0]  == dr0y/.initcon[[1]],
r[0]   == r0,
 
θ''[τ] == Re[-1/(2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^4) ((2 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3 Sin[θ[τ]] r'[τ]^2)/(a^2+r[τ]^2)-2 a^2 q ℧ (℧^2-2 r[τ]) r[τ] Sin[2 θ[τ]] t'[τ]+
a^2 q ℧ r[τ] (a^2+2 ℧^2+a^2 Cos[2 θ[τ]]-4 r[τ]+2 r[τ]^2) Sin[2 θ[τ]] t'[τ]+
2 a^2 Cos[θ[τ]] (℧^2-2 r[τ]) (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]] t'[τ]^2+
4 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^3 r'[τ] θ'[τ]-2 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3 Sin[θ[τ]] θ'[τ]^2-4 a^3 q ℧ Cos[θ[τ]] (℧^2-2 r[τ]) r[τ] Sin[θ[τ]]^3 φ'[τ]+
4 a q ℧ Cot[θ[τ]] r[τ] (-(a^2+r[τ]^2)^2 Sin[θ[τ]]^2+a^2 (a^2+℧^2+(-2+
r[τ]) r[τ]) Sin[θ[τ]]^4) φ'[τ]+(2 a Sqrt[-℧^2+2 r[τ]] (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3 Sin[2 θ[τ]] r'[τ] φ'[τ])/Sqrt[a^2+r[τ]^2]-2 a (℧^2-2 r[τ]) (a^2+
r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[2 θ[τ]] t'[τ] φ'[τ]+(a^2 Cos[θ[τ]]^2+
r[τ]^2) (2 a^2 Cos[θ[τ]] Sin[θ[τ]]^3 (-(a^2+r[τ]^2)^2+a^2 (a^2+℧^2+(-2+
r[τ]) r[τ]) Sin[θ[τ]]^2)+(a^2 Cos[θ[τ]]^2+r[τ]^2) (4 a^2 Cos[θ[τ]] (a^2+℧^2+
(-2+r[τ]) r[τ]) Sin[θ[τ]]^3-(a^2+r[τ]^2)^2 Sin[2 θ[τ]])) φ'[τ]^2)],
 
θ'[0]  == dθ0y/.initcon[[1]],
θ[0]   == θi,
 
φ''[τ] == Re[1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3) ((4 a q ℧ (-a^4 Cos[θ[τ]]^4+r[τ]^4) r'[τ])/(a^2+
r[τ]^2)+(4 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 r'[τ]^2)/(Sqrt[-℧^2+2 r[τ]] (a^2+r[τ]^2)^(3/2))+(4 a q ℧ Sqrt[-℧^2+
2 r[τ]] (-a^2 Cos[θ[τ]]^2+r[τ]^2) t'[τ])/Sqrt[a^2+r[τ]^2]+(8 a (a^2 Cos[θ[τ]]^2+(℧^2-
r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+r[τ]^2) r'[τ] t'[τ])/(a^2+r[τ]^2)+(4 a Sqrt[-℧^2+
2 r[τ]] (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) t'[τ]^2)/Sqrt[a^2+r[τ]^2]-
8 a q ℧ Cot[θ[τ]] r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2) θ'[τ]+(8 a Cot[θ[τ]] Sqrt[-℧^2+
2 r[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ])/Sqrt[a^2+r[τ]^2]-8 a Cot[θ[τ]] (℧^2-
2 r[τ]) (a^2 Cos[θ[τ]]^2+r[τ]^2) t'[τ] θ'[τ]+(4 a r[τ] Sqrt[-℧^2+2 r[τ]] (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 θ'[τ]^2)/Sqrt[a^2+r[τ]^2]-(4 a^2 q ℧ Sqrt[-℧^2+2 r[τ]] (-a^2 Cos[θ[τ]]^2+
r[τ]^2) Sin[θ[τ]]^2 φ'[τ])/Sqrt[a^2+r[τ]^2]+(8 (a^2 Cos[θ[τ]]^2+
r[τ]^2) (a^4 Cos[θ[τ]]^2 (-1+r[τ])-r[τ] (a^2 (a^2+℧^2-r[τ])+2 a^2 Cot[θ[τ]]^2 (a^2+
r[τ]^2)+Csc[θ[τ]]^2 (-a^4+r[τ]^4))) Sin[θ[τ]]^2 r'[τ] φ'[τ])/(a^2+r[τ]^2)-
(8 a^2 Sqrt[-℧^2+2 r[τ]] (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-
r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ])/Sqrt[a^2+r[τ]^2]-Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+
r[τ]^2) (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) θ'[τ] φ'[τ]+
(4 a Sqrt[-℧^2+2 r[τ]] Sin[θ[τ]]^2 (r[τ] (-a^4+r[τ]^4+a^2 (a^2+℧^2-r[τ]) Sin[θ[τ]]^2)+
Cos[θ[τ]]^2 (2 a^2 r[τ] (a^2+r[τ]^2)-a^4 (-1+r[τ]) Sin[θ[τ]]^2)) φ'[τ]^2)/Sqrt[a^2+r[τ]^2])],
 
φ'[0]  == dφ0y/.initcon[[1]],
φ[0]   == φ0,
 
WhenEvent[Abs[r[τ]] == N[R0]||Abs[r[τ]] > R0 &&
NumericQ[plunge6] == False,
(plunge6=τ) &&
(tt6=t[τ]) && (rt6=r[τ]) && (θt6=θ[τ]) && (φt6=φ[τ]);
"StopIntegration"],
 
WhenEvent[Abs[r[τ]] == N[R1] &&
NumericQ[plunge5] == False,
(plunge5=τ) && (tt5=t[τ]) && (rt5=r[τ]) && (θt5=θ[τ]) && (φt5=φ[τ]);
If[r0>R1, "StopIntegration"]],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && τ < -0.05 &&
NumericQ[plunge0] == False,
(plunge0=τ) &&
(tt0=t[τ]) && (rt0=r[τ]) && (θt0=θ[τ]) && (φt0=φ[τ]) &&
(δt0=t'[τ]) && (δr0=r'[τ]) && (δθ0=θ'[τ]) && (δφ0=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && Re[θ'[τ]]>0 && τ < -0.05 &&
NumericQ[plunge1] == False,
(plunge1=τ) &&
(tt1=t[τ]) && (rt1=r[τ]) && (θt1=θ[τ]) && (φt1=φ[τ]) &&
(δt1=t'[τ]) && (δr1=r'[τ]) && (δθ1=θ'[τ]) && (δφ1=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && Re[θ'[τ]]<0 && τ < -0.05 &&
NumericQ[plunge2] == False,
(plunge2=τ) &&
(tt2=t[τ]) && (rt2=r[τ]) && (θt2=θ[τ]) && (φt2=φ[τ]) &&
(δt2=t'[τ]) && (δr2=r'[τ]) && (δθ2=θ'[τ]) && (δφ2=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && Re[θ'[τ]]>0 &&
τ < plunge1-0.05 &&
NumericQ[plunge3]==False && NumericQ[plunge1] == True,
(plunge3=τ) &&
(tt3=t[τ]) && (rt3=r[τ]) && (θt3=θ[τ]) && (φt3=φ[τ]) &&
(δt3=t'[τ]) && (δr3=r'[τ]) && (δθ3=θ'[τ]) && (δφ3=φ'[τ])],
 
WhenEvent[Mod[180/π Re[θ[τ]], 180] == 90.0 && r[τ]>N[si] &&
r[τ]<N[sr] && Re[θ'[τ]]<0 &&
τ < plunge2-0.05 &&
NumericQ[plunge4] == False && NumericQ[plunge2] == True,
(plunge4=τ) &&
(tt4=t[τ]) && (rt4=r[τ]) && (θt4=θ[τ]) && (φt4=φ[τ]) &&
(δt4=t'[τ]) && (δr4=r'[τ]) && (δθ4=θ'[τ]) && (δφ4=φ'[τ])]
 
}];
                                                                            (* Integrator *)
sol = TimeConstrained[If[NumericQ[dt0y /. initcon[[1]]] == False, {},
NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
Method-> mta,
MaxSteps-> Infinity,
InterpolationOrder-> All,
StepMonitor :> (laststep=plunge7; plunge7=τ;
stepsize=plunge7-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<1*^-4, 0, 1])}
]], 10];
 
ft5h=If[NumericQ[plunge5], {φt5-π, θt5+π/2}, {0, -π/2}];
ft5b=If[NumericQ[plunge6], {φt6-π, θt6+π/2}, {0, -π/2}];
ft5f={0, Min[fmax, ς shF]};
ft5x=If[NumericQ[plunge6], {0, 1}, {0, 0}];
ft5y=If[NumericQ[plunge5], {0, 1}, {0, 0}];
 
ft0s=If[NumericQ[plunge0], {φt0, rt0}, {π, sr}];
ft1s=If[NumericQ[plunge1], {φt1, rt1}, {π, sr}];
ft2s=If[NumericQ[plunge2], {φt2, rt2}, {π, sr}];
ft3s=If[NumericQ[plunge3], {φt3, rt3}, {π, sr}];
ft4s=If[NumericQ[plunge4], {φt4, rt4}, {π, sr}];
 
ft0v=If[NumericQ[plunge0], {-dφv[rt0, tt0+t0+rT], 0}, {0, 0}];
ft1v=If[NumericQ[plunge1], {-dφv[rt1, tt1+t0+rT], 0}, {0, 0}];
ft2v=If[NumericQ[plunge2], {-dφv[rt2, tt2+t0+rT], 0}, {0, 0}];
ft3v=If[NumericQ[plunge3], {-dφv[rt3, tt3+t0+rT], 0}, {0, 0}];
ft4v=If[NumericQ[plunge4], {-dφv[rt4, tt4+t0+rT], 0}, {0, 0}];
ft5v=If[NumericQ[plunge5], {-(tt5+t0+rT) ωR1, 0}, {0, 0}];
 
ft0f=If[NumericQ[plunge0], {0, Min[fmax, shF shf[rt0, δt0, δr0, δθ0, δφ0]]}, {0, 0}];
ft1f=If[NumericQ[plunge1], {0, Min[fmax, shF shf[rt1, δt1, δr1, δθ1, δφ1]]}, {0, 0}];
ft2f=If[NumericQ[plunge2], {0, Min[fmax, shF shf[rt2, δt2, δr2, δθ2, δφ2]]}, {0, 0}];
ft3f=If[NumericQ[plunge3], {0, Min[fmax, shF shf[rt3, δt3, δr3, δθ3, δφ3]]}, {0, 0}];
ft4f=If[NumericQ[plunge4], {0, Min[fmax, shF shf[rt4, δt4, δr4, δθ4, δφ4]]}, {0, 0}];
{
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0s+ft0v, ft1s+ft1v, ft2s+ft2v, ft3s+ft3v, ft4s+ft4v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h+ft5v, ft5b, ft5f, ft5x, ft5y
}]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 10) Memoryfunktion ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
mem : raytrace[{Ф_, ϑ_}] := mem = Quiet[Re[raytracer[{Ф, ϑ}]]];
 
ray01[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[01]];
ray02[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[02]];
ray03[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[03]];
ray04[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[04]];
ray05[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[05]];
 
ray06[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[06]];
ray07[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[07]];
ray08[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[08]];
ray09[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[09]];
ray10[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[10]];
 
ray11[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[11]];
ray12[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[12]];
ray13[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[13]];
ray14[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[14]];
ray15[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[15]];
 
ray16[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[16]];
ray17[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[17]];
ray18[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[18]];
ray19[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[19]];
ray20[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[20]];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 11) Proportionen ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
width1  = ImageDimensions[pic1][[1]]; height1 = ImageDimensions[pic1][[2]];
width2  = ImageDimensions[pic2][[1]]; height2 = ImageDimensions[pic2][[2]];
width3  = ImageDimensions[pic3][[1]]; height3 = ImageDimensions[pic3][[2]];
width4  = ImageDimensions[pic4][[1]]; height4 = ImageDimensions[pic4][[2]];
width5  = ImageDimensions[pic5][[1]]; height5 = ImageDimensions[pic5][[2]];
width6  = ImageDimensions[pic7][[1]]; height6 = ImageDimensions[pic7][[2]];
hzoom   = If[breite>2 hoehe, 1/zoom, 1/zoom/2/hoehe*breite];
vzoom   = If[breite>2 hoehe, 1/zoom*2 hoehe/breite, 1/zoom];
 
"approximate time remaining" -> 1.2 (AbsoluteTiming[Do[raytracer[{
RandomReal[{-π, π}]/zoom, RandomReal[{-π/2, π/2}]/zoom
}], {ü, 1, 50}]][[1]])/50*hoehe*breite/kernels/60 "minutes"
 
FOV->{360.0 hzoom "degree", 180.0 vzoom "degree"}
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 12) Output ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
img = ParallelTable[{
 
(* 1 Hintergrundpanorama *)
Quiet@ImageTransformation[pct1, ray17, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-2π/width1},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+vvs/vzoom+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 2 Sphäre *)
Quiet@ImageTransformation[pct2, ray16, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width2},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+π/kernels/grain+vvs/vzoom} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 3 Scheibe Textur komplett *)
Quiet@ImageTransformation[pic3, ray01, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 4 Scheibe Textur Front *)
Quiet@ImageTransformation[pic3, ray02, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 5 Scheibe Textur Echo 1 *)
Quiet@ImageTransformation[pic3, ray03, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 6 Scheibe Textur Echo 2 *)
Quiet@ImageTransformation[pic3, ray04, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 7 Scheibe Textur Echo 3 *)
Quiet@ImageTransformation[pic3, ray05, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width3},
{si, sr+(sr-si)/height3}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 8 Scheibe Geometrie still komplett *)
Quiet@ImageTransformation[pic4, ray01, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 9 Scheibe Geometrie still Front *)
Quiet@ImageTransformation[pic4, ray02, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 10 Scheibe Geometrie still Echo 1 *)
Quiet@ImageTransformation[pic4, ray03, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 11 Scheibe Geometrie still Echo 2 *)
Quiet@ImageTransformation[pic4, ray04, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 12 Scheibe Geometrie still Echo 3 *)
Quiet@ImageTransformation[pic4, ray05, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 13 Scheibe Geometrie rotierend komplett *)
Quiet@ImageTransformation[pic4, ray06, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 14 Scheibe Geometrie rotierend Front *)
Quiet@ImageTransformation[pic4, ray07, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 15 Scheibe Geometrie rotierend Echo 1 *)
Quiet@ImageTransformation[pic4, ray08, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 16 Scheibe Geometrie rotierend Echo 2 *)
Quiet@ImageTransformation[pic4, ray09, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 17 Scheibe Geometrie rotierend Echo 3 *)
Quiet@ImageTransformation[pic4, ray10, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width4},
{si, sr}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 18 Frequenzverschiebung komplett SW *)
Quiet@ImageTransformation[pic5, ray11, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-1, 1},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 19 Frequenzverschiebung Front SW *)
Quiet@ImageTransformation[pic5, ray12, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 20 Frequenzverschiebung Echo 1 SW *)
Quiet@ImageTransformation[pic5, ray13, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 21 Frequenzverschiebung Echo 2 SW *)
Quiet@ImageTransformation[pic5, ray14, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 22 Frequenzverschiebung Echo 3 SW *)
Quiet@ImageTransformation[pic5, ray15, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 23 Frequenzverschiebung komplett Farbe *)
Quiet@ImageTransformation[pic6, ray11, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 24 Frequenzverschiebung Front Farbe *)
Quiet@ImageTransformation[pic6, ray12, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 25 Frequenzverschiebung Echo 1 Farbe *)
Quiet@ImageTransformation[pic6, ray13, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 26 Frequenzverschiebung Echo 2 Farbe *)
Quiet@ImageTransformation[pic6, ray14, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 27 Frequenzverschiebung Echo 3 Farbe *)
Quiet@ImageTransformation[pic6, ray15, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],
 
(* 28 Hintergrundpanorama *)
ImageTransformation[pct3, ray17, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-2π/width6},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+vvs/vzoom+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 29 Sphäre *)
ImageTransformation[pct3, ray16, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-2π/width6},
{-π/2, 3π/2}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+x+vvs/vzoom, -π/2+x+π/kernels/grain+vvs/vzoom} vzoom
},
Resampling->rsp, Padding->"Periodic"],
 
(* 30 Scheibe Textur komplett *)
ImageTransformation[pic7, ray01, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π-2π/width6},
{si, sr+(sr-si)/height6}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Periodic"],

(* 31 Frequenzverschiebung Hintergrund Farbe *)
Quiet@ImageTransformation[pic6, ray18, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],

(* 32 Frequenzverschiebung Hintergrund SW *)
Quiet@ImageTransformation[pic5, ray18, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{fmin, fmax}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],

(* 33 Kontur BH *)
Quiet@ImageTransformation[pic5, ray19, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 1}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"],

(* 34 Kontur IH *)
Quiet@ImageTransformation[pic5, ray20, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 1}
},
PlotRange->{
{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
{-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom
},
Resampling->rsp, Padding->"Fixed"]
 
}, {x, 0, π-π/kernels/grain, π/kernels/grain}];
 
image[num_] := ImageAssemble[{Table[{img[[x, num]]}, {x, kernels grain, 1, -1}]}];
fi[x_] := ColorNegate[x];
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 13) Composite |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ParallelTable[image[n], {n, 1, 34, 1}]
 
ds1 = image[8]
ds2 = ImageMultiply[fi@image[09], fi@image[10]];
ds3 = ImageMultiply[fi@image[11], fi@image[12]];
ds4 = fi@ImageMultiply[ds2, ds3]
 
ds5 = image[13]
ds6 = ImageMultiply[fi@image[14], fi@image[15]];
ds7 = ImageMultiply[fi@image[16], fi@image[17]];
ds8 = fi@ImageMultiply[ds6, ds7]

bf1 = image[31]
bf2 = image[32]
bf3 = image[33]
bf4 = image[34]
 
bl1 = image[23]
bl2 = ImageMultiply[fi@image[24], fi@image[25]];
bl3 = ImageMultiply[fi@image[26], fi@image[27]];
bl4 = fi@ImageMultiply[bl2, bl3]
 
gr1 = image[18]
gr2 = ImageMultiply[fi@image[19], fi@image[20]];
gr3 = ImageMultiply[fi@image[21], fi@image[22]];
gr4 = fi@ImageMultiply[gr2, gr3]
 
sea = ImageMultiply[image[04], image[19]];
seb = ImageMultiply[image[05], image[20]];
sec = ImageMultiply[image[06], image[21]];
sed = ImageMultiply[image[07], image[22]];
 
se1 = image[03]
se2 = ImageMultiply[fi@image[04], fi@image[05]];
se3 = ImageMultiply[fi@image[06], fi@image[07]];
se4 = fi@ImageMultiply[se2, se3]
se5 = ImageMultiply[fi@sea, fi@seb];
se6 = ImageMultiply[fi@sec, fi@sed];
se7 = fi@ImageMultiply[se5, se6]
se8 = fi@ImageMultiply[fi@se7, fi@se7]
se9 = fi@ImageMultiply[fi@se8, fi@se7]
 
bkg = image[1]
hrz = image[2]
 
cb1 = image[28]
cb2 = image[29]
cb3 = image[30]
 
Quit[]