Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
- (* > 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}]}]
RAW Paste Data