











ungelinste (oben) und gravitationsgelinste Ansicht (unten) eines rotierenden und geladenen Kerr Newman SL mit Akkretionsscheibe, FOV=90°×45°:



Codes (Syntax: Wolfram)


Code: Alles auswählen
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 - 29.07.2019 | Version 8A | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
ClearAll["Global`*"];
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
kernels = 6; (* Parallelisierung *)
grain = 5; (* Subparallelisierung auf kernels*grain Streifen *)
rsp = "Nearest"; (* Resampling *)
breite = 120; (* Zielabmessungen in Pixeln *)
hoehe = 120; (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom = 15; (* doppelter Zoom ergibt halben Sichtwinkel *)
LaunchKernels[kernels]
wp = MachinePrecision; (* Genauigkeit *)
st = 0.02; (* Auflösung des Gradienten *)
pic1 = Import["http://yukterez.net/mw/2/flip70.png"]; (* Hintergrundpanorama laden *)
pic2 = Import["https://i.stack.imgur.com/5ARxo.jpg"]; (* Sphärenoberfläche laden *)
pic3 = Import["http://yukterez.net/mw/akkretionsscheibe.jpg"]; (* Scheibentextur laden *)
pic4 = Import["http://yukterez.net/mw/disk.png"]; (* Scheibengeometrie laden *)
pic5 = Import["http://yukterez.net/mw/gradient1.png"]; (* Helligkeitsgradient laden *)
pic6 = Import["http://yukterez.net/mw/gradient2.png"]; (* Farbgradient laden *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
r0 = 100; (* Radialkoordinate des Beobachters *)
R0 = 500; (* Radius des umspannenden Kugelschalenpanoramas *)
R1 = Max[1, Re[1.01 rA]]; (* Sphäre *)
si = isco; (* Akkretionsscheibe Innenradius *)
sr = 7; (* Akkretionsscheibe Außenradius *)
θ0 = 70 π/180; (* Breitengrad *)
φ0 = 0; (* Längengrad *)
tmax =-3 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 = 0; (* 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 *)
hvs = ArcSin[vφ]; (* horizontaler Versatz in Radianten *)
vvs = ArcSin[vϑ]; (* vertikaler Versatz in Radianten *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 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}]}];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 3) Funktionen |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
gtt = (2r0-℧^2)/Σ-1;
grr = Σ/Δ;
gθθ = Σ;
gφφ = Χ/Σ Sin[θ0]^2;
gtφ =-a (2r0-℧^2) Sin[θ0]^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];
ς = Sqrt[Χ/Δ/Σ];
j[v_] := Sqrt[1-v^2];
Ы[rs_] := Sqrt[Χs[rs]/Σs[rs]];
ωs[rs_] := (a (2 rs - ℧^2))/Χs[rs];
ε[rs_] := Sqrt[Δs[rs] Σs[rs]/Χs[rs]]/j[vt]+Lz[rs] ωs[rs];
Lz[rs_] := vt Ы[rs]/j[vt];
nq[x_] := If[NumericQ[x], If[Element[x, Reals], x, 0], 0];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 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_, δθ_, δφ_] := ς/(1-vs[rt] vφt[rt, δt, δr, δθ, δφ])/δt;
dφ[rt_, tt_] := (tt+r0) φs[rt]/ts[rt];
vθ =-vϑ;
θs = π/2;
θi =-θ0+π;
rA = 1+Sqrt[1-a^2-℧^2];
rE = 1+Sqrt[1-℧^2];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 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" -> [email protected], "r ergosphere" -> [email protected], "r isco" -> [email protected],
"r photon pro" -> [email protected]Min[rp], "r photon ret" -> [email protected][rp], "r disk" -> [email protected],
"r observer" -> [email protected], "θ observer" -> [email protected]θ0 180/π}
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 6) Geschwindigkeit und Zeitdilatation auf der Scheibe |||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
red[rs_] := Quiet[Reduce[
dt == (Lz[rs] (-a (a^2+rs^2)+Δs[rs] κs[rs])+ε[rs] ((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] (-a^2+Δs[rs])+ε[rs] (a (a^2+rs^2)-Δs[rs] κs[rs]))/(Δs[rs] Σs[rs])
&&
vt > 0,
{vt, dΦ, dt},
Reals]];
vs = Interpolation[ParallelTable[{rr, If[[email protected][red[rr][[1, 2]]],
red[rr][[1, 2]], 0]}, {rr, 0, sr+st, st}]];
φs = Interpolation[ParallelTable[{rr, If[[email protected][red[rr][[2, 2]]],
red[rr][[2, 2]], 0]}, {rr, 0, sr+st, st}]];
ts = Interpolation[ParallelTable[{rr, If[[email protected][red[rr][[3, 2]]],
red[rr][[3, 2]], 0]}, {rr, 0, sr+st, st}]];
plot[func_, label_] := Plot[func, {rr, rP, sr},
GridLines -> {{nq[Min[rp]], nq[Max[rp]], nq[rA], nq[si], nq[isco], nq[rE], nq[sr]}, {}},
Frame -> True, ImagePadding -> {{40, 1}, {12, 1}}, ImageSize -> 340,
PlotLabel -> label, PlotRange->{{0, sr}, Automatic}]
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)]);
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,
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, mta, ft, fτ, varb,
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0v, ft1v, ft2v, ft3v, ft4v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b},
vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
(* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
vr0a = vw[[3]] Sqrt[grr];
vφ0a = vw[[2]] Sqrt[gφφ]/r0/Sin[θi];
vθia = vw[[1]] Sqrt[gθθ]/r0;
(* 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]];
DGL = { (* Kerr Newman Bewegungsgleichungen *)
t''[τ]==-(((r'[τ] ((a^2+r[τ]^2) (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+
2 (-℧^2+r[τ]) t'[τ]))+a (2 a^4 Cos[θ[τ]]^2+a^2 ℧^2 (3+Cos[2 θ[τ]]) r[τ]-
a^2 (3+Cos[2 θ[τ]]) r[τ]^2+4 ℧^2 r[τ]^3-6 r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))/(a^2+℧^2+(-2+
r[τ]) r[τ])+a^2 θ'[τ] (Sin[2 θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])-2 a Cos[θ[τ]] (℧^2-
2 r[τ]) Sin[θ[τ]]^3 φ'[τ]))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^2),
t'[0]==-((a (2 r0-℧^2) Sin[θi]^2 (vφ0i (-2 r0+r0^2+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+
r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-
2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)]+
(a vφ0i (2 r0-℧^2) Sin[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+
a^2 Cos[θi]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2))))/((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θi]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θi]^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[θi]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θi]^2)+
(a^2 (-2 r0+℧^2)^2 Sin[θi]^4 (vφ0i (-2 r0+r0^2+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+r0^2)^2-
a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-2 r0+
r0^2+℧^2) (r0^2+a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)]+
(a vφ0i (2 r0-℧^2) Sin[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+
a^2 Cos[θi]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)))^2)/((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θi]^2)^2)-((2 r0-r0^2-℧^2-a^2 Cos[θi]^2) Sin[θi]^2 ((a^2+r0^2)^2-
a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2) (vφ0i (-2 r0+r0^2+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+
r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-
2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)]+
(a vφ0i (2 r0-℧^2) Sin[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+
a^2 Cos[θi]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)))^2)/((a^2-2 r0+r0^2+
℧^2) (r0^2+a^2 Cos[θi]^2)^2))/((a^2-2 r0+r0^2+℧^2) (-2 r0+r0^2+℧^2+a^2 Cos[θi]^2)^2)),
t[0]==0,
r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
(a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θi]^2)/(a^2+(-2+r0) r0+℧^2)],
r[0]==r0,
θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
10 a^4 r[τ]+11 a^4 r[τ]^2-8 a^2 ℧^2 r[τ]^2+16 a^2 r[τ]^3+16 a^2 r[τ]^4+8 r[τ]^6+
a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2],
θ[0]==θi,
φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),
φ'[0]==(vφ0i ((-2+r0) r0+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+
a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θi]^2)]+(a vφ0i (2 r0-
℧^2) Sin[θi])/((r0^2+a^2 Cos[θi]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θi]^2)),
φ[0]==φ0,
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr &&
NumericQ[plunge0]==False,
(plunge0=τ) &&
(tt0=t[τ]) && (rt0=r[τ]) && (θt0=θ[τ]) && (φt0=φ[τ]) &&
(δt0=t'[τ]) && (δr0=r'[τ]) && (δθ0=θ'[τ]) && (δφ0=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0 &&
NumericQ[plunge1]==False,
(plunge1=τ) &&
(tt1=t[τ]) && (rt1=r[τ]) && (θt1=θ[τ]) && (φt1=φ[τ]) &&
(δt1=t'[τ]) && (δr1=r'[τ]) && (δθ1=θ'[τ]) && (δφ1=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0 &&
NumericQ[plunge2]==False,
(plunge2=τ) &&
(tt2=t[τ]) && (rt2=r[τ]) && (θt2=θ[τ]) && (φt2=φ[τ]) &&
(δt2=t'[τ]) && (δr2=r'[τ]) && (δθ2=θ'[τ]) && (δφ2=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0 && τ<plunge1-0.1 &&
NumericQ[plunge3]==False,
(plunge3=τ) &&
(tt3=t[τ]) && (rt3=r[τ]) && (θt3=θ[τ]) && (φt3=φ[τ]) &&
(δt3=t'[τ]) && (δr3=r'[τ]) && (δθ3=θ'[τ]) && (δφ3=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0 && τ<plunge2-0.1 &&
NumericQ[plunge4]==False,
(plunge4=τ) &&
(tt4=t[τ]) && (rt4=r[τ]) && (θt4=θ[τ]) && (φt4=φ[τ]) &&
(δt4=t'[τ]) && (δr4=r'[τ]) && (δθ4=θ'[τ]) && (δφ4=φ'[τ])],
WhenEvent[r[τ]==R1||r[τ]<R1 &&
NumericQ[plunge5]==False,
(plunge5=τ) && (tt5=t[τ]) && (rt5=r[τ]) && (θt5=θ[τ]) && (φt5=φ[τ])],
WhenEvent[r[τ]==R0||r[τ]>R0 &&
NumericQ[plunge6]==False,
(plunge6=τ) &&
(tt6=t[τ]) && (rt6=r[τ]) && (θt6=θ[τ]) && (φt6=φ[τ]);
"StopIntegration"]
};
(* Integrator *)
sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
InterpolationOrder-> All];
ft0s=If[NumericQ[plunge0], If[rt0>sr, {π, sr},
If[rt0<si, {π, sr}, {φt0, rt0}]], {π, sr}];
ft1s=If[NumericQ[plunge1], If[rt1>sr, {π, sr},
If[rt1<si, {π, sr}, {φt1, rt1}]], {π, sr}];
ft2s=If[NumericQ[plunge2], If[rt2>sr, {π, sr},
If[rt2<si, {π, sr}, {φt2, rt2}]], {π, sr}];
ft3s=If[NumericQ[plunge3], If[rt3>sr, {π, sr},
If[rt3<si, {π, sr}, {φt3, rt3}]], {π, sr}];
ft4s=If[NumericQ[plunge4], If[rt4>sr, {π, sr},
If[rt4<si, {π, sr}, {φt4, rt4}]], {π, sr}];
ft0v=If[NumericQ[plunge0], If[rt0>sr, {0, 0},
If[rt0<si, {0, 0}, {-dφ[rt0, tt0], 0}]], {0, 0}];
ft1v=If[NumericQ[plunge1], If[rt1>sr, {0, 0},
If[rt1<si, {0, 0}, {-dφ[rt1, tt1], 0}]], {0, 0}];
ft2v=If[NumericQ[plunge2], If[rt2>sr, {0, 0},
If[rt2<si, {0, 0}, {-dφ[rt2, tt2], 0}]], {0, 0}];
ft3v=If[NumericQ[plunge3], If[rt3>sr, {0, 0},
If[rt3<si, {0, 0}, {-dφ[rt3, tt3], 0}]], {0, 0}];
ft4v=If[NumericQ[plunge4], If[rt4>sr, {0, 0},
If[rt4<si, {0, 0}, {-dφ[rt4, tt4], 0}]], {0, 0}];
ft0f=If[NumericQ[plunge0], If[rt0>sr, {0, 0},
If[rt0<si, {0, 0}, {0, Min[2, shf[rt0, δt0, δr0, δθ0, δφ0]]}]], {0, 0}];
ft1f=If[NumericQ[plunge1], If[rt1>sr, {0, 0},
If[rt1<si, {0, 0}, {0, Min[2, shf[rt1, δt1, δr1, δθ1, δφ1]]}]], {0, 0}];
ft2f=If[NumericQ[plunge2], If[rt2>sr, {0, 0},
If[rt2<si, {0, 0}, {0, Min[2, shf[rt2, δt2, δr2, δθ2, δφ2]]}]], {0, 0}];
ft3f=If[NumericQ[plunge3], If[rt3>sr, {0, 0},
If[rt3<si, {0, 0}, {0, Min[2, shf[rt3, δt3, δr3, δθ3, δφ3]]}]], {0, 0}];
ft4f=If[NumericQ[plunge4], If[rt4>sr, {0, 0},
If[rt4<si, {0, 0}, {0, Min[2, shf[rt4, δt4, δr4, δθ4, δφ4]]}]], {0, 0}];
ft5h=If[NumericQ[plunge5], {φt5-(tt5+r0) ωs[R1], θ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}];
{
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0s+ft0v, ft1s+ft1v, ft2s+ft2v, ft3s+ft3v, ft4s+ft4v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b
}]];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 10) Memoryfunktion ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
mem : raytrace[{Ф_, ϑ_}] := mem = 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]];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 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]];
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 *)
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 *)
ImageTransformation[pct2, ray16, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
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 *)
ImageTransformation[pic5, ray11, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-1, 1},
{0, 2}
},
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 *)
ImageTransformation[pic5, ray12, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic5, ray13, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic5, ray14, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic5, ray15, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic6, ray11, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic6, ray12, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic6, ray13, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic6, ray14, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 *)
ImageTransformation[pic6, ray15, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{0, 2π},
{0, 2}
},
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 |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
Table[image[n], {n, 1, 27, 1}]
ds1 = image[8]
ds2 = ImageMultiply[[email protected][09], [email protected][10]];
ds3 = ImageMultiply[[email protected][11], [email protected][12]];
ds4 = [email protected][ds2, ds3]
ds5 = image[13]
ds6 = ImageMultiply[[email protected][14], [email protected][15]];
ds7 = ImageMultiply[[email protected][16], [email protected][17]];
ds8 = [email protected][ds6, ds7]
bl1 = image[23]
bl2 = ImageMultiply[[email protected][24], [email protected][25]];
bl3 = ImageMultiply[[email protected][26], [email protected][27]];
bl4 = [email protected][bl2, bl3]
gr1 = image[18]
gr2 = ImageMultiply[[email protected][19], [email protected][20]];
gr3 = ImageMultiply[[email protected][21], [email protected][22]];
gr4 = [email protected][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[[email protected][04], [email protected][05]];
se3 = ImageMultiply[[email protected][06], [email protected][07]];
se4 = [email protected][se2, se3]
se5 = ImageMultiply[[email protected], [email protected]];
se6 = ImageMultiply[[email protected], [email protected]];
se7 = [email protected][se5, se6]
se8 = [email protected][[email protected], [email protected]]
se9 = [email protected][[email protected], [email protected]]
bkg = image[1]
hrz = image[2]
Output (in diesem Beispiel mit a=℧=0.3, θ=85°, r0=50, ri=isco=4.826, ra=10, v=0, zoom=6):


☑ Kontrollmodul Minkowski (all in one)
Code: Alles auswählen
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* | yukterez.net | Minkowski Raytracer | 29.07.2019 | Version 8M | Simon Tyran, Vienna | *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
kernels = 6; (* Parallelisierung *)
grain = 5; (* Subparallelisierung auf kernels*grain Streifen *)
rsp = "Nearest"; (* Resampling *)
breite = 120; (* Zielabmessungen in Pixeln *)
hoehe = 120; (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom = 15; (* doppelter Zoom ergibt halben Sichtwinkel *)
LaunchKernels[kernels]
wp = MachinePrecision; (* Genauigkeit *)
pic1 = Import["http://yukterez.net/mw/2/flip70.png"]; (* Hintergrundpanorama laden *)
pic2 = Import["https://i.stack.imgur.com/5ARxo.jpg"]; (* Sphärenoberfläche laden *)
pic3 = Import["http://yukterez.net/mw/akkretionsscheibe.jpg"]; (* Scheibentextur laden *)
pic4 = Import["http://yukterez.net/mw/disk.png"]; (* Scheibengeometrie laden *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
r0 = 100; (* Radialkoordinate des Beobachters *)
R0 = 500; (* Radius des umspannenden Kugelschalenpanoramas *)
R1 = 3; (* Sphäre *)
si = 3; (* Akkretionsscheibe Innenradius *)
sr = 7; (* Akkretionsscheibe Außenradius *)
θ0 = 70 π/180; (* Breitengrad *)
φ0 = 0; (* Längengrad *)
tmax =-3 R0; (* zeitlicher Integrationsbereich *)
vr = 0; (* Radiale Geschwindigkeit des Beobachters *)
vφ = 0; (* Azimutale Geschwindigkeit des Beobachters: 0 für ZAMO, -й0 für stationär *)
hvs = ArcSin[vφ]; (* Aberrationswinkel *)
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}}, 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}}, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct2 = ImageAssemble[{Table[{pcr2[[x]]}, {x, 2 kernels, 1, -1}]}];
a = 0; ℧ = 0; v0 = 1; q = 0; st = 0.1; vϑ = 0; vvs = ArcSin[vϑ];
vθ =-vϑ;
θs = π/2;
θi =-θ0+π;
θ1 = θi;
dφ[rt_, tt_] := 0;
shf[rt_, δt_, δr_, δθ_, δφ_] := 1;
ωφ = -vφ Csc[θ1]/r0;
ωθ = -vϑ/r0;
gtt = -1;
grr = +1;
gθθ = +r0^2;
gφφ = +r0^2 Sin[θ1]^2;
gtφ = +0;
ς = 1;
j[v_] := Sqrt[1-v^2];
nq[x_] := If[NumericQ[x], If[Element[x, Reals], x, 0], 0];
U={+vr, +vθ, +vφ};
γ=1/Sqrt[1-Norm[U]^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[ψ]};
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,
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, mta, ft, fτ, varb,
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0v, ft1v, ft2v, ft3v, ft4v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b},
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]];
DGL = { (* Kerr Newman Bewegungsgleichungen *)
t''[τ]==0,
t'[0]==1,
t[0]==0,
r''[τ]==r[τ](θ'[τ]^2+Sin[θ[τ]]^2 φ'[τ]^2),
r'[0]==vr0i,
r[0]==r0,
θ''[τ]==Sin[θ[τ]] Cos[θ[τ]] φ'[τ]^2-2 θ'[τ] r'[τ]/r[τ],
θ'[0]==vθ0i/r0,
θ[0]==θi,
φ''[τ]==-2 φ'[τ] (r'[τ]+r[τ] θ'[τ] Cot[θ[τ]])/r[τ],
φ'[0]==vφ0i Csc[θ1]/r0,
φ[0]==φ0,
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr &&
NumericQ[plunge0]==False,
(plunge0=τ) &&
(tt0=t[τ]) && (rt0=r[τ]) && (θt0=θ[τ]) && (φt0=φ[τ]) &&
(δt0=t'[τ]) && (δr0=r'[τ]) && (δθ0=θ'[τ]) && (δφ0=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0 &&
NumericQ[plunge1]==False,
(plunge1=τ) &&
(tt1=t[τ]) && (rt1=r[τ]) && (θt1=θ[τ]) && (φt1=φ[τ]) &&
(δt1=t'[τ]) && (δr1=r'[τ]) && (δθ1=θ'[τ]) && (δφ1=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0 &&
NumericQ[plunge2]==False,
(plunge2=τ) &&
(tt2=t[τ]) && (rt2=r[τ]) && (θt2=θ[τ]) && (φt2=φ[τ]) &&
(δt2=t'[τ]) && (δr2=r'[τ]) && (δθ2=θ'[τ]) && (δφ2=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0 && τ<plunge1-0.1 &&
NumericQ[plunge3]==False,
(plunge3=τ) &&
(tt3=t[τ]) && (rt3=r[τ]) && (θt3=θ[τ]) && (φt3=φ[τ]) &&
(δt3=t'[τ]) && (δr3=r'[τ]) && (δθ3=θ'[τ]) && (δφ3=φ'[τ])],
WhenEvent[Mod[θ[τ], π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0 && τ<plunge2-0.1 &&
NumericQ[plunge4]==False,
(plunge4=τ) &&
(tt4=t[τ]) && (rt4=r[τ]) && (θt4=θ[τ]) && (φt4=φ[τ]) &&
(δt4=t'[τ]) && (δr4=r'[τ]) && (δθ4=θ'[τ]) && (δφ4=φ'[τ])],
WhenEvent[r[τ]==R1||r[τ]<R1 &&
NumericQ[plunge5]==False,
(plunge5=τ) && (tt5=t[τ]) && (rt5=r[τ]) && (θt5=θ[τ]) && (φt5=φ[τ]);
"StopIntegration"],
WhenEvent[r[τ]==R0||r[τ]>R0 &&
NumericQ[plunge6]==False,
(plunge6=τ) &&
(tt6=t[τ]) && (rt6=r[τ]) && (θt6=θ[τ]) && (φt6=φ[τ]);
"StopIntegration"]
};
(* Integrator *)
sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
InterpolationOrder-> All];
ft0s=If[NumericQ[plunge0], If[rt0>sr, {π, sr},
If[rt0<si, {π, sr}, {φt0+(tt0+r0) ωφ, rt0}]], {π, sr}];
ft1s=If[NumericQ[plunge1], If[rt1>sr, {π, sr},
If[rt1<si, {π, sr}, {φt1+(tt1+r0) ωφ, rt1}]], {π, sr}];
ft2s=If[NumericQ[plunge2], If[rt2>sr, {π, sr},
If[rt2<si, {π, sr}, {φt2+(tt2+r0) ωφ, rt2}]], {π, sr}];
ft3s=If[NumericQ[plunge3], If[rt3>sr, {π, sr},
If[rt3<si, {π, sr}, {φt3+(tt3+r0) ωφ, rt3}]], {π, sr}];
ft4s=If[NumericQ[plunge4], If[rt4>sr, {π, sr},
If[rt4<si, {π, sr}, {φt4+(tt4+r0) ωφ, rt4}]], {π, sr}];
ft0v=If[NumericQ[plunge0], If[rt0>sr, {0, 0},
If[rt0<si, {0, 0}, {-dφ[rt0, tt0], 0}]], {0, 0}];
ft1v=If[NumericQ[plunge1], If[rt1>sr, {0, 0},
If[rt1<si, {0, 0}, {-dφ[rt1, tt1], 0}]], {0, 0}];
ft2v=If[NumericQ[plunge2], If[rt2>sr, {0, 0},
If[rt2<si, {0, 0}, {-dφ[rt2, tt2], 0}]], {0, 0}];
ft3v=If[NumericQ[plunge3], If[rt3>sr, {0, 0},
If[rt3<si, {0, 0}, {-dφ[rt3, tt3], 0}]], {0, 0}];
ft4v=If[NumericQ[plunge4], If[rt4>sr, {0, 0},
If[rt4<si, {0, 0}, {-dφ[rt4, tt4], 0}]], {0, 0}];
ft0f=If[NumericQ[plunge0], If[rt0>sr, {0, 0},
If[rt0<si, {0, 0}, {0, Min[2, shf[rt0, δt0, δr0, δθ0, δφ0]]}]], {0, 0}];
ft1f=If[NumericQ[plunge1], If[rt1>sr, {0, 0},
If[rt1<si, {0, 0}, {0, Min[2, shf[rt1, δt1, δr1, δθ1, δφ1]]}]], {0, 0}];
ft2f=If[NumericQ[plunge2], If[rt2>sr, {0, 0},
If[rt2<si, {0, 0}, {0, Min[2, shf[rt2, δt2, δr2, δθ2, δφ2]]}]], {0, 0}];
ft3f=If[NumericQ[plunge3], If[rt3>sr, {0, 0},
If[rt3<si, {0, 0}, {0, Min[2, shf[rt3, δt3, δr3, δθ3, δφ3]]}]], {0, 0}];
ft4f=If[NumericQ[plunge4], If[rt4>sr, {0, 0},
If[rt4<si, {0, 0}, {0, Min[2, shf[rt4, δt4, δr4, δθ4, δφ4]]}]], {0, 0}];
ft5h=If[NumericQ[plunge5], {φt5+(tt5+r0) ωφ, θt5+π/2+(tt5+r0) ωθ}, {0, -π/2}];
ft5b=If[NumericQ[plunge6], If[rt6<If[Element[R1, Reals], R1, 0], {0, -π/2},
If[rt6>4 R0, {0, -π/2}, {φt6-π, θt6+π/2}]], {0, -π/2}];
{
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0s+ft0v, ft1s+ft1v, ft2s+ft2v, ft3s+ft3v, ft4s+ft4v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b
}]];
mem : raytrace[{Ф_, ϑ_}] := mem = 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]];
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]];
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"}
img = ParallelTable[{
(* 1 Hintergrundpanorama *)
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 *)
ImageTransformation[pct2, ray16, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{
{-π, π-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 *)
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 Geometrie *)
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"]
}, {x, 0, π-π/kernels/grain, π/kernels/grain}];
image[num_] := ImageAssemble[{Table[{img[[x, num]]}, {x, kernels grain, 1, -1}]}];
fi[x_] := ColorNegate[x];
Grid[{Table[image[n], {n, 1, 4, 1}]}]
Output (in diesem Beispiel mit unbewegter Scheibe und ruhendem Beobachter, r=100, θ=70°, rk=3, ri=3, ra=7):

Aberration mit vr=vθ=0, vφ=0.95 (ArcSin[vφ]=71.8°) auf θ=70°; rk=3, ri=4, ra=7; 1. Zeile: r=40, r=30, 2. Zeile: r=20, r=10:


◎ zusätzliche Rot/Blauverschiebung durch die Geschwindigkeit und Bewegungsrichtung des Beobachters:

Code: Alles auswählen
pic=Import["http://yukterez.net/mw/gradient2.png"]; (* Frequenzgradient *)
vr = -0.5; (* radiale Geschwindigkeitskomponente *)
vφ = +0.5; (* azimutale Geschwindigkeitskomponente *)
vθ = +0.5; (* polare Geschwindigkeitskomponente *)
ς = +1; (* Lapse *)
width = 120; (* Bildhöhe *)
height = 60; (* Bildbreite *)
v = +Sqrt[vr^2+vθ^2+vφ^2];
If[v>1, Print["Error: v > 1: v" -> v], Print["v" -> v]];
dΘ = +Limit[ArcCos[-vθ/Sqrt[vR^2+vθ^2+vφ^2]], vR->vr];
dθ = If[vr>0, -dΘ, +dΘ];
dΦ = +Limit[ArcTan[Abs[vφ]/vR], vR->vr];
dφ = If[vr!=0, -1, 1] If[vφ<0, +1, -1] If[NumericQ[dΦ], dΦ, π/2];
dψ = 0 If[vr<0, 0, π];
f[{φ_,θ_}] := {0, ς/(1-v Cos[θ-π/2]) Sqrt[1-v^2]}; (* winkelabhängige Frequenz *)
img = ImageTransformation[pic, f, {width, height}, (* Frequenzbereich *)
DataRange->{{-1, 1}, {0, 2}},
PlotRange->{{-π, π}, {-π/2, π/2}},
Padding->"Fixed"];
xyz[{x_,y_}] := {Sin[y] Cos[x],Sin[y] Sin[x],Cos[y]}; (* Rotationsmatrix *)
Xyz[{x_,y_,z_},dφ_] := {x Cos[dφ]-y Sin[dφ],x Sin[dφ]+y Cos[dφ],z};
xYz[{x_,y_,z_},dθ_] := {x Cos[dθ]+z Sin[dθ],y,z Cos[dθ]-x Sin[dθ]};
xyZ[{x_,y_,z_},dψ_] := {x,y Cos[dψ]-z Sin[dψ],y Sin[dψ]+z Cos[dψ]};
xy[{x_,y_,z_}] := {ArcTan[x,y],ArcCos[z]};
fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]}; (* Range *)
pct = ImageTransformation[img, fpt, DataRange->{{-1, 1}, {0, 1}},
PlotRange->{{-1, 1}, {-1, 1}}, Padding->"Periodic"];
rm[pct_,dφ_,dθ_,dψ_] := xy[xyZ[xYz[Xyz[xyz[pct], dφ], dθ], dψ]]; (* Rotation *)
RM[{x_,y_}] := rm[{x,y}, dφ, dθ, dψ];
red = ImageTransformation[pct, RM, (* Output im 360°x180° Plattkartenformat *)
DataRange->{{-π, π}, {-π, π}},
PlotRange->{{-π, π}, {0, π}},
Padding->"Periodic"]

◎ Solver für ISCO und Kreisbahngeschwindigkeit:

Code: Alles auswählen
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* || Orbital Velocity & ISCO Solver | yukterez.net | 16.07.2019 | Simon Tyran, Vienna || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
ClearAll["Global`*"];
a = 0.7; (* Spinparameter *)
℧ = 0.7; (* spezifische Ladung des schwarzen Lochs *)
si = rP; (* untere Grenze, prograder Photonenkreis *)
sr = 10; (* obere Grenze *)
st = 0.02; (* Interpolationsintervall für den Plot *)
Σ[я_] := я^2; (* Komponenten für die äquatoriale Ebene *)
Δ[я_] := я^2-2 я+a^2+℧^2;
Χ[я_] := (я^2+a^2)^2-a^2 Δ[я];
κ[я_] := a;
rA = 1+Sqrt[1-a^2-℧^2]; (* Horizont *)
rE = 1+Sqrt[1-℧^2]; (* Ergosphäre *)
R0 = If[Element[rA, Reals], rA, 0]; (* Mindestradius *)
rp = rf/.Solve[4 a^2 (rf-℧^2)-(rf^2-3 rf+2 ℧^2)^2 == 0 && rf >= R0, rf];
rP = Min[rp]; Rp = Max[rp]; (* prograder und retrograder Photonenkreis *)
isco = (* innermost stable circular orbit *)
Quiet[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 >= R0, RI]];
{"r horizon" -> [email protected], "r ergosphere" -> [email protected], "r isco" -> [email protected][isco],
"r photon pro" -> [email protected][rp], "r photon ret" -> [email protected][rp], "r disk" -> [email protected]}
j[v_] := Sqrt[1-v^2]; (* inverser Lorentzfaktor *)
Ы[я_] := Sqrt[Χ[я]/Σ[я]]; (* axialer Gyrationsradius *)
ωs[я_] := (a (2 я - ℧^2))/Χ[я]; (* Frame Dragging *)
ε[я_] := Sqrt[Δ[я] Σ[я]/Χ[я]]/j[v]+Lz[я] ωs[я]; (* Gesamtenergie *)
Lz[я_] := v Ы[я]/j[v]; (* Bahndrehimpuls*)
red[я_] := Quiet[Reduce[ (* Gleichungen *)
dt == (Lz[я] (-a (a^2+я^2)+Δ[я] κ[я])+ε[я] ((a^2+я^2)^2-Δ[я] κ[я]^2))/(Δ[я] Σ[я])
&&
0 == ((a^2+(-2+я) я+℧^2) (16 a dt dΦ я (я-℧^2)+8 dt^2 я (-я+℧^2)+dΦ^2 я (8 я (-a^2+
я^3)+a^2 (4 a^2+4 ℧^2-4 (a-℧) (a+℧)))))/(8 я^6)
&&
dΦ == (Lz[я] (-a^2+Δ[я])+ε[я] (a (a^2+я^2)-Δ[я] κ[я]))/(Δ[я] Σ[я])
&&
v > 0,
{v, dΦ, dt},
Reals]];
(* Lösung nach Radius *)
sol[x_, я_] := If[[email protected][red[я][[x, 2]]], red[я][[x, 2]], 0]
(* Interpolationstabelle *)
vs = Interpolation[Table[{я, sol[1, я]}, {я, si, sr, st}]];
φs = Interpolation[Table[{я, sol[2, я]}, {я, si, sr, st}]];
ts = Interpolation[Table[{я, sol[3, я]}, {я, si, sr, st}]];
(* Plotfunktion *)
plot[func_, label_] := Plot[func, {я, rP, sr},
GridLines -> {{Min[rp], Max[rp], rA, si, Min[isco], Max[isco], rE, sr}, {}},
Frame -> True, ImagePadding -> {{40, 1}, {12, 1}}, ImageSize -> 340, PlotLabel -> label,
PlotRange->{{0, sr}, Automatic}]
(* Plots *)
plot[Sqrt[Χ[я]/Δ[я]/Σ[я]], "Gravitational time dilation: y=dt/dт, x=r"]
plot[ts[я], "Total time dilation: y=dt/dτ, x=r"]
plot[(a (2 я-℧^2))/((a^2+я^2)^2-a^2 (a^2-2 я+я^2+℧^2)), "Frame dragging: y=dφ/dт, x=r"]
plot[φs[я]/ts[я], "Shapirodelayed angular velocity: y=dφ/dt, x=r"]
plot[φs[я], "Coordinate speed: y=dφ/dτ, x=r"]
plot[vs[я], "Local velocity: y=v=dl/dτ, x=r"]
r0 = Min[isco]; (* Radialkoordinate *)
"r isco" -> r0 "GM/c²" (* r isco *)
"dt/dτ " -> sol[3, r0] "sec/sec" (* totale ZD dt/dτ *)
"dt/dт " -> Sqrt[Χ[r0]/Δ[r0]/Σ[r0]] "sec/sec" (* gravitative ZD dt/dт *)
"dφ/dτ " -> sol[2, r0] "c³/G/M" (* Koordinatenschnelligkeit dφ/dτ *)
"dφ/dt " -> sol[2, r0]/sol[3, r0] "c³/G/M" (* shapiroverzögerte Winkelgeschw. dφ/dt *)
"v " -> sol[1, r0] "c" (* prograde Kreisbahngeschwindigkeit v lokal *)

⍟ Stereographische Projektion:

Code: Alles auswählen
Eq2St[{x_,y_}]:={0+ArcTan[-y,x],-2ArcTan[2,Sqrt[x^2+y^2]]+π/2};
St=ImageTransformation[pic,Eq2St,DataRange->{{-π,π},{-π/2,π/2}},PlotRange->{{-2π,2π},{-2π,2π}}]

Code Syntax: Mathematica. Archiv der alten Versionen, Testmodule und Extras:


verwandte Themen:
