Yukterez

Raytracer (Horizon Surface)

Jun 9th, 2019 (edited)
145
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* > raytracing.yukterez.net | 29.04.2018 - 09.06.2019 | Version 7H | Simon Tyran, Vienna *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. kernels = 5; (* Parallelisierung *)
  6. grain = 4; (* Subparallelisierung auf kernels*grain Streifen *)
  7. breite = 120; (* Zielabmessungen in Pixeln *)
  8. hoehe = 120; (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
  9. zoom = 15; (* doppelter Zoom ergibt halben Sichtwinkel *)
  10.  
  11. LaunchKernels[kernels]
  12. wp = MachinePrecision; (* Genauigkeit *)
  13. (* Planetenoberfläche laden *)
  14. pic =Import["https://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/Equirectangular-projection.jpg/320px-Equirectangular-projection.jpg"];
  15. snc = 3; (* Kugelsynchronisierung: 1 = statisch, 2 = bei Empfang, 3 = bei Emission *)
  16.  
  17. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  18. (* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
  19. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  20.  
  21. r0 = 100; (* Radialkoordinate des Beobachters *)
  22. θ0 = 70 π/180; (* Breitengrad *)
  23. φ0 = π/4; (* Längengrad *)
  24.  
  25. R0 = 1.01 rA; (* Radius des umspannenden Kugelschalenpanoramas *)
  26. tmax =-5/3 r0; (* zeitlicher Integrationsbereich *)
  27.  
  28. a = 0.7; (* Spinparameter *)
  29. ℧ = 0.7; (* spezifische Ladung des schwarzen Lochs *)
  30. v0 = 1; (* Geschwindigkeit des Photons *)
  31.  
  32. vr = 0; (* Radiale Geschwindigkeit des Beobachters *)
  33. vϑ = 0; (* Polare Geschwindigkeit des Beobachters *)
  34. vφ = 0; (* Azimutale Geschwindigkeit des Beobachters: 0 für ZAMO, -й0 für stationär *)
  35.  
  36. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  37. (* 2) Metrische Koeffizienten und Formeln ||||||||||||||||||||||||||||||||||||||||||||||| *)
  38. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  39.  
  40. Σ = r0^2+a^2 Cos[θ0]^2;
  41. Δ = r0^2-2 r0+a^2+℧^2;
  42. Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
  43. μ = 0; q = 0;
  44.  
  45. gtt = (2r0-℧^2)/Σ-1;
  46. grr = Σ/Δ;
  47. gθθ = Σ;
  48. gφφ = Χ/Σ Sin[θ0]^2;
  49. gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;
  50.  
  51. vθ =-vϑ;
  52. θi =-θ0+π;
  53. rA = 1+Sqrt[1-a^2-℧^2];
  54. й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-
  55. 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-
  56. 2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]);
  57.  
  58. U={+vr, +vθ, +vφ};
  59. γ=1/Sqrt[1-Norm[U]^2];
  60.  
  61. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  62. (* 3) Rotationsmatrix für die auf der Sichtebene eintreffenden Strahlen ||||||||||||||||| *)
  63. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  64.  
  65. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
  66. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  67. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  68.  
  69. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  70. (* 4) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  71. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  72.  
  73. raytracer[{Ф_, ϑ_}] :=
  74.  
  75. Quiet[Module[{V, W, vw, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
  76. DGL, sol, ε, Lz, pθi, pr0, Q, k, t10, r10, Θ10, Φ10, т0, т1, R1, t, r, θ, φ, τ, plunge,
  77. X, Y, Z, tt, rt, θt, φt, т, ξ, stepsize, laststep, mtl, mta, ft, fτ, Σj, Δj, Χj, ωj, dθ, dφ},
  78.  
  79. vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
  80. (* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
  81. vr0a = vw[[3]];
  82. vφ0a = vw[[2]];
  83. vθia = vw[[1]];
  84. (* Betrag *)
  85. vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2];
  86. (* Normierung *)
  87. vr0n = vr0a/vt0a;
  88. vφ0n = vφ0a/vt0a;
  89. vθ0n = vθia/vt0a;
  90. (* Relativistische Geschwindigkeitsaddition *)
  91. V={vr0n, vθ0n, vφ0n};
  92. W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V);
  93. (* Aberration *)
  94. vr0i = W[[1]];
  95. vθ0i = W[[2]];
  96. vφ0i = W[[3]];
  97. (* Integrationsende *)
  98. mtl={"EventLocator", "Event"->If[(r[τ]==R0||r[τ]<R0) == True, 0, 1]};
  99. mta=mtl;
  100.  
  101. DGL = { (* Kerr Newman Bewegungsgleichungen *)
  102.  
  103. t''[τ]==-(((r'[τ] ((a^2+r[τ]^2) (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+
  104. 2 (-℧^2+r[τ]) t'[τ]))+a (2 a^4 Cos[θ[τ]]^2+a^2 ℧^2 (3+Cos[2 θ[τ]]) r[τ]-
  105. a^2 (3+Cos[2 θ[τ]]) r[τ]^2+4 ℧^2 r[τ]^3-6 r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))/(a^2+℧^2+(-2+
  106. r[τ]) r[τ])+a^2 θ'[τ] (Sin[2 θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])-2 a Cos[θ[τ]] (℧^2-
  107. 2 r[τ]) Sin[θ[τ]]^3 φ'[τ]))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^2),
  108.  
  109. t'[0]==ς,
  110. t[0]==0,
  111.  
  112. r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
  113. (a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
  114. r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
  115. r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
  116. 8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
  117. Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
  118. 8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
  119.  
  120. r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θi]^2)/(a^2+(-2+r0) r0+℧^2)],
  121. r[0]==r0,
  122.  
  123. θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
  124. r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
  125. r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
  126. r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
  127. 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+
  128. a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
  129. r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
  130.  
  131. θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2],
  132. θ[0]==θi,
  133.  
  134. φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
  135. 8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
  136. Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
  137. 16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
  138. θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
  139. 2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),
  140.  
  141. φ'[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+
  142. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+
  143. 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-
  144. ℧^2) Sin[θi])/((r0^2+a^2 Cos[θi]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
  145. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θi]^2)),
  146. φ[0]==φ0,
  147.  
  148. WhenEvent[r[τ]==R0||r[τ]<R0,
  149. (plunge=τ) && (tt=If[snc==1, 0, t[τ]]) && (rt=r[τ]) && (θt=θ[τ]) && (φt=φ[τ]);
  150. "StopIntegration"]
  151. };
  152. (* Integrator *)
  153. sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
  154. WorkingPrecision-> wp,
  155. MaxSteps-> Infinity,
  156. InterpolationOrder-> All];
  157.  
  158. Σj = rt^2+a^2 Cos[θt]^2;
  159. Δj = rt^2-2 rt+a^2+℧^2;
  160. Χj = (rt^2+a^2)^2-a^2 Sin[θt]^2 Δj;
  161. ωj = (a(2 rt-℧^2))/Χj;
  162.  
  163. т0 = If[NumericQ[plunge], plunge, tmax];
  164. т1 = If[NumericQ[plunge], tt, 0];
  165.  
  166. dφ = If[snc==1, 0, If[snc==2, т1 ωj, If[snc==3, (т1+r0) ωj, -(т1+r0) ωj]]];
  167.  
  168. т0 = If[NumericQ[plunge], plunge, tmax];
  169. R1 = rt;
  170. (* Berechung der Ursprungskoordinaten der Photonen *)
  171. If[т0<tmax+1, {0, -π/2}, If[R1<rA, {0, -π/2}, If[R1>4 R0, {0, -π/2}, {φt-dφ, θt+π/2}]]]]]
  172.  
  173. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  174. (* 5) Testbild laden und transformieren ||||||||||||||||||||||||||||||||||||||||||||||||| *)
  175. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  176.  
  177. fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]}
  178.  
  179. pcr = ParallelTable[
  180. ImageTransformation[pic, fpt, DataRange->{{-1, 1}, {0, 1}},
  181. PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Padding->"Periodic"],
  182. {x, 1, 2 kernels}];
  183. pct = ImageAssemble[Table[{pcr[[x]]}, {x, 2 kernels, 1, -1}]];
  184.  
  185. width = ImageDimensions[pic][[1]]; height = ImageDimensions[pic][[2]];
  186. hzoom = If[breite>2 hoehe, 1/zoom, 1/zoom/2/hoehe*breite];
  187. vzoom = If[breite>2 hoehe, 1/zoom*2 hoehe/breite, 1/zoom];
  188.  
  189. FOV -> {360.0 hzoom "degree", 180.0 vzoom "degree"}
  190.  
  191. img = ParallelTable[
  192. ImageTransformation[pct, raytracer, {breite, Ceiling[hoehe/kernels/grain]},
  193. DataRange->{{-π, π-2π/width}, {-π/2, 3π/2}},
  194. PlotRange->{{-π, π} hzoom, {-π/2+x, -π/2+x+π/kernels/grain} vzoom},
  195. Padding->"Periodic"],
  196. {x, 0, π-π/kernels/grain, π/kernels/grain}]
  197.  
  198. source = ImageResize[pic, {2 hoehe, hoehe}]
  199. image = ImageAssemble[Table[{img[[x]]}, {x, kernels grain, 1, -1}]]
RAW Paste Data