Relativistischer Raytracer

[Yukterez]

|Kerr-Newman || Kerr || Schwarzschild || Geodätengleichung || Gummituchmodell || Gravitationslinsen || Projektionseffekt || Kartographie

   

Das ist die deutschsprachige Version.   English versions will be available on en.yukterez.net and yukipedia.

|Inhalt: Codes || Leinwand || Milchstraße || Landschaft || Aberrationseffekt || Nackte Singularitäten || Probe || Anhang

☑ 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:

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 G: 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 | 29.04.2018 - 15.06.2018 | Version 6G | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
ClearAll["Global`*"]
                                       
kernels = 6;                                                          (* Parallelisierung *)
grain   = 18;                           (* Subparallelisierung auf kernels*grain Streifen *)
breite  = 216;                                               (* Zielabmessungen in Pixeln *)
hoehe   = 108;          (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom    = 1;                                  (* doppelter Zoom ergibt halben Sichtwinkel *)
 
LaunchKernels[kernels]
wp = MachinePrecision;                                                     (* Genauigkeit *)

pic = Import["http://yukterez.net/mw/2/flip90.png"];         (* Hintergrundpanorama laden *)
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
r0   = 5;                                             (* Radialkoordinate des Beobachters *)
θ0   = π/2;                                                                (* Breitengrad *)
φ0   = 0;                                                                   (* Längengrad *)
 
R0   = 500;                              (* Radius des umspannenden Kugelschalenpanoramas *)
tmax =-5/3 R0;                                          (* zeitlicher Integrationsbereich *)
 
a    = 0.7;                                                              (* Spinparameter *)
℧    = 0.7;                                     (* 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 ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
raytracer[{Ф_, ϑ_}] :=
 
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, mtl, 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 *)
                                                   
mtl=If[a^2+℧^2>1,
{"EventLocator", "Event"->r[τ]-R0-1.0},
{"EventLocator", "Event"->If[(r[τ]==1.01rA || r[τ] == R0+1.0) == True, 0, 1]}];
mta=mtl;
 
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, {0, -π/2}, If[R1<10r0, {0, -π/2}, If[R1>10R0, {0, -π/2}, {φt[т0], θt[т0]}]]]]]
 
mem : raytrace[{Ф_, ϑ_}] := mem = raytracer[{Ф, ϑ}]
 
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 5) Testbild laden und transformieren ||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
 
fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]}
 
pcr = ParallelTable[
ImageTransformation[pic, fpt, DataRange->{{-1, 1}, {0, 1}}, PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Padding->"Periodic"],
{x, 1, 2 kernels}];
pct = ImageAssemble[{Table[{pcr[[x]]}, {x, 2 kernels, 1, -1}]}];
 
width = ImageDimensions[pic][[1]]; height = ImageDimensions[pic][[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];

FOV -> {360.0 hzoom "degree", 180.0 vzoom "degree"}
 
img = ParallelTable[
ImageTransformation[pct, raytrace, {breite, Ceiling[hoehe/kernels/grain]},
DataRange->{{-π, π-2π/width}, {-π/2, 3π/2}},
PlotRange->{{-π, π-2π/width} hzoom, {-π/2+x, -π/2+x+π/kernels/grain} vzoom},
Padding->"Periodic"],
{x, 0, π-π/kernels/grain, π/kernels/grain}]
 
source = ImageResize[pic, {2 hoehe, hoehe}]
image  = ImageAssemble[{Table[{img[[x]]}, {x, kernels grain, 1, -1}]}]













◈ 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->{{-π,π-2π/width}, {0,π}}]

⍟ 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π}}]

⧆ 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, Archiv der alten Version, Testmodule und Extras:

Akrretionscheibe (noch ungetestet):

Code: Alles auswählen

rt[τ_] := Evaluate[r[τ]/.sol][[1]];
θt[τ_] := Mod[Evaluate[θ[τ]/.sol][[1]]-Pi/2, Pi];
φt[τ_] := Evaluate[φ[τ]/.sol][[1]];
                                                        (* Affiner Parameter bei Emission *)
т[coord_,dist_] := Quiet[ξ/.FindRoot[coord[ξ]-dist, {ξ,-r0,tMax,-r0/2}]];
т0 = т[θt,0];
θ1 = θt[т0];                              (* Check ob die Photonen von der Scheibe kommen *)
                                       (* Berechung der Ursprungskoordinaten der Photonen *)
Quiet[If[rt[т0]>4+Pi, {Pi,4+Pi}, If[rt[т0]<4, {Pi,4+Pi}, {φt[т0], rt[т0]}]]]]]

Aufteilung in 4 Quadranten und äquatoriale Umrundung in 1° Schritten:

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* Quadrant 1 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Do[
pic=Import["F00" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["RO" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{0,Pi*860/3000},{0,512/1500*Pi/2}},Padding->"Periodic"]
]   
],
{tk, 1, 9, 1}]

Do[
pic=Import["F0" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["RO" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{0,Pi*860/3000},{0,512/1500*Pi/2}},Padding->"Periodic"]
]   
],
{tk, 10, 99, 1}]

Do[
pic=Import["F" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["RO" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{0,Pi*860/3000},{0,512/1500*Pi/2}},Padding->"Periodic"]
]   
],
{tk, 100, 360, 1}]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* Quadrant 2 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Do[
pic=Import["F00" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["LO" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{-Pi*860/3000,0},{0,512/1500*Pi/2}},Padding->"Periodic"]
]   
],
{tk, 1, 9, 1}]

Do[
pic=Import["F0" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["LO" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{-Pi*860/3000,0},{0,512/1500*Pi/2}},Padding->"Periodic"]
]   
],
{tk, 10, 99, 1}]

Do[
pic=Import["F" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["LO" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{-Pi*860/3000,0},{0,512/1500*Pi/2}},Padding->"Periodic"]
]   
],
{tk, 100, 360, 1}]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* Quadrant 3 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Do[
pic=Import["F00" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["RU" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{0,Pi*860/3000},{-512/1500*Pi/2,0}},Padding->"Periodic"]
]   
],
{tk, 1, 9, 1}]

Do[
pic=Import["F0" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["RU" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{0,Pi*860/3000},{-512/1500*Pi/2,0}},Padding->"Periodic"]
]   
],
{tk, 10, 99, 1}]

Do[
pic=Import["F" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["RU" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{0,Pi*860/3000},{-512/1500*Pi/2,0}},Padding->"Periodic"]
]   
],
{tk, 100, 360, 1}]

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* Quadrant 4 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

Do[
pic=Import["F00" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["LU" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{-Pi*860/3000,0},{-512/1500*Pi/2,0}},Padding->"Periodic"]
]   
],
{tk, 1, 9, 1}]

Do[
pic=Import["F0" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["LU" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{-Pi*860/3000,0},{-512/1500*Pi/2,0}},Padding->"Periodic"]
]   
],
{tk, 10, 99, 1}]

Do[
pic=Import["F" <> ToString[tk] <> ".png"]; fpt[{x_,y_}]:={If[y<0,x+1,x],If[y<0,-y,y]};
width=ImageDimensions[pic][[1]];
pct=ImageTransformation[pic,fpt,DataRange->{{-1,1},{0,1}},PlotRange->{{-1,1},{-1,1}},Padding->"Periodic"];
Export["LU" <> ToString[tk] <> ".png",
Rasterize[
ImageTransformation[pct,raytrace,DataRange->{{-Pi,Pi-2π/width},{-Pi/2,3Pi/2}},PlotRange->{{-Pi*860/3000,0},{-512/1500*Pi/2,0}},Padding->"Periodic"]
]   
],
{tk, 100, 360, 1}]













Kontrollmodul, Minkowski

Code: Alles auswählen

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 07.04.2018 | Minkowski Kontrollmodul | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

ClearAll["Global`*"]
wp = MachinePrecision;

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

r0   = 1;                                                                  (* Startradius *)
θ0   = π/2;                                                                (* Breitengrad *)
φ0   = 0;                                                                   (* Längengrad *)
v0   = 1;                                           (* Anfangsgeschwindigkeit des Photons *)
vφ   = 0;                                   (* Transversalgeschwindigkeit des Beobachters *)
μ    = 0;                                                                       (* Photon *)
R0   = 1000;                                               (* Ebene des verzerrten Bildes *)
tmax =-4/3 R0;                                          (* zeitlicher Integrationsbereich *)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 2) 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[ψ]};

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* 3) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

raytrace[{yy_,zz_}] :=

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

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

mta=Automatic;

vr0n = vw[[3]];
vφ0n = vw[[2]];
vθ0n = vw[[1]];

vφ0i = (vφ0n-vφ)/(1-vφ0n vφ);
vr0i = vr0n Sqrt[1-vφ^2]/(1-vφ0n vφ);
vθ0i = vθ0n Sqrt[1-vφ^2]/(1-vφ0n vφ);

DGL= {

t''[τ]==0,
t'[0]==Sqrt[φ0^2 r0^2 Sin[θ0]^2+vw[[3]]^2+(vw[[1]]/r0)^2 0^2],
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]==θ0,
 
φ''[τ]==-2 φ'[τ] (r'[τ]+r[τ] θ'[τ] Cot[θ[τ]])/r[τ],
φ'[0]==vφ0n/Csc[θ0]/r0,
φ[0]==φ0};

sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
Method-> mta];

tMax = tmax;

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

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

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

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

pic=Import["http://666kb.com/i/dsfr4u76175v8q277.png"]
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->{{-π,π-2π/width},{-π/2,π/2}},Padding->"Periodic"]













Sicherungskopie des alten 5er-Codes:

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 ||||||||||||||||||||||||||| *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)













Rohmaterial: Kartographie, Kerr-Newman Metrik, Gravitationslinsen, Sichtebene.nb

[Yukterez]

◩ Code Version A, rechtwinkelige Leinwand


Konstruktion der Sichtebene:

 
Die türkis gestrichelte Kurve zeigt was man dort wo man im flachen Raum die linke untere Ecke des Bildes sehen würde sieht (dadurch erscheint das gravitationsgelinste Bild größer als das unverzerrte). Vorläufige Testbilder mit M=Q=a=0 (Vakuum), a=Q=0 (Schwarzschild) und a=M (extremal Kerr), erstellt mit auf 500 begrenzten Integrationsschritten pro Pixel; Perspektive: äquatorial.

Das schwarze Loch befindet sich relativ zum Betrachter in einer kartesischen Entfernung von 20GM/c². Der Mittelpunkt des betrachteten Bilds das sich auf einer rigiden stationären Leinwand die rechtwinkelig zur Sichtachse des Beobachters steht befindet ist ebensoweit vom schwarzen Loch, und damit doppelt so weit vom Betrachter entfernt. Die projizierte Sichtebene beträgt 60x60(GM/c²)²; alle Lichtstrahlen die von wo anders als der Leinwand kämen werden grau dargestellt.

Die projizierte Sichtebene in der Entfernung des schwarzen Lochs ist wegen der halben Entfernung damit 30x30(GM/c²)². Der Durchmesser des Schattens des schwarzen Lochs beträgt ca. 1/3 und damit 10GM/c², was einem Radius von ca. 5GM/c² entspricht (zum Vergleich siehe die semianalytische Konturlösung auf gravitylense.yukterez.net):



In der unteren (also der mittleren) Reihe befindet der Beobachter sich 100 Mal weiter vom schwarzen Loch entfernt, in einem kartesischen Abstand von 2000GM/c². Die Leinwand befindet sich ebensoweit in die andere Richtung vom SL entfernt; durch den nun spitzeren Sichtwinkel und den kleineren Impact Parameter der Photonen ergibt sich jetzt eine sehr viel stärkere Verzerrung.

Ohne schwarzes Loch im Vordergrund würde man innerhalb des Sichtwinkels im Zoom das gleiche Bild wie oben links, bzw. unten links in der Mitte (die 4 bekannten Kästchen) sehen; das heißt die höhere Entfernung wird durch höheren Zoom kompensiert, so dass das schwarze Loch den gleichen Winkeldurchmesser hat wie im oberen Beispiel (die Kamera ist immer noch auf das nun 100 mal weiter entfernte und rot umrahmte Quadrat von 30GM/c² Kantenlänge fokussiert):

Der graue Bereich symbolisiert wieder den Bereich wo keine Photonen von der Leinwand, aber theoretisch welche von seitwärts derselben beim Betrachter ankommen können (siehe nächstes Kapitel). Im komplett schwarzen Bereich würden weder Photonen von der Leinwand noch von sonst irgendwo im Raum das Auge des Betrachters treffen.

Die gerade Leinwand wird auf den y,z-Achsen bis in die Unendlichkeit fortgesetzt und das Bild wie unten links auf dieselbe gekachelt, da man ansonsten nicht viel sehen würde weil das Bild in der Entfernung vom Schatten des schwarzen Lochs verdeckt werden würde. Unterste Reihe: nackte Singularitäten mit Q=2M, a²+Q²=2M² und a=2M, Abstände wie in der mittleren Reihe:

[Yukterez]

◐ Code Version B, Panorama


Konstruktion der Sichtebene:



Also im Prinzip wie oben, aber mit Mapping des gesamten Umfelds um auch Licht das von außerhalb der Leinwand ins Sichtfeld der Kamera gebogen wird einzufangen. Fokus auf einen kleinen Bereich simuliert eine Lochkamera, während eine Plot-Range von φ=-π..π und θ=-π/2..π/2 ein 360° Kugelpanorama ergibt. Das gesamte Kugelsphärenpanorama wird dann ins Plattkartenformat exportiert und kann danach in jedes beliebige andere Format transformiert oder anderweitig projiziert werden. In den folgenden Beispielen wird der Beobachter auf r=50GM/c² platziert und das Himmelszelt auf R=1000GM/c² aufgespannt.



Sphärischer Hintergund: Milchstraßenpanorama (360°-Aufnahme im Plattkartenformat, Fotocredit ESO/José Francisco). Zoom auf den Bereich vor dem das schwarze Loch platziert wird; ungelinste Ansicht ohne das schwarze Loch im Vordergrund, FOV=120°×60°:





Schwarzschild: M=1, a=0, ℧=0 → ℳ=1 - einfaches schwarzes Loch, Beobachter auf r=50GM/c², ϑ=egal (kugelsymmetrisch), FOV=120°×60°:





Reissner-Nordström: M=1, a=0, ℧=1 → ℳ=½ - geladenes SL, Beobachter auf r=50GM/c², ϑ=egal (kugelsymmetrisch), FOV=120°×60°:





Kerr-Newman: M=1, a=√½, ℧=√½ → ℳ=√(3/8) - Geladen und rotierend, Beobachter auf r=50GM/c², ϑ=90° (äquatorial), FOV=120°×60°:





Kerr-Newman: M=1, a=√½, ℧=√½ → ℳ=√(3/8) - Geladen und rotierend, Beobachter auf r=50GM/c², ϑ=45° (schräg), FOV=120°×60°:





Kerr-Newman: M=1, a=√½, ℧=√½ → ℳ=√(3/8) - Geladen und rotierend, Beobachter auf r=50GM/c², ϑ=1° (polar), FOV=120°×60°:





Kerr: M=1, a=1, ℧=0 → ℳ=√½ - Beobachter auf r=50GM/c², ϑ=90° (äquatorial), FOV=120°×60°:





Kerr: M=1, a=1, ℧=0 → ℳ=√½ - Beobachter auf r=50GM/c², ϑ=45° (schräg), FOV=120°×60°:





Kerr: M=1, a=1, ℧=0 → ℳ=√½ - Beobachter auf r=50GM/c², ϑ=1° (polare Ausrichtung, Rotation gegen den Uhrzeigersinn), FOV=120°×60°:





Einheiten: G=M=c=k=1. M=Massenäquivalent der Gesamtenergie, ℳ=irreduzible Masse, a=Spinparameter, ℧=Ladung G=Gravitationskonstante, k=Coulombkonstante, c=Lichtgeschwindigkeit. Die Gesamtenergie Mc² des schwarzen Lochs ist in allen obigen Szenarien gleich, aber unterschiedlich auf irreduzible Masse, Rotationsenergie und elektromagnetische Feldenergie aufgeteilt.



Polare und äquatoriale Ansicht eines Kerr-SL mit a=0.99, im unteren Bild mit eingeblendeten Horizonten und Ergosphären (ein Beobachter würde diese natürlich nicht sehen). Die θ-Ausrichtung des SL läuft von 1° (polare Ansicht) bis 90° (äquatoriale Ansicht):



Hochaufgelöste Videos: https://pewtube.com/user/Yukterez_Net

[Yukterez]

Rohmaterial: Panoramaphoto geschossen von Masato OHTA (heiwa4126), Source: Commons, Ort: Tokyo. Assumptions: alles im Panorama ist sehr viel weiter vom Beobachter entfernt als das schwarze Loch. Hier wird nur der Gravitationslinseneffekt gezeigt, nicht aber die Zerstörung die ein SL in so einem Szenario anrichten würde.


SL: M=1, a=1, ℧=0. FOV=360°×180°, ungelinstes Panorama:



Zoom auf die Stelle an der das schwarze Loch platziert wird, Ф=0°, ϑ=0°, FOV=103.2°×61.4°:



Das schwarze Loch ist r=50GM/c² vom Beobachter entfernt. Position, Fokus und Zoom auf Ф=0°, ϑ=0°, FOV=103.2°×61.4°:



Entfernung vom SL: 20GM/c². Position, Fokus und Zoom auf Ф=0°, ϑ=0°, FOV=103.2°×61.4°:





Blick in die entgegengesetzte Richtung, FOV=360°×180°:



Zoom auf die Stelle an der das schwarze Loch platziert wird, Ф=0°, ϑ=180°, FOV=103.2°×61.4°:



Entfernung vom SL: 50GM/c². Position, Fokus und Zoom auf Ф=180°, ϑ=0°, FOV=103.2°×61.4°:



Entfernung vom SL: 20GM/c². Position, Fokus und Zoom auf Ф=180°, ϑ=0°, FOV=103.2°×61.4°:





Hier wird das schwarze Loch auf halber Höhe über der Skyline platziert; Ansicht ohne SL, FOV=103.2°×61.4°:



Entfernung vom SL: 50GM/c². Position, Fokus und Zoom auf Ф=180°, ϑ=-45°, FOV=103.2°×61.4°:



Entfernung vom SL: 20GM/c². Position, Fokus und Zoom auf Ф=180°, ϑ=-45°, FOV=103.2°×61.4°:





Video: 360° Korotation eines ZAMO (bei r=20 benötigt ein voller Umlauf 8022πGM/c³ Eigenzeit), Video: ogg & mp4, FOV=103.2°×61.4°:





Ausrichtung des SL, relative Perspektive des Beobachters: Edge On (äquatoriale Ansicht: ϑ=90°). Die linke Seite des SL rotiert auf den Betrachter zu, und die rechte von ihm weg. Siehe auch den dazugehörigen Beitrag im Uwudl-Forum.

[Yukterez]

Rohmaterial: Double Negative's Interstellar Plattkarte. Auf den folgenden Bildern befindet sich der Beobachter auf r=50GM/c² vom mit a/M=1 rotierenden SL entfernt. Der Sichtwinkel ist immer FOV=151.2°×75.6°.


Ungelinst, flache Raumzeit (Minkowski):



Mit schwarzem Loch, lokale Geschwindigkeit v=0 (Inertialsystem eines ZAMO):



Radialgeschwindigkeit v=c/2 auf das SL zu:



Radialgeschwindigkeit v=c/2 vom SL weg:



Azimutalgeschwindigkeit v=-c/2 retrograd:



Azimutalgeschwindigkeit v=+c/2 prograd:



Polargeschwindigkeit v=-c/2:



Polargeschwindigkeit v=+c/2::



Gesamtgeschwindigkeit v=√(vr²+vθ²+vφ²)=c/2 (Bewegung in alle 3 Richtungen gleichzeitig):



Im letzten Bild bewegt sich der Beobachter radial auf das SL zu während er eine retrograde Axial- und eine in Richtung Südpol gehende Polargeschwindigkeit hat. Bei vorhandener Transversalbewegung sieht der bewegte Beobachter das Objekt nicht dort wo es sich tatsächlich befindet, sondern in die andere Richtung als die in die er selber sich bewegt verschoben (das schwarze Loch befindet sich in allen Bildern genau in der Mitte).




Kerr-Newman SL mit a=0.7, ℧=0.7, Beobachter auf r=5. Sichtwinkel: 360°×180° Vollpanorama im Plattkartenformat:


Geschwindigkeit: lokal ruhender ZAMO, v=0 (v relativ zu den Fixsternen: 0.0666)



Axialgeschwindigkeit: vφ=-0.45



Axialgeschwindigkeit: vφ=-0.675



Axialgeschwindigkeit: vφ=-0.8



Axialgeschwindigkeit: vφ=-0.9



Axialgeschwindigkeit: vφ=-0.999 (retrograd)



Axialgeschwindigkeit: vφ=+0.999 (prograd)



Radialgeschwindigkeit: vr=-0.78273 (Freifallgeschwindigkeit)



Radialgeschwindigkeit: vr=+0.78273 (Fluchtgeschwindigkeit)

[Yukterez]

Die folgenden überextremen Lösungen beschreiben keine schwarzen Löcher, und kommen in der Natur vermutlich nicht so vor. Da sie aber dennoch von theoretischem Interesse sind bekommen sie ebenfalls eine Gallerie. Wie im oberen Beitrag mit den schwarzen Löchern ist auch hier die Gesamtenergie Mc² der nackten Singularitäten in allen Beispielen gleich, aber verschieden aufgeteilt.

Standbilder (a²+℧²=2², zum Vergrößern auf die Bilder klicken)

relativistic raytracing of the shadow and the gravitational lensing of naked singularities
Reissner-Nordström: M=1, a=0, ℧=2 → ℳ=½+i√¾ - nackte Singularität. Beobachter: r=50GM/c², ϑ=egal (kugelsymmetrisch), FOV=120°×60°:





Kerr-Newman: M=1, a=√2, ℧=√2 → ℳ=∜(3/16)·(1+i) - nackte Singularität. Beobachter: r=50GM/c², ϑ=90° (äquatorial), FOV=120°×60°:





Kerr-Newman: M=1, a=√2, ℧=√2 → ℳ=∜(3/16)·(1+i) - nackte Singularität. Beobachter: r=50GM/c², ϑ=45° (schräg), FOV=120°×60°:





Kerr-Newman: M=1, a=√2, ℧=√2 → ℳ=∜(3/16)·(1+i) - nackte Singularität. Beobachter: r=50GM/c², ϑ=1° (polar), FOV=120°×60°:





Kerr: M=1, a=2, ℧=0 → ℳ=√¾+i/2 - nackte Singularität. Beobachter: r=50GM/c², ϑ=90° (äquatorial), FOV=120°×60°:





Kerr: M=1, a=2, ℧=0 → ℳ=√¾+i/2 - nackte Singularität. Beobachter: r=50GM/c², ϑ=45° (schräg), FOV=120°×60°:





Kerr: M=1, a=2, ℧=0 → ℳ=√¾+i/2 - nackte Singularität. Beobachter: r=50GM/c², ϑ=1° (polar), FOV=120°×60°:

Videos: äquatoriale Orbits (60 fps, a²+℧²=1.01²)


Reissner Nordström, M=1, a=0, ℧=1.01 - nackte Singularität. Beobachter: r=20GM/c², ϑ=90°, FOV=154.8°×77.4°:





Kerr, M=1, a=1.01, ℧=0 - nackte Singularität. Beobachter: r=20GM/c², ϑ=90°, FOV=154.8°×77.4°:





Kerr Newman, M=1, a=0.7141778, ℧=0.7141778 - nackte Singularität. Beobachter: r=20GM/c², ϑ=90°, FOV=154.8°×77.4°:





Die Videos sind auch in Full HD verfügbar (dann brauchen sie aber auch entsprechend länger um zu laden): pewtube.com/user/Yukterez_Net. Vergleich: Waseda, S.14 und DeVries, S.20

[Yukterez]

Probe: Vergleich mit einem unabhängigen Code
Vorlage aus Kip Thorne et al's Lösungsheft arXiv:1502.03808 für (a) v=0.546c, (b) v=0 (ZAMO) und (c) v=-0.481c (stationär) bei a=0.999:





Vergleich mit dem entsprechenden Output des hiesigen Raytracers; nachempfundenes Hintergrundmuster:



(a) Prograde Orbitalgeschwindigkeit, v=0.546c auf r=2.6GM/c²:



(b) Lokal drehimpulsfreier ZAMO, v=0 auf r=2.6GM/c²:



(c) Relativ zu den Fixsternen stationärer Beobachter mit negativer Frame Dragging Geschwindigkeit, v=-0.481c auf r=2.6GM/c²:



Bei den obigen Bildern handelt es sich um 360°×180° Vollpanoramas im Plattkartenformat.

[Yukterez]

Mathematischer Anhang (im Aufbau):

1) Geschwindigkeitskomponenten der beim Betrachter eintreffenden Strahlen


Das schwarze Loch befindet sich im lokalen Koordinatensystem des Beobachters auf {x,z}={0,0}. Der {x,y,z}-Vektor eines geradewegs von vorne kommenden Photons ist daher



Rotation entlang der x-Achse:



Rotation entlang der z-Achse:



Normierter Geschwindigkeitsvektor der eintreffenden Strahlen im euklidschen Referenzsystem:



Übersetzung in den lokalen Tetrad eines ZAMO:



Normierung:



Aberration bei lokaler (relativ zum ZAMO) Beobachtergeschwindigkeit {ur,uθ,uФ}:



mit den Komponenten



Radiale Komponente:



Polare Komponente:



Azimutale Komponente:



Das sind die Endbedingungen der im Auge des Betrachters auf der {α,β}-Sichtebene eintreffenden Strahlen. Diese werden hernach rückwärts integriert bis sie die Ebene des Interesses, in dem Fall ein auf R>>r aufgespanntes Kugelpanorama, treffen.


2) Geodäten der Photonen


Zusammenhang der lokalen Geschwindigkeitskomponenten mit den ersten Ableitungen:



Zweite Ableitungen:



mit den Christoffelsymbolen



Für die aufsummierten Terme in expliziter Form siehe den Faden über die Kerr Newman Metrik.


3) Bildtransformation


Nun wird numerisch nach der Eigenzeit bzw. dem affinen Parameter bei dem das Photon in der Region des Interesses (in dem Fall auf der auf R aufgespannten Kugelschale) war gesolved:



Die für τ berechneten {Ф,θ} Koordinaten des Lichtstrahls werden dann vom Rohmaterial im Plattkartenformat auf die lokale Sichtebene gemappt:




4) Technische Details


An sich steht bereits alles im Code, aber eine Übersetzung ins Latex folgt dennoch demnächst (in den nächsten Wochen).

5) Nächster Schritt


Akkretionsscheibe (Mapping des Rohmaterials auf die {r,Ф} statt auf die {Ф,θ} Ebene)

All work here as well as the 3rd party source material (the raw panorama photos) is published under the Creative Commons License. Simon Tyran (Yukterez), 2018