




☑ Version A: rechtwinkelige Leinwand; Strahlen die ihren Ursprung außerhalb der Leinwand haben werden nicht gerendert, da diese als die einzige Lichtquelle in einer ansonsten dunklen Umgebung behandelt wird:

☑ Version B: sphärischer Hintergrund um auch Photonen die von jenseits der Leinwand kommen einzufangen. Als Layer für die umspannende Kugelschale können Bilder im Plattkartenformat verwendet werden; andere Formate müssen zuerst in dieses Format transformiert werden. Für die Aberration kann eine Orbitalgeschwindigkeit des Beobachters eingegeben werden. Deren Betrag muss kleiner als +1 sein, und wird relativ zum lokal drehimpulsfreien ZAMO bemessen. Ein positives vr steht für eine radiale Bewegung auf das schwarze Loch zu und ein negatives 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 Nordpol, und ein negatives in Richtung Südpol (und umgekehrt in der gegenüberliegenden Kugelschalenhälfte):

◈ Hintergrundpanorama: wenn der Beobachter über dem Nordpol des schwarzen Lochs (θ0=0±epsilon) platziert wird blickt er auf den Südpol der umspannenden Kugelschale. Wenn der Nordpol des SL vor der äquatorialen Ebene des Hintergrundbilds betrachtet werden soll muss dieses zuvor kartographisch transformiert werden. Rotation um θ:

⍟ Performance Boost: wenn verschiedene Bilder mit den selben Einstellungen geraytraced werden sollen (beispielsweise bei einem Orbit auf konstantem θ) kann der Vorgang mithilfe eines Lookup-Tables extrem beschleunigt werden. Das erste Bild dauert dann normal lang, während alle weiteren sehr schnell gehen. Außerdem kann jeder Kernel einen eigenen Teil des Bildes rendern, und die verschiedenen Teile danach mit ImageAssemble zusammenfügen:

◬ Der Code ist zwar noch nicht für GPU-Beschleunigung optimiert, kann aber auch über die CPU parallelisiert werden indem die PlotRange beispielsweise geviertelt und jedes Viertel parallel von einem eigenen Kernel gerendert wird, siehe Screenshot. Bei 4 vorhandenen Kernels kann die Rechendauer damit schon einmal geviertelt werden, sofern man bereit ist 100% CPU-Auslastung zu akzeptieren. Wenn viele Kernels vorhanden sind kann der Vorgang auch mit ParallelSubmit[...] und Parallelize[...] beschleunigt werden. Die Genauigkeit wird über die Parameter mta und wp sowie über den EventLocator geregelt.

Code Syntax: Mathematica. Prototyp der nächsten Version, Testmodule, Archiv und Extras:
Code: Alles auswählen
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 | 3A | Kerr Newman Metrik | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
ClearAll["Global`*"]
rA = 1+Sqrt[1-a^2-℧^2];
wp = MachinePrecision;
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
r0 = Sqrt[X0^2-a^2]; (* Startradius *)
θ0 = π/2; (* Breitengrad *)
φ0 = 0; (* Längengrad *)
a = 9/10; (* Spinparameter *)
℧ = 2/5; (* spezifische Ladung des schwarzen Lochs *)
v0 = 1; (* Anfangsgeschwindigkeit *)
μ = 0; (* Photon *)
X0 =-20; (* Ebene des verzerrten Bildes *)
tmax = 4/3(X0-r0); (* zeitlicher Integrationsbereich *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 2) Geschwindigkeitskomponenten der auf der Sichtebene eintreffenden Strahlen ||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
fPlaneA[{x_,y_,z_}] := Solve[1==(x/y vφ0)^2+(z/y vφ0)^2+vφ0^2 && vr0==x/y vφ0 && vθ0==z/y vφ0 && vr0>0, {vr0,vθ0,vφ0}];
fPlaneB[{x_,y_,z_}] := Solve[1==(x/z vθ0)^2+(y/z vθ0)^2+vθ0^2 && vr0==x/z vθ0 && vφ0==y/z vθ0 && vr0>0, {vr0,vθ0,vφ0}];
fPlaneC[{x_,y_,z_}] := Solve[1==vr0 && 0==vφ0 && 0==vθ0, {vr0,vθ0,vφ0}];
fPlaneD[{x_,y_,z_}] := If[y==0, If[z==0, fPlaneC[{x,y,z}], fPlaneB[{x,y,z}]], fPlaneA[{x,y,z}]];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 3) Metrische Koeffizienten ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
Σ = r0^2+a^2 Cos[θ0]^2;
Δ = r0^2-2 r0+a^2+℧^2;
Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
gtt = (2r0-℧^2)/Σ-1;
grr = Σ/Δ;
gθθ = Σ;
gφφ = Χ/Σ Sin[θ0]^2;
gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;
fPlane0[{y_,z_}] := fPlaneD[{2X0,y,z}];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 4) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
raytrace[{ϒ_,Ζ_}] :=
Quiet[Module[{tMax, vr0i, vθ0i, vφ0i, vr0a, vθ0a, vφ0a, vt0a, DGL, sol, ε, Lz, pθ0, pr0, Q, k, t10, r10, Θ10, Φ10, т0, plunge, R1, X1, X, Y, Z, rt, θt, φt, т, τ, t, r, θ, φ, stepsize, laststep, mta},
vr0a = (vr0/.fPlane0[{ϒ,Ζ}][[1,1]]) Sqrt[grr];
vθ0a = (vθ0/.fPlane0[{ϒ,Ζ}][[1,2]]) Sqrt[gθθ]/r0;
vφ0a = (vφ0/.fPlane0[{ϒ,Ζ}][[1,3]]) Sqrt[gφφ]/r0 Sin[θ0];
vt0a = Sqrt[vr0a^2+vφ0a^2+vθ0a^2];
vr0i = vr0a/vt0a;
vφ0i = vφ0a/vt0a;
vθ0i = vθ0a/vt0a;
mta = {"EventLocator","Event"->(r[τ]-101/100 rA)};
DGL = {
t''[τ]==(4 (((a^2+a^2 Cos[2 θ[τ]]+2 (℧^2-r[τ]) r[τ]) (a^2+r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+a^2 (-℧^2+2 r[τ]) Sin[2 θ[τ]] t'[τ] θ'[τ]-(1/(a^2+℧^2-2 r[τ]+r[τ]^2))a (a^4+4 ℧^2 r[τ]^3-6 r[τ]^4-3 a^2 r[τ] (-℧^2+r[τ])+a^2 Cos[2 θ[τ]] (a^2+(℧^2-r[τ]) r[τ])) Sin[θ[τ]]^2 r'[τ] φ'[τ]-2 a^3 Cos[θ[τ]] (-℧^2+2 r[τ]) Sin[θ[τ]]^3 θ'[τ] φ'[τ]))/(a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^2,
t'[0]==-((a (2 r0-℧^2) Sin[θ0]^2 (vφ0i (-2 r0+r0^2+a^2 Cos[θ0]^2) Csc[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a vφ0i (2 r0-℧^2) Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2))))/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θ0]^2)))+\[Sqrt](((vr0i^2 (a^2-2 r0+r0^2+℧^2)+vθ0i^2 (a^2-2 r0+r0^2+℧^2)) (r0^2+a^2 Cos[θ0]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θ0]^2)+(a^2 (-2 r0+℧^2)^2 Sin[θ0]^4 (vφ0i (-2 r0+r0^2+a^2 Cos[θ0]^2) Csc[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a vφ0i (2 r0-℧^2) Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)))^2)/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)^2)-((2 r0-r0^2-℧^2-a^2 Cos[θ0]^2) Sin[θ0]^2 ((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2) (vφ0i (-2 r0+r0^2+a^2 Cos[θ0]^2) Csc[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)]+(a vφ0i (2 r0-℧^2) Sin[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θ0]^2)))^2)/((a^2-2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θ0]^2)^2))/((a^2-2 r0+r0^2+℧^2) (-2 r0+r0^2+℧^2+a^2 Cos[θ0]^2)^2)),
t[0]==0,
r''[τ]==(1/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((8 (a^2 Cos[θ[τ]]^2+(a^2+℧^2-a^2 Cos[θ[τ]]^2) r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) (t'[τ])^2+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] r'[τ] θ'[τ]+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (a^2+℧^2-2 r[τ]+r[τ]^2) (θ'[τ])^2-16 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+(a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 (a^2 (3 a^2+4 ℧^2+4 (a^2-℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]+16 a^2 Cos[θ[τ]]^2 r[τ]^3+8 r[τ]^5-8 a^2 r[τ]^2 Sin[θ[τ]]^2+2 a^4 Sin[2 θ[τ]]^2) (φ'[τ])^2),
r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θ0]^2)/(a^2+(-2+r0) r0+℧^2)],
r[0]==r0,
θ''[τ]==(1/(16 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))-8 a^2 (℧^2-2 r[τ]) Sin[2 θ[τ]] (t'[τ])^2-32 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (θ'[τ])^2+16 a (℧^2-2 r[τ]) (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ] φ'[τ]+(a^4 (3 a^2-5 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+(a^2+℧^2) Cos[4 θ[τ]])+a^2 (11 a^2-8 ℧^2+4 (3 a^2+2 ℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]^2+8 a^2 (2+Cos[2 θ[τ]]) r[τ]^4+8 r[τ]^6+8 a^4 (3+Cos[2 θ[τ]]) r[τ] Sin[θ[τ]]^2+32 a^2 r[τ]^3 Sin[θ[τ]]^2) Sin[2 θ[τ]] (φ'[τ])^2),
θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θ0]^2],
θ[0]==θ0,
φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))(-((8 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 a Cot[θ[τ]] (℧^2-2 r[τ]) t'[τ] θ'[τ]+((a^2 (3 a^2+8 ℧^2+4 a^2 Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) r'[τ] φ'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+Cot[θ[τ]] (a^2 (3 a^2-4 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+16 a^2 Cos[θ[τ]]^2 r[τ]^2+8 r[τ]^4+16 a^2 r[τ] Sin[θ[τ]]^2) θ'[τ] φ'[τ]),
φ'[0]==(vφ0i ((-2+r0) r0+a^2 Cos[θ0]^2) Csc[θ0] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)]+(a vφ0i (2 r0-℧^2) Sin[θ0])/((r0^2+a^2 Cos[θ0]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θ0]^2)),
φ[0]==φ0
};
sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
Method-> mta,
InterpolationOrder-> All,
StepMonitor :> (laststep=plunge; plunge=τ;
stepsize=plunge-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<2*^-2, 0, 1])}];
tMax = Max[tmax, plunge+1/10];
X[τ_] := Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
Y[τ_] := Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
Z[τ_] := Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];
т[coord_,dist_] := Quiet[ξ/.FindRoot[coord[ξ]-dist, {ξ,tMax 9/10,tmax,-1}]];
т0 = т[X,X0];
X1 = X[т0];R1=Evaluate[r[т0]/.sol][[1]];
Quiet[If[Round[X1] == Round[X0],{+Y[т0],-Z[т0]},If[R1>3,{0,0},{0,30}]]]]]
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 5) Testbild laden und transformieren ||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
Eq=Import["http://666kb.com/i/ds8dprdq40nslgzsk.png"] (* Testbild *)
ImageTransformation[Eq,raytrace,DataRange->{{-30,30},{-30,30}},PlotRange->{{-30,30},{-30,30}},Padding->"Reflected"]
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* Code stabil, aber noch nicht auf Geschwindigkeit optimiert ||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

☑ Version B: sphärischer Hintergrund um auch Photonen die von jenseits der Leinwand kommen einzufangen. Als Layer für die umspannende Kugelschale können Bilder im Plattkartenformat verwendet werden; andere Formate müssen zuerst in dieses Format transformiert werden. Für die Aberration kann eine Orbitalgeschwindigkeit des Beobachters eingegeben werden. Deren Betrag muss kleiner als +1 sein, und wird relativ zum lokal drehimpulsfreien ZAMO bemessen. Ein positives vr steht für eine radiale Bewegung auf das schwarze Loch zu und ein negatives 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 Nordpol, und ein negatives in Richtung Südpol (und umgekehrt in der gegenüberliegenden Kugelschalenhälfte):
Code: Alles auswählen
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 25.04.2018 | 5B | Kerr Newman Metrik | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
ClearAll["Global`*"]
wp = MachinePrecision;
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
r0 = 50; (* Radialkoordinate des Beobachters *)
θ0 = π/2; (* Breitengrad *)
φ0 = 0; (* Längengrad *)
R0 = 1000; (* Radius des umspannenden Kugelschalenpanoramas *)
tmax =-4/3 R0; (* zeitlicher Integrationsbereich *)
a = 1; (* Spinparameter *)
℧ = 0; (* spezifische Ladung des schwarzen Lochs *)
v0 = 1; (* Geschwindigkeit des Photons *)
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 *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 2) Metrische Koeffizienten und Formeln ||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
Σ = r0^2+a^2 Cos[θ0]^2;
Δ = r0^2-2 r0+a^2+℧^2;
Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
μ = 0;
gtt = (2r0-℧^2)/Σ-1;
grr = Σ/Δ;
gθθ = Σ;
gφφ = Χ/Σ Sin[θ0]^2;
gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;
θi =-θ0+π;
rA = 1+Sqrt[1-a^2-℧^2];
й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];
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 3) Rotationsmatrix für die auf der Sichtebene eintreffenden Strahlen ||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 4) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
raytrace[{Ф_,ϑ_}] :=
Quiet[Module[{V, W, vw, tMax, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a, DGL, sol, ε, Lz, pθi, pr0, Q, k, t10, r10, Θ10, Φ10, т0, R1, t, r, θ, φ, τ, plunge, X, Y, Z, rt, θt, φt, т, ξ, stepsize, laststep, mta},
vw=xyZ[Xyz[{0,1,0},ϑ],Ф+π/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]];
(* Integrationsende *)
mta={"EventLocator","Event"->If[(r[τ]==1.01rA || r[τ] == R0+1.0) == True, 0, 1]};
DGL = { (* Kerr Newman Bewegungsgleichungen *)
t''[τ]==(4 (((a^2+a^2 Cos[2 θ[τ]]+2 (℧^2-r[τ]) r[τ]) (a^2+r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+a^2 (-℧^2+2 r[τ]) Sin[2 θ[τ]] t'[τ] θ'[τ]-(1/(a^2+℧^2-2 r[τ]+r[τ]^2))a (a^4+4 ℧^2 r[τ]^3-6 r[τ]^4-3 a^2 r[τ] (-℧^2+r[τ])+a^2 Cos[2 θ[τ]] (a^2+(℧^2-r[τ]) r[τ])) Sin[θ[τ]]^2 r'[τ] φ'[τ]-2 a^3 Cos[θ[τ]] (-℧^2+2 r[τ]) Sin[θ[τ]]^3 θ'[τ] φ'[τ]))/(a^2+a^2 Cos[2 θ[τ]]+2 r[τ]^2)^2,
t'[0]==-((a (2 r0-℧^2) Sin[θ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/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((8 (a^2 Cos[θ[τ]]^2+(a^2+℧^2-a^2 Cos[θ[τ]]^2) r[τ]-r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) (t'[τ])^2+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] r'[τ] θ'[τ]+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 (a^2+℧^2-2 r[τ]+r[τ]^2) (θ'[τ])^2-16 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) (a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+(a^2+℧^2-2 r[τ]+r[τ]^2) Sin[θ[τ]]^2 (a^2 (3 a^2+4 ℧^2+4 (a^2-℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]+16 a^2 Cos[θ[τ]]^2 r[τ]^3+8 r[τ]^5-8 a^2 r[τ]^2 Sin[θ[τ]]^2+2 a^4 Sin[2 θ[τ]]^2) (φ'[τ])^2),
r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θi]^2)/(a^2+(-2+r0) r0+℧^2)],
r[0]==r0,
θ''[τ]==(1/(16 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3))(-((16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (r'[τ])^2)/(a^2+℧^2-2 r[τ]+r[τ]^2))-8 a^2 (℧^2-2 r[τ]) Sin[2 θ[τ]] (t'[τ])^2-32 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+16 a^2 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]] (θ'[τ])^2+16 a (℧^2-2 r[τ]) (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ] φ'[τ]+(a^4 (3 a^2-5 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+(a^2+℧^2) Cos[4 θ[τ]])+a^2 (11 a^2-8 ℧^2+4 (3 a^2+2 ℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]^2+8 a^2 (2+Cos[2 θ[τ]]) r[τ]^4+8 r[τ]^6+8 a^4 (3+Cos[2 θ[τ]]) r[τ] Sin[θ[τ]]^2+32 a^2 r[τ]^3 Sin[θ[τ]]^2) Sin[2 θ[τ]] (φ'[τ])^2),
θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2],
θ[0]==θi,
φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))(-((8 a (a^2 Cos[θ[τ]]^2+℧^2 r[τ]-r[τ]^2) r'[τ] t'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2))+8 a Cot[θ[τ]] (℧^2-2 r[τ]) t'[τ] θ'[τ]+((a^2 (3 a^2+8 ℧^2+4 a^2 Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]]) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) r'[τ] φ'[τ])/(a^2+℧^2-2 r[τ]+r[τ]^2)+Cot[θ[τ]] (a^2 (3 a^2-4 ℧^2+4 (a^2+℧^2) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+16 a^2 Cos[θ[τ]]^2 r[τ]^2+8 r[τ]^4+16 a^2 r[τ] Sin[θ[τ]]^2) θ'[τ] φ'[τ]),
φ'[0]==(vφ0i ((-2+r0) r0+a^2 Cos[θ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
};
(* Integrator *)
sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
Method-> mta,
InterpolationOrder-> All,
StepMonitor :> (laststep=plunge; plunge=τ;
stepsize=plunge-laststep;), Method->{"EventLocator",
"Event" :> (If[stepsize<2*^-2, 0, 1])}];
(* Integrationszeit *)
tMax = Max[tmax, plunge+1/10];
(* Raumkoordinaten *)
rt[τ_] := Evaluate[r[τ]/.sol][[1]];
θt[τ_] := Evaluate[θ[τ]/.sol][[1]]+π/2;
φt[τ_] :=-Evaluate[φ[τ]/.sol][[1]]-π;
(* Affiner Parameter bei Emission *)
т[coord_,dist_] := ξ/.FindRoot[coord[ξ]-dist, {ξ,tMax 9/10,tmax,-1}];
т0 = т[rt,R0];
R1 = rt[т0]; (* Check ob die Photonen von der Hemisphäre kommen *)
(* Berechung der Ursprungskoordinaten der Photonen *)
If[т0>0, {π,π/2}, If[Round[R1] == Round[R0],{φt[т0],θt[т0]},If[R1>3,{-π,-π/2},{π,π/2}]]]]]
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 5) Testbild laden und transformieren ||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
pic=Import["https://cdn.eso.org/images/thumb700x/eso0932a.jpg"]
width=ImageDimensions[pic][[1]];
fpt[{x_,y_}] := {If[y<0,x+1,x],If[y<0,-y,y]};
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
ImageTransformation[pct,raytrace,DataRange->{{-π,π-2π/width},{-π/2,3π/2}},PlotRange->{{-π/5,π/5},{-π/10,π/10}},Padding->"Periodic"]
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* Code stabil, aber noch nicht auf Geschwindigkeit optimiert ||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

◈ Hintergrundpanorama: wenn der Beobachter über dem Nordpol des schwarzen Lochs (θ0=0±epsilon) platziert wird blickt er auf den Südpol der umspannenden Kugelschale. Wenn der Nordpol des SL vor der äquatorialen Ebene des Hintergrundbilds betrachtet werden soll muss dieses zuvor kartographisch transformiert werden. Rotation um θ:
Code: Alles auswählen
Ep=Import["http://666kb.com/i/dsfr4u76175v8q277.png"]
width=ImageDimensions[Ep][[1]];
RM[{x_,y_}]:={ArcTan[Cos[x] Cos[ϑ] Sin[y]+Cos[y] Sin[ϑ], Sin[x] Sin[y]], ArcCos[Cos[y] Cos[ϑ]-Cos[x] Sin[y] Sin[ϑ]]};
ϑ=π/2-θ0; pic=ImageTransformation[Ep, RM, DataRange->{{-π,π-2π/width}, {0,π}}, PlotRange->{{-π,π}, {0,π}}]

⍟ Performance Boost: wenn verschiedene Bilder mit den selben Einstellungen geraytraced werden sollen (beispielsweise bei einem Orbit auf konstantem θ) kann der Vorgang mithilfe eines Lookup-Tables extrem beschleunigt werden. Das erste Bild dauert dann normal lang, während alle weiteren sehr schnell gehen. Außerdem kann jeder Kernel einen eigenen Teil des Bildes rendern, und die verschiedenen Teile danach mit ImageAssemble zusammenfügen:
Code: Alles auswählen
1) raytrace[{Ф_,ϑ_}] := Module[...] -> mem : raytrace[{Ф_,ϑ_}] := mem = Module[...]
2) ParallelDo[ImageTransformation[pct,raytrace,
DataRange->{{-π ,π },{-π/2 ,3π/2 }},PlotRange->plotrange,Padding->"Periodic"],
{plotrange,{{{-π,0},{-π/2,0}},{{-π,0},{0,π/2}},{{0,π},{-π/2,0}},{{0,π},{0,π/2}}}}]

◬ Der Code ist zwar noch nicht für GPU-Beschleunigung optimiert, kann aber auch über die CPU parallelisiert werden indem die PlotRange beispielsweise geviertelt und jedes Viertel parallel von einem eigenen Kernel gerendert wird, siehe Screenshot. Bei 4 vorhandenen Kernels kann die Rechendauer damit schon einmal geviertelt werden, sofern man bereit ist 100% CPU-Auslastung zu akzeptieren. Wenn viele Kernels vorhanden sind kann der Vorgang auch mit ParallelSubmit[...] und Parallelize[...] beschleunigt werden. Die Genauigkeit wird über die Parameter mta und wp sowie über den EventLocator geregelt.

Code Syntax: Mathematica. Prototyp der nächsten Version, Testmodule, Archiv und Extras: