Yukterez

Raindrop Raytracer (Background Sky)

May 22nd, 2020 (edited)
95
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* > raytracing.yukterez.net | 23.05.2020 | Raindrop coordinates, B | Simon Tyran, Vienna *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. kernels = 6; (* Parallelisierung *)
  6. grain = 10; (* Subparallelisierung auf kernels*grain Streifen *)
  7. breite = 240; (* Zielabmessungen in Pixeln *)
  8. hoehe = 120; (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
  9. zoom = 1; (* doppelter Zoom ergibt halben Sichtwinkel *)
  10.  
  11. LaunchKernels[kernels]
  12. wp = MachinePrecision; (* Genauigkeit *)
  13.  
  14. pic = Import["http://yukterez.net/mw/2/flip90.png"];
  15.  
  16. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  17. (* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
  18. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  19.  
  20. r0 = 10; (* Radialkoordinate des Beobachters *)
  21. θ0 = π/2; (* Breitengrad *)
  22. φ0 = 0; (* Längengrad *)
  23.  
  24. R0 = 50; (* Radius des umspannenden Kugelschalenpanoramas *)
  25. tmax =-50 R0; (* zeitlicher Integrationsbereich *)
  26.  
  27. v0 = 1; (* Geschwindigkeit des Photons *)
  28.  
  29. vr = Sqrt[2/r0]; (* Radiale Geschwindigkeit des Beobachters, relativ zum Raindrop *)
  30. vϑ = 0; (* Polare Geschwindigkeit des Beobachters *)
  31. vφ = 0; (* Azimutale Geschwindigkeit des Beobachters *)
  32.  
  33. hvs = 0; (* horizontaler Versatz in Radianten *)
  34. vvs = 0; (* vertikaler Versatz in Radianten *)
  35.  
  36. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  37. (* 2) Metrische Koeffizienten und Formeln ||||||||||||||||||||||||||||||||||||||||||||||| *)
  38. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  39.  
  40. gtt = 2/r0-1;
  41. grr = 1;
  42. gθθ = r0^2;
  43. gφφ = r0^2 Sin[θ0]^2;
  44. gtr = Sqrt[2 r0^3]/r0^2;
  45.  
  46. ς = Sqrt[-gtt];
  47. j[v_] := Sqrt[1-v^2];
  48.  
  49. vθ =-vϑ;
  50. θi =-θ0+π;
  51. rA = 2;
  52.  
  53. U={+vr, +vθ, +vφ};
  54. γ=1/Sqrt[1-Norm[U]^2];
  55. μ=0;
  56.  
  57. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  58. (* 3) Rotationsmatrix für die auf der Sichtebene eintreffenden Strahlen ||||||||||||||||| *)
  59. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  60.  
  61. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
  62. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  63. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  64.  
  65. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  66. (* 4) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  67. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  68.  
  69. raytracer[{Ф_, ϑ_}] :=
  70.  
  71. Quiet[Module[{V, W, vw, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
  72. DGL, sol, ε, Lz, pθi, pr0, Q, k, t10, r10, Θ10, Φ10,
  73. т0, т1, R1, t, r, θ, φ, τ, plunge, plunge2,
  74. X, Y, Z, rt, θt, φt, т, ξ, stepsize, laststep, mtl, mta, ft, fτ,
  75. dt0, dr0, dθ0, dφ0, initcon},
  76.  
  77. vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
  78. (* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
  79. vr0a = vw[[3]] Sqrt[grr];
  80. vφ0a = vw[[2]] Sqrt[gφφ]/r0/Sin[θi];
  81. vθia = vw[[1]] Sqrt[gθθ]/r0;
  82. (* Betrag *)
  83. vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2];
  84. (* Normierung *)
  85. vr0n = vr0a/vt0a;
  86. vφ0n = vφ0a/vt0a;
  87. vθ0n = vθia/vt0a;
  88. (* Relativistische Geschwindigkeitsaddition *)
  89. V={vr0n, vθ0n, vφ0n};
  90. W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V);
  91. (* Aberration *)
  92. vr0i = W[[1]];
  93. vθ0i = W[[2]];
  94. vφ0i = W[[3]];
  95.  
  96. initcon = NSolve[
  97. dr0 == vr0i-Sqrt[2 vr0i^2+
  98. 2 r0^2 (vθ0i/r0)^2+
  99. r0^2 (vφ0i/Sqrt[r0^2 Sin[θ0]^2])^2-
  100. Cos[2 θ0] r0^2 (vφ0i/Sqrt[
  101. r0^2 Sin[θ0]^2])^2]/Sqrt[r0]
  102. &&
  103. dθ0 == vθ0i/r0
  104. &&
  105. dφ0 == vφ0i/Sqrt[r0^2 Sin[θ0]^2]
  106. &&
  107. -μ ==
  108. N[-(((0^2+2 r0^2+0^2 Cos[2 θ0]) (dr0)^2)/(2 (0^2+r0^2)))-(2 Sqrt[2 r0-
  109. 0^2] dr0 dt0)/Sqrt[0^2+r0^2]+(1+(-4 r0+2 0^2)/(0^2+2 r0^2+0^2 Cos[2 θ0])) (dt0)^2+
  110. (-r0^2-0^2 Cos[θ0]^2) dθ0^2+(2 0 Sqrt[2 r0-0^2] Sin[θ0]^2 dr0 dφ0)/Sqrt[0^2+
  111. r0^2]+(2 0 (2 r0-0^2) Sin[θ0]^2 dt0 dφ0)/(r0^2+0^2 Cos[θ0]^2)+((-(0^2+r0^2)^2 Sin[θ0]^2+
  112. 0^2 (0^2+(-2+r0) r0+0^2) Sin[θ0]^4) (dφ0)^2)/(r0^2+0^2 Cos[θ0]^2)]
  113. &&
  114. dt0 == 1,
  115. {dθ0, dr0, dt0, dφ0}, Reals];
  116.  
  117. DGL = { (* Kerr Newman Bewegungsgleichungen *)
  118.  
  119. t''[τ] == -(Derivative[1][r][τ]^2/(Sqrt[2] Sqrt[r[τ]^3]))-
  120. (2 Derivative[1][r][τ] Derivative[1][t][τ])/r[τ]^2+
  121. (Sqrt[2] Sqrt[
  122. r[τ]^3] (-Derivative[1][t][τ]^2+
  123. r[τ]^3 (Derivative[1][θ][τ]^2+
  124. Sin[θ[τ]]^2 Derivative[1][φ][τ]^2)))/r[τ]^4,
  125.  
  126. t'[0] == dt0/.initcon[[1]],
  127. t[0] == 0,
  128.  
  129. r''[τ]==(1/(r[τ]^4))(r[τ]^2 Derivative[1][r][τ]^2+
  130. 2 Sqrt[2] Sqrt[r[τ]^3]
  131. Derivative[1][r][τ] Derivative[1][t][τ]+
  132. (-2+r[τ]) r[τ] (-Derivative[1][t][τ]^2+
  133. r[τ]^3 (Derivative[1][θ][τ]^2+
  134. Sin[θ[τ]]^2 Derivative[1][φ][τ]^2))),
  135.  
  136. r'[0] == dr0/.initcon[[1]],
  137. r[0] == r0,
  138.  
  139. θ''[τ] == -((2 Derivative[1][r][τ] Derivative[1][θ][τ])/
  140. r[τ])+Cos[θ[τ]] Sin[θ[τ]] Derivative[1][φ][τ]^2,
  141.  
  142. θ'[0] == dθ0/.initcon[[1]],
  143. θ[0] == θi,
  144.  
  145. φ''[τ] == -((2 (Derivative[1][r][τ]+
  146. Cot[θ[τ]] r[τ] Derivative[1][θ][τ]) Derivative[1][φ][τ])/r[τ]),
  147.  
  148. φ'[0] == dφ0/.initcon[[1]],
  149. φ[0] == φ0,
  150.  
  151. WhenEvent[r[τ]>1.0 R0,
  152. (plunge=τ) && (rt=r[τ]) && (θt=θ[τ]) && (φt=φ[τ]);"StopIntegration"]
  153.  
  154. };
  155. (* Integrator *)
  156. sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
  157. WorkingPrecision-> wp,
  158. MaxSteps-> Infinity,
  159. InterpolationOrder-> All];
  160. (* Affiner Parameter bei Emission *)
  161. т0 = plunge;
  162. R1 = rt;
  163. (* Berechung der Ursprungskoordinaten der Photonen *)
  164. If[NumericQ[plunge],
  165. If[т0>0, {0, -π/2},
  166. If[R1<r0, {0, -π/2},
  167. If[R1>5 R0, {0, -π/2},
  168. {φt-π, θt+π/2}]]],
  169. {0, -π/2}]]]
  170.  
  171. mem : raytrace[{Ф_, ϑ_}] := mem = raytracer[{Ф, ϑ}]
  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->{{-π+hvs/hzoom, π+hvs/hzoom} hzoom,
  195. {-π/2+x+vvs/vzoom, -π/2+x+vvs/vzoom+π/kernels/grain} vzoom},
  196. Padding->"Periodic"],
  197. {x, 0, π-π/kernels/grain, π/kernels/grain}]
  198.  
  199. image = ImageAssemble[Table[{img[[x]]}, {x, kernels grain, 1, -1}]]
  200.  
  201. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  202. (* 6) Alternative initcon ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  203. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  204.  
  205. initcon = NSolve[
  206. dr0 == vr0i-Sqrt[2 vr0i^2+
  207. 2 r0^2 (vθ0i/r0)^2+
  208. r0^2 (vφ0i/Sqrt[r0^2 Sin[θ0]^2])^2-
  209. Cos[2 θ0] r0^2 (vφ0i/Sqrt[
  210. r0^2 Sin[θ0]^2])^2]/Sqrt[r0]
  211. &&
  212. dθ0 == vθ0i/r0
  213. &&
  214. dφ0 == vφ0i/Sqrt[r0^2 Sin[θ0]^2]
  215. &&
  216. -μ ==
  217. N[-(((0^2+2 r0^2+0^2 Cos[2 θ0]) (dr0)^2)/(2 (0^2+r0^2)))-(2 Sqrt[2 r0-
  218. 0^2] dr0 dt0)/Sqrt[0^2+r0^2]+(1+(-4 r0+2 0^2)/(0^2+2 r0^2+0^2 Cos[2 θ0])) (dt0)^2+
  219. (-r0^2-0^2 Cos[θ0]^2) dθ0^2+(2 0 Sqrt[2 r0-0^2] Sin[θ0]^2 dr0 dφ0)/Sqrt[0^2+
  220. r0^2]+(2 0 (2 r0-0^2) Sin[θ0]^2 dt0 dφ0)/(r0^2+0^2 Cos[θ0]^2)+((-(0^2+r0^2)^2 Sin[θ0]^2+
  221. 0^2 (0^2+(-2+r0) r0+0^2) Sin[θ0]^4) (dφ0)^2)/(r0^2+0^2 Cos[θ0]^2)]
  222. &&
  223. dt0 > 0,
  224. {dθ0, dr0, dt0, dφ0}, Reals];
RAW Paste Data