(* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* > raytracing.yukterez.net | 07.04.2018 - 09.09.2020 | Version 20 | Simon Tyran, Vienna *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) Pause[1] (* Minkowski Koordinaten *) wp = MachinePrecision; mt0 = Automatic; mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}}; mta = mt0; (* mt0 für Geschwindigkeit, mt1 für Genauigkeit *) kernels = 6; (* Parallelisierung *) grain = 5; (* Subparallelisierung auf kernels*grain Streifen *) rsp = "Nearest"; (* Resampling *) breite = 240; (* Zielabmessungen in Pixeln *) hoehe = 120; (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *) zoom = 1; (* doppelter Zoom ergibt halben Sichtwinkel *) LaunchKernels[kernels] wp = MachinePrecision; (* Genauigkeit *) pic1 = Import["http://www.yukterez.net/mw/1/flip70.png"]; (* Hintergrundpanorama laden *) pic2 = Import["http://www.yukterez.net/mw/1/worldmap.png"]; (* Sphärenoberfläche laden *) pic3 = Import["http://www.yukterez.net/mw/akkretionsscheibe.jpg"];(* Scheibentextur laden *) pic4 = Import["http://www.yukterez.net/mw/disk.png"]; (* Scheibengeometrie laden *) pic7 = Import["http://www.yukterez.net/mw/1/bw.png"]; (* Checkerboard laden *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) (* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *) (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *) r0 = 10; (* Radialkoordinate des Beobachters *) R0 = 500; (* Radius des umspannenden Kugelschalenpanoramas *) R1 = 1+1/(5 Sqrt[2]); (* Sphäre *) si = 1.6367824; (* Akkretionsscheibe Innenradius *) sr = 7; (* Akkretionsscheibe Außenradius *) θ0 = 70 π/180; (* Breitengrad *) φ0 = 0; (* Längengrad *) t0 = 0; (* Eigenzeit des Beobachters *) tmax =-3 R0; (* zeitlicher Integrationsbereich *) vr = 0; (* Radiale Geschwindigkeit des Beobachters *) vϑ = 0; (* Polare Geschwindigkeit des Beobachters *) vφ = 0; (* Azimutale Geschwindigkeit des Beobachters *) vL = Sqrt[vr^2+vϑ^2+vφ^2]; hvs = 0; (* ArcSin[vφ/vL] *) (* horizontaler Versatz in Radianten *) vvs = 0; (* ArcSin[vϑ/vL] *) (* vertikaler Versatz in Radianten *) fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]} pcr1 = ParallelTable[ ImageTransformation[pic1, fpt, DataRange->{{-1, 1}, {0, 1}}, PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Padding->"Periodic"], {x, 1, 2 kernels}]; pct1 = ImageAssemble[Table[{pcr1[[x]]}, {x, 2 kernels, 1, -1}]]; pcr2 = ParallelTable[ ImageTransformation[pic2, fpt, DataRange->{{-1, 1}, {0, 1}}, PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Padding->"Periodic"], {x, 1, 2 kernels}]; pct2 = ImageAssemble[Table[{pcr2[[x]]}, {x, 2 kernels, 1, -1}]]; pcr3 = ParallelTable[ ImageTransformation[pic7, fpt, DataRange->{{-1, 1}, {0, 1}}, PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Resampling->rsp, Padding->"Periodic"], {x, 1, 2 kernels}]; pct3 = ImageAssemble[Table[{pcr3[[x]]}, {x, 2 kernels, 1, -1}]]; a = 0; ℧ = 0; v0 = 1; q = 0; st = 0.01; vϑ = 0; vθ =-vϑ; θs = π/2; θi =-θ0+π; θ1 = θi; dφv[rt_, tt_] := 0; shf[rt_, δt_, δr_, δθ_, δφ_] := 1; ωφ = -vφ Csc[θ1]/r0; ωθ = -vϑ/r0; gtt = -1; grr = +1; gθθ = +r0^2; gφφ = +r0^2 Sin[θ1]^2; gtφ = +0; ς = 1; j[v_] := Sqrt[1-v^2]; nq[x_] := If[NumericQ[x], If[Element[x, Reals], x, 0], 0]; U = {+vr, +vθ, +vφ}; γ = 1/Sqrt[1-Norm[U]^2]; Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z}; xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]}; xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]}; raytracer[{Ф_, ϑ_}] := Quiet[Module[{DGL, sol, εj, pθi, pr0, Q, k, V, W, vw, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a, t10, r10, Θ10, Φ10, t, r, θ, φ, τ, т, т0, т1, т2, т3, т4, т5, plunge, plunge0, plunge1, plunge2, plunge3, plunge4, plunge5, plunge6, plunge7, dθ0, dφ0, δφ0, δθ0, δr0, δt0, tt0, rt0, θt0, φt0, dθ1, dφ1, δφ1, δθ1, δr1, δt1, tt1, rt1, θt1, φt1, dθ2, dφ2, δφ2, δθ2, δr2, δt2, tt2, rt2, θt2, φt2, dθ3, dφ3, δφ3, δθ3, δr3, δt3, tt3, rt3, θt3, φt3, dθ4, dφ4, δφ4, δθ4, δr4, δt4, tt4, rt4, θt4, φt4, dθ5, dφ5, δφ5, δθ5, δr5, δt5, tt5, rt5, θt5, φt5, dθ6, dφ6, δφ6, δθ6, δr6, δt6, tt6, rt6, θt6, φt6, X, Y, Z, ξ, stepsize, laststep, mtl, ft, fτ, varb, ft0s, ft1s, ft2s, ft3s, ft4s, ft0v, ft1v, ft2v, ft3v, ft4v, ft5v, ft0f, ft1f, ft2f, ft3f, ft4f, ft5h, ft5b, ft5f}, vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2]; (* Übersetzung des Einfallswinkels in den lokalen Tetrad *) vr0a = vw[[3]]; vφ0a = vw[[2]]; vθia = vw[[1]]; (* Betrag *) vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2]; (* Normierung *) vr0n = vr0a/vt0a; vφ0n = vφ0a/vt0a; vθ0n = vθia/vt0a; (* Relativistische Geschwindigkeitsaddition *) V={vr0n, vθ0n, vφ0n}; W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V); (* Aberration *) vr0i = W[[1]]; vθ0i = W[[2]]; vφ0i = W[[3]]; DGL = { (* Kerr Newman Bewegungsgleichungen *) t''[τ] == 0, t'[0] == 1, t[0] == 0, r''[τ] == r[τ](θ'[τ]^2+Sin[θ[τ]]^2 φ'[τ]^2), r'[0] == vr0i, r[0] == r0, θ''[τ] == Sin[θ[τ]] Cos[θ[τ]] φ'[τ]^2-2 θ'[τ] r'[τ]/r[τ], θ'[0] == vθ0i/r0, θ[0] == θi, φ''[τ] == -2 φ'[τ] (r'[τ]+r[τ] θ'[τ] Cot[θ[τ]])/r[τ], φ'[0] == vφ0i Csc[θ1]/r0, φ[0] == φ0, WhenEvent[Mod[θ[τ], π] == π/2.0 && r[τ]>si && r[τ]si && r[τ]0 && NumericQ[plunge1] == False, (plunge1=τ) && (tt1=t[τ]) && (rt1=r[τ]) && (θt1=θ[τ]) && (φt1=φ[τ]) && (δt1=t'[τ]) && (δr1=r'[τ]) && (δθ1=θ'[τ]) && (δφ1=φ'[τ])], WhenEvent[Mod[θ[τ], π] == π/2.0 && r[τ]>si && r[τ]si && r[τ]0 && τ < plunge1-0.05 && NumericQ[plunge3]==False && NumericQ[plunge1] == True, (plunge3=τ) && (tt3=t[τ]) && (rt3=r[τ]) && (θt3=θ[τ]) && (φt3=φ[τ]) && (δt3=t'[τ]) && (δr3=r'[τ]) && (δθ3=θ'[τ]) && (δφ3=φ'[τ])], WhenEvent[Mod[θ[τ], π] == π/2.0 && r[τ]>si && r[τ]R1, "StopIntegration"]], WhenEvent[Abs[r[τ]] == R0||r[τ]>R0 && NumericQ[plunge6] == False, (plunge6=τ) && (tt6=t[τ]) && (rt6=r[τ]) && (θt6=θ[τ]) && (φt6=φ[τ]); "StopIntegration"] }; (* Integrator *) sol = TimeConstrained[NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax}, WorkingPrecision-> wp, Method-> mta, MaxSteps-> Infinity, StepMonitor :> (laststep=plunge7; plunge7=τ; stepsize=plunge7-laststep;), Method->{"EventLocator", "Event" :> (If[stepsize<1*^-4, 0, 1])} ], 10]; ft0s=If[NumericQ[plunge0], If[rt0>sr, {π, sr}, If[rt0sr, {π, sr}, If[rt1sr, {π, sr}, If[rt2sr, {π, sr}, If[rt3sr, {π, sr}, If[rt4sr, {0, 0}, If[rt0sr, {0, 0}, If[rt1sr, {0, 0}, If[rt2sr, {0, 0}, If[rt3sr, {0, 0}, If[rt4sr, {0, 0}, If[rt0sr, {0, 0}, If[rt1sr, {0, 0}, If[rt2sr, {0, 0}, If[rt3sr, {0, 0}, If[rt44 R0, {0, -π/2}, {φt6-π, θt6+π/2}]], {0, -π/2}]; { ft0s, ft1s, ft2s, ft3s, ft4s, ft0s+ft0v, ft1s+ft1v, ft2s+ft2v, ft3s+ft3v, ft4s+ft4v, ft0f, ft1f, ft2f, ft3f, ft4f, ft5h, ft5b }]]; mem : raytrace[{Ф_, ϑ_}] := mem = Quiet[Re[raytracer[{Ф, ϑ}]]]; ray01[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[01]]; ray02[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[02]]; ray03[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[03]]; ray04[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[04]]; ray05[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[05]]; ray06[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[06]]; ray07[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[07]]; ray08[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[08]]; ray09[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[09]]; ray10[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[10]]; ray11[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[11]]; ray12[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[12]]; ray13[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[13]]; ray14[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[14]]; ray15[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[15]]; ray16[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[16]]; ray17[{Ф_, ϑ_}] := raytrace[{Ф, ϑ}][[17]]; width1 = ImageDimensions[pic1][[1]]; height1 = ImageDimensions[pic1][[2]]; width2 = ImageDimensions[pic2][[1]]; height2 = ImageDimensions[pic2][[2]]; width3 = ImageDimensions[pic3][[1]]; height3 = ImageDimensions[pic3][[2]]; width4 = ImageDimensions[pic4][[1]]; height4 = ImageDimensions[pic4][[2]]; width5 = ImageDimensions[pic7][[1]]; height5 = ImageDimensions[pic7][[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]; "Geschätzte Rechenzeit" -> 1.2 (AbsoluteTiming[Do[raytracer[{ RandomReal[{-π, π}]/zoom, RandomReal[{-π/2, π/2}]/zoom }], {ü, 1, 50}]][[1]])/50*hoehe*breite/kernels/60 "Minuten" FOV->{360.0 hzoom "degree", 180.0 vzoom "degree"} img = ParallelTable[{ (* 1 Hintergrundpanorama *) ImageTransformation[pct1, ray17, {breite, Ceiling[hoehe/kernels/grain]}, DataRange->{ {-π, π-2π/width1}, {-π/2, 3π/2} }, PlotRange->{ {-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+x+vvs/vzoom, -π/2+x+vvs/vzoom+π/kernels/grain} vzoom }, Resampling->rsp, Padding->"Periodic"], (* 2 Sphäre *) ImageTransformation[pct2, ray16, {breite, Ceiling[hoehe/kernels/grain]}, DataRange->{ {-π, π-2π/width2}, {-π/2, 3π/2} }, PlotRange->{ {-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+x+vvs/vzoom, -π/2+x+π/kernels/grain+vvs/vzoom} vzoom }, Resampling->rsp, Padding->"Periodic"], (* 3 Scheibe Textur *) ImageTransformation[pic3, ray01, {breite, Ceiling[hoehe/kernels/grain]}, DataRange->{ {0, 2π-2π/width3}, {si, sr+(sr-si)/height3} }, PlotRange->{ {-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom }, Resampling->rsp, Padding->"Periodic"], (* 4 Scheibe Geometrie *) ImageTransformation[pic4, ray01, {breite, Ceiling[hoehe/kernels/grain]}, DataRange->{ {0, 2π-2π/width4}, {si, sr} }, PlotRange->{ {-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom }, Resampling->rsp, Padding->"Periodic"], (* 5 Hintergrundpanorama *) ImageTransformation[pct3, ray17, {breite, Ceiling[hoehe/kernels/grain]}, DataRange->{ {-π, π-2π/width5}, {-π/2, 3π/2} }, PlotRange->{ {-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+x+vvs/vzoom, -π/2+x+vvs/vzoom+π/kernels/grain} vzoom }, Resampling->rsp, Padding->"Periodic"], (* 6 Sphäre *) ImageTransformation[pct3, ray16, {breite, Ceiling[hoehe/kernels/grain]}, DataRange->{ {-π, π-2π/width5}, {-π/2, 3π/2} }, PlotRange->{ {-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+x+vvs/vzoom, -π/2+x+π/kernels/grain+vvs/vzoom} vzoom }, Resampling->rsp, Padding->"Periodic"], (* 7 Scheibe Textur komplett *) ImageTransformation[pic7, ray01, {breite, Ceiling[hoehe/kernels/grain]}, DataRange->{ {0, 2π-2π/width5}, {si, sr+(sr-si)/height5} }, PlotRange->{ {-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom }, Resampling->rsp, Padding->"Periodic"] }, {x, 0, π-π/kernels/grain, π/kernels/grain}]; image[num_] := ImageAssemble[Table[{img[[x, num]]}, {x, kernels grain, 1, -1}]]; fi[x_] := ColorNegate[x]; Grid[{Table[image[n], {n, 1, 7, 1}]}] Quit[]