Raytracer (Accretion Disk)

(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
(* > raytracing.yukterez.net | 29.04.2018 - 20.06.2019 | Version 7T | Simon Tyran, Vienna *)
(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)

kernels = 5;                                                          (* Parallelisierung *)
grain   = 4;                            (* Subparallelisierung auf kernels*grain Streifen *)
breite  = 120;                                               (* Zielabmessungen in Pixeln *)
hoehe   = 120;          (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
zoom    = 15;                                 (* doppelter Zoom ergibt halben Sichtwinkel *)
round   = 0;      (* 0: undurchsichtig, 1: Vorderseite, 2: Rückseite, 3 und 4: Lichtechos *)

LaunchKernels[kernels]
wp = MachinePrecision;                                                     (* Genauigkeit *)

pic = Import["http://yukterez.net/mw/scheibe.png"]                (* Scheibentextur laden *)
snc = 1;    (* Scheibensynchronisierung: 1 = statisch, 2|3 = Empfang|Emission, 4 = vOrbit *)

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

r0   = 100;                                           (* Radialkoordinate des Beobachters *)
si   = isco;                                             (* Akkretionsscheibe Innenradius *)
sr   = 7;                                                (* Akkretionsscheibe Außenradius *)
θ0   = 70 π/180;                                                           (* Breitengrad *)
φ0   = 0;                                                                   (* Längengrad *)

tmax =-3 r0;                                            (* zeitlicher Integrationsbereich *)

a    = 0.7;                                                              (* Spinparameter *)
℧    = 0.7;                                     (* spezifische Ladung des schwarzen Lochs *)
v0   = 1;                                                  (* Geschwindigkeit 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 *)

hvs  = 0;                                            (* horizontaler Versatz in Radianten *)
vvs  = 0;                                              (* vertikaler Versatz in Radianten *)

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;

vθ  =-vϑ;
θs  = π/2;
θi  =-θ0+π;
q   = 0;
rA  = 1+Sqrt[1-a^2-℧^2];
rE  = 1+Sqrt[1-℧^2];

rp = rf/.Solve[4 a^2 (rf-℧^2)-(rf^2-3 rf+2 ℧^2)^2 == 0 && rf >= 1, rf];
rP = 1.01 Min[rp]; Rp = 1.01 Max[rp];
isco = Quiet[Min[RI/.NSolve[0 == RI (6 RI-RI^2-9 ℧^2+3 a^2)+4 ℧^2 (℧^2-a^2)-8 a (RI-℧^2)^(3/2) && RI>=If[Element[rA, Reals], rA, 0], RI]]];
{"r horizon" -> N@rA, "r ergosphere" -> N@rE, "r isco" -> N@isco, "r photon pro" -> N@Min[rp], "r photon ret" -> N@Max[rp], "r disk" -> N@sr, "r observer" -> N@r0, "θ observer" -> N@θ0}

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];

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[Table[{rr, If[Quiet@NumericQ[red[rr][[1, 2]]], red[rr][[1, 2]], 0]}, {rr, 0, sr, 0.02}]];
φs = Interpolation[Table[{rr, If[Quiet@NumericQ[red[rr][[2, 2]]], red[rr][[2, 2]], 0]}, {rr, 0, sr, 0.02}]];
ts = Interpolation[Table[{rr, If[Quiet@NumericQ[red[rr][[3, 2]]], red[rr][[3, 2]], 0]}, {rr, 0, sr, 0.02}]];

plot[func_, label_] := Plot[func, {rr, rP, sr}, GridLines -> {{Min[rp], Max[rp], rA, si, isco, rE, sr}, {}},
Frame -> True, ImagePadding -> {{40,1}, {12,1}}, ImageSize -> 340, PlotLabel -> label, PlotRange->{{0,sr}, Automatic}]

plot[Sqrt[Χs[rr]/Δs[rr]/Σs[rr]],  "Gravitational time dilation: y=dt/dт, x=r"]
plot[ts[rr],  "Total time dilation: y=dt/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[φs[rr]/ts[rr], "Shapirodelayed angular velocity: y=dφ/dt, x=r"]
plot[φs[rr],  "Coordinate speed: y=dφ/dτ, x=r"]
plot[vs[rr],  "Local velocity: y=v=dl/dτ, x=r"]

й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, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
DGL, sol, εj, pθi, pr0, Q, k, t10, r10, Θ10, Φ10, т0, т1, R1, t, r, θ, φ, τ, plunge, plunge2,
X, Y, Z, tt, rt, θt, φt, т, ξ, stepsize, laststep, mtl, mta, ft, fτ, varb, Σj, Δj, Χj,
ωj, dθ, dφ, dphi, dtheta, drad, dtau, vrj, vθj, vφj, vtj, vφt, shf},

vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
                                 (* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
vr0a = vw[[3]];
vφ0a = vw[[2]];
vθia = vw[[1]];
                                                                                (* Betrag *)
vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2];
                                                                            (* Normierung *)
vr0n = vr0a/vt0a;
vφ0n = vφ0a/vt0a;
vθ0n = vθia/vt0a;
                                              (* Relativistische Geschwindigkeitsaddition *)
V={vr0n, vθ0n, vφ0n};
W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V);
                                                                            (* Aberration *)
vr0i = W[[1]];
vθ0i = W[[2]];
vφ0i = W[[3]];

rnd = If[round == 2,
{WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0,
(plunge=τ) && (tt=If[snc==1, 0, t[τ]]) && (rt=r[τ]) && (θt=θ[τ]) &&
(φt=φ[τ]) && (dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) &&
(dphi=φ'[τ]); "StopIntegration"]},
If[round == 1,
{WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0,
(plunge=τ) && (tt=If[snc==1, 0, t[τ]]) && (rt=r[τ]) && (θt=θ[τ]) &&
(φt=φ[τ]) && (dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) &&
(dphi=φ'[τ]); "StopIntegration"]},
If[round == 4,
{WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0,
(plunge2=τ)],
WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0
&&τ<plunge2-0.1,(plunge=τ)&&(tt=If[snc==1,0,t[τ]])&&(rt=r[τ])&&(θt=θ[τ])&&(φt=φ[τ]) &&
(dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) && (dphi=φ'[τ]); "StopIntegration"]},
If[round == 3,
{WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0,
(plunge2=τ)],
WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0
&&τ<plunge2-0.1,(plunge=τ)&&(tt=If[snc==1,0,t[τ]])&&(rt=r[τ])&&(θt=θ[τ])&&(φt=φ[τ]) &&
(dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) && (dphi=φ'[τ]); "StopIntegration"]},
{WhenEvent[Mod[θ[τ],π]==Pi/2.0 && r[τ]>si && r[τ]<sr, (plunge=τ) &&
(tt=If[snc==1, 0, t[τ]]) && (rt=r[τ]) && (θt=θ[τ]) && (φt=φ[τ]) && (dtau=t'[τ]) &&
(drad=r'[τ]) && (dtheta=θ'[τ]) && (dphi=φ'[τ]); "StopIntegration"]}
]]]];

DGL = {                                               (* Kerr Newman Bewegungsgleichungen *)

t''[τ]==-(((r'[τ] ((a^2+r[τ]^2) (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+
2 (-℧^2+r[τ]) t'[τ]))+a (2 a^4 Cos[θ[τ]]^2+a^2 ℧^2 (3+Cos[2 θ[τ]]) r[τ]-
a^2 (3+Cos[2 θ[τ]]) r[τ]^2+4 ℧^2 r[τ]^3-6 r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))/(a^2+℧^2+(-2+
r[τ]) r[τ])+a^2 θ'[τ] (Sin[2 θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])-2 a Cos[θ[τ]] (℧^2-
2 r[τ]) Sin[θ[τ]]^3 φ'[τ]))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^2),

t'[0]==ς,
t[0]==0,

r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
(a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),

r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θi]^2)/(a^2+(-2+r0) r0+℧^2)],
r[0]==r0,

θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
10 a^4 r[τ]+11 a^4 r[τ]^2-8 a^2 ℧^2 r[τ]^2+16 a^2 r[τ]^3+16 a^2 r[τ]^4+8 r[τ]^6+
a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),

θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2],
θ[0]==θi,

φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),

φ'[0]==(vφ0i ((-2+r0) r0+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+
a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+℧^2) Sin[θi]^2)]+(a vφ0i (2 r0-
℧^2) Sin[θi])/((r0^2+a^2 Cos[θi]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θi]^2)),
φ[0]==φ0,

rnd};
                                                                            (* Integrator *)
sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
WorkingPrecision-> wp,
MaxSteps-> Infinity,
InterpolationOrder-> All];

Σj = rt^2+a^2 Cos[θt]^2;
Δj = rt^2-2 rt+a^2+℧^2;
Χj = (rt^2+a^2)^2-a^2 Sin[θt]^2 Δj;
ωj = (a(2 rt-℧^2))/Χj;

т0 = If[NumericQ[plunge], plunge, tmax];
т1 = If[NumericQ[plunge], tt, 0];

dφ = If[snc==1, 0, If[snc==2, т1 ωj, If[snc==3, (т1+r0) ωj, (т1+r0) φs[rt]/ts[rt]]]];

ft=If[т0>tmax+1,If[rt>sr,{π,sr},If[rt<si,{π,sr},{φt-dφ,rt}]],{π,sr}];
ft]];

(* Plot[{
raytracer[{u, 0}][[1]],
raytracer[{u, 0}][[2]]
},{u, -π/2, π/2},
PlotPoints -> 40,
PlotStyle -> {Red, Blue},
Frame->True] *)

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[pic,raytracer,{breite,Ceiling[hoehe/kernels/grain]},
DataRange->{{0,2π-2π/width},{si,sr+(sr-si)/height}},
PlotRange->{{-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom},
Padding->"Periodic"],
{x, 0, π-π/kernels/grain, π/kernels/grain}]

image  = ImageAssemble[{Table[{img[[x]]}, {x, kernels grain, 1, -1}]}]