Das ist die deutschsprachige Version. 
English versions will be available on en.yukterez.net and yukipedia. 
Last update: 30.7.2019Inhalt: Codes || Leinwand || Milchstraße || Landschaft || Aberration || Nackte Singularitäten || Akkretionsscheibe || Anhang
ungelinste (oben) und verzerrte Ansicht (unten) eines Kerr Newman SL mit Akkretionsscheibe (hier die geometrische Verzerrung ohne Rotverschiebung), FOV=90°×45°:
Codes (Syntax: Wolfram)
Aktuelle Version (all in one): Output der verzerrten Scheibengeometrie mit und ohne Kreisgeschwindigkeit derselben, Frequenzverschiebung, Hintergrund- und Horizontverzerrung. Als Layer für die umspannende Kugelschale werden Bilder im Plattkartenformat verwendet, andere Formate müssen zuerst in dieses Format transformiert werden. Für die Aberration kann die Lokalgeschwindigkeit des Beobachters eingegeben werden. Deren Betrag muss kleiner als 1 sein, und wird relativ zum ZAMO bemessen. Ein negatives vr steht für eine radiale Bewegung auf das schwarze Loch zu und ein positives von ihm weg, ein positives vφ für eine prograde und ein negatives für eine retrograde, und ein positives vθ im Bereich φ=-π/2..+π/2 für eine Bewegung in Richtung Südpol, und ein negatives in Richtung Nordpol (und umgekehrt in der gegenüberliegenden Kugelschalenhälfte). Falls in gewissen Konfigurationen (insbesondere mit dem Beobachter sehr nahe am Horizont) Artefakte auftreten kann R0 verringert und der Betrag von tmax erhöht werden. Boyer Lindquist Koordinaten (in der Nähe des Horizonts numerisch instabil, geeignet für Beobachter die sich in weiter Entfernung vom schwarzen Loch befinden):
Wenn der Beobachter nahe oder hinter dem Horizont ist empfehlen sich Raindrop Doran Koordinaten; ein relativ zum schwarzen Loch stationärer ZAMO hat in diesen Koordinaten eine lokale radiale Auswärtsgeschwindigkeit von vr=+√((2r-℧)/(r²+a²)):
Relativistische Geschwindigkeitsaddition:
Raytracing (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)
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:
Solver für ISCO und Kreisbahngeschwindigkeit:
Stereographische Projektion:
verwandte Themen:Kartographie, Kerr-Newman Metrik, Gravitationslinsen
Code: Alles auswählen
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 - 04.06.2020 | Version 12 | 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 = 240; (* Zielabmessungen in Pixeln *)
hoehe = 120; (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom = 1; (* 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["http://yukterez.net/mw/worldmap1.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 = 10; (* Radialkoordinate des Beobachters *)
R0 = 500; (* Radius des umspannenden Kugelschalenpanoramas *)
R1 = Max[1, Re[1.0001 rA]]; (* Sphäre *)
si = 3; (* 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 = 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[Χ/Δ/Σ];
μ = If[Abs[v0]==1, 0, -1]; (* Partikel: μ=-1, Photon: μ=0 *)
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][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}]];
"Accretion disk:"
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, ft, fτ, varb,
ft0s, ft1s, ft2s, ft3s, ft4s,
ft0v, ft1v, ft2v, ft3v, ft4v,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b,
dt0y, dr0y, dθ0y, dφ0y},
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]];
initcon = FindInstance[
dt0y ==
-((2 a ℧^2 Sin[θi]^2 dφ0y-4 a r0 Sin[θi]^2 dφ0y+a^2 Cos[θi]^2 \[Sqrt]((1/((a^2 Cos[θi]^2+
r0^2)^2))(4 a^2 (℧^2-2 r0)^2 Sin[θi]^4 dφ0y^2-(1/(a^2+℧^2-2 r0+r0^2))(-a^2-℧^2+2 r0-
r0^2+a^2 Sin[θi]^2) ((a^2+a^2 Cos[2 θi]+2 r0^2)^2 dr0y^2+(a^2+℧^2-2 r0+
r0^2) (-4 μ (a^2 Cos[θi]^2+r0^2)+(a^2+a^2 Cos[2 θi]+2 r0^2)^2 dθ0y^2+4 Sin[θi]^2 ((a^2+
r0^2)^2-a^2 (a^2+℧^2-2 r0+r0^2) Sin[θi]^2) dφ0y^2))))+r0^2 \[Sqrt]((1/((a^2 Cos[θi]^2+
r0^2)^2))(4 a^2 (℧^2-2 r0)^2 Sin[θi]^4 dφ0y^2-(1/(a^2+℧^2-2 r0+r0^2))(-a^2-℧^2+2 r0-
r0^2+a^2 Sin[θi]^2) ((a^2+a^2 Cos[2 θi]+2 r0^2)^2 dr0y^2+(a^2+℧^2-
2 r0+r0^2) (-4 μ (a^2 Cos[θi]^2+r0^2)+(a^2+a^2 Cos[2 θi]+2 r0^2)^2 dθ0y^2+
4 Sin[θi]^2 ((a^2+r0^2)^2-a^2 (a^2+℧^2-2 r0+r0^2) Sin[θi]^2) dφ0y^2)))))/(-2 (a^2+
℧^2-2 r0+r0^2)+2 a^2 Sin[θi]^2))
&&
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];
DGL = { (* Kerr Newman Bewegungsgleichungen *)
t''[τ]==-(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, 1]],
t[0]==0,
r''[τ]==-(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, 2]],
r[0]==r0,
θ''[τ]==-(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, 3]],
θ[0]==θi,
φ''[τ]==-(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, 4]],
φ[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[τ]==N[R1]||r[τ]<R1 &&
NumericQ[plunge5]==False,
(plunge5=τ) && (tt5=t[τ]) && (rt5=r[τ]) && (θt5=θ[τ]) && (φt5=φ[τ])],
WhenEvent[r[τ]==N[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,
Method-> mta,
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 = 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]];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 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]
Quit[]
Wenn der Beobachter nahe oder hinter dem Horizont ist empfehlen sich Raindrop Doran Koordinaten; ein relativ zum schwarzen Loch stationärer ZAMO hat in diesen Koordinaten eine lokale radiale Auswärtsgeschwindigkeit von vr=+√((2r-℧)/(r²+a²)):
Code: Alles auswählen
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 - 04.06.2020 | Version 12 | 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 = 240; (* Zielabmessungen in Pixeln *)
hoehe = 120; (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom = 1; (* 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["http://yukterez.net/mw/worldmap1.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 = 10; (* Radialkoordinate des Beobachters *)
R0 = 500; (* Radius des umspannenden Kugelschalenpanoramas *)
R1 = Max[1, Re[1.0001 rA]]; (* Sphäre *)
si = 3; (* 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 = vR;(* Radiale Geschwindigkeit relativ zum Raindrop, 0 für Fluss, +vR für stationär *)
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[
[email protected][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[
[email protected][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 |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
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;
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];
ς = 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][rp], "r photon ret" -> [email protected][rp], "r disk" -> [email protected],
"r observer" -> N[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}]];
"Accretion disk:"
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)]);
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,
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,
ft0f, ft1f, ft2f, ft3f, ft4f,
ft5h, ft5b,
dt0y, dr0y, dθ0y, dφ0y,
initcon},
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]];
initcon = FindInstance[
dr0y == -((Sqrt[2] (a^2+r0^2) Sqrt[(a^2+2 r0^2+a^2 Cos[2 θi])/(a^2+r0^2)] (-vr0i+
(Sqrt[2] Sqrt[2 r0-℧^2] (1-a dφ0y Sin[θi]^2))/(Sqrt[a^2+r0^2] Sqrt[(a^2+2 r0^2+
a^2 Cos[2 θi])/(a^2+r0^2)])))/(a^2+2 r0^2+a^2 Cos[2 θi]))
&&
dθ0y == vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2]
&&
dφ0y == -(((r0^2+a^2 Cos[θi]^2) Csc[θi]^2 Sqrt[(Sin[θi]^2 ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+
℧^2) Sin[θi]^2))/(r0^2+a^2 Cos[θi]^2)] (-vφ0i+(a (-2 r0+℧^2) Sin[θi]^2)/((r0^2+
a^2 Cos[θi]^2) Sqrt[(Sin[θi]^2 ((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 dr0y Sqrt[2 r0-℧^2] (a^2+2 r0^2+a^2 Cos[2 θi]) Sin[θi]^2)/(2 Sqrt[a^2+
r0^2] (r0^2+a^2 Cos[θi]^2) Sqrt[(Sin[θi]^2 ((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)),
{dr0y, dθ0y, dφ0y}];
DGL = { (* Kerr Newman Bewegungsgleichungen *)
t''[τ] == 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] == 1,
t[0] == 0,
r''[τ] == 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, 1]],
r[0] == r0,
θ''[τ] == -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, 2]],
θ[0] == θi,
φ''[τ] == 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, 3]],
φ[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[τ]==N[R1] &&
NumericQ[plunge5]==False,
(plunge5=τ) && (tt5=t[τ]) && (rt5=r[τ]) && (θt5=θ[τ]) && (φt5=φ[τ])],
WhenEvent[r[τ]==N[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,
Method-> mta,
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 = [email protected][{Ф, ϑ}];
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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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 *)
[email protected][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]
Quit[]
Relativistische Geschwindigkeitsaddition:
Code: Alles auswählen
V={vr, vθ, vφ}; (* Input: Geschwindigkeit relativ zum ZAMO *)
U={-Sqrt[(2 r0-℧^2)/(a^2+r0^2)], 0, 0}; (* Raumflussgeschwindigkeit *)
γ=1/Sqrt[1-Norm[U]^2]; (* Gamma *)
W=(-U+V+γ/(1+γ)(-U\[Cross](-U\[Cross]V)))/(1-U.V) (* Output: Doran Geschwindigkeitsvektor *)Raytracing (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 | *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
wp = MachinePrecision;
mt0 = Automatic;
mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
mta = mt0; (* mt0 for speed, mt1 for accuracy *)
Pause[1] (* Minkowski Koordinaten *)
kernels = 6; (* Parallelisierung *)
grain = 5; (* Subparallelisierung auf kernels*grain Streifen *)
rsp = "Nearest"; (* Resampling *)
breite = 240; (* 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, 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,
Method-> mta,
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}]}]
Quit[]
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 || *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
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]Q[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π}}]verwandte Themen:
Kartographie, Kerr-Newman Metrik, Gravitationslinsen
