Yukterez

Raytracer (Old Version)

Apr 7th, 2018 (edited)
344
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* > raytracing.yukterez.net | 29.04.2018 - 12.04.2019 | Version 7F | Simon Tyran, Vienna *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. kernels = 6; (* Parallelisierung *)
  6. grain = 18; (* Subparallelisierung auf kernels*grain Streifen *)
  7. breite = 216; (* Zielabmessungen in Pixeln *)
  8. hoehe = 108; (* 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"]; (* Hintergrundpanorama laden *)
  15.  
  16. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  17. (* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
  18. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  19.  
  20. r0 = 5; (* Radialkoordinate des Beobachters *)
  21. θ0 = π/2; (* Breitengrad *)
  22. φ0 = 0; (* Längengrad *)
  23.  
  24. R0 = 500; (* Radius des umspannenden Kugelschalenpanoramas *)
  25. tmax =-5/3 R0; (* zeitlicher Integrationsbereich *)
  26.  
  27. a = 0.7; (* Spinparameter *)
  28. ℧ = 0.7; (* spezifische Ladung des schwarzen Lochs *)
  29. v0 = 1; (* Geschwindigkeit des Photons *)
  30.  
  31. vr = 0; (* Radiale Geschwindigkeit des Beobachters *)
  32. vθ = 0; (* Polare Geschwindigkeit des Beobachters *)
  33. vφ = 0; (* Azimutale Geschwindigkeit des Beobachters: 0 für ZAMO, -й0 für stationär *)
  34.  
  35. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  36. (* 2) Metrische Koeffizienten und Formeln ||||||||||||||||||||||||||||||||||||||||||||||| *)
  37. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  38.  
  39. Σ = r0^2+a^2 Cos[θ0]^2;
  40. Δ = r0^2-2 r0+a^2+℧^2;
  41. Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
  42. μ = 0; q = 0;
  43.  
  44. gtt = (2r0-℧^2)/Σ-1;
  45. grr = Σ/Δ;
  46. gθθ = Σ;
  47. gφφ = Χ/Σ Sin[θ0]^2;
  48. gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;
  49.  
  50. θi =-θ0+π;
  51. rA = 1+Sqrt[1-a^2-℧^2];
  52. й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-
  53. 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-
  54. 2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]);
  55.  
  56. U={+vr, +vθ, +vφ};
  57. γ=1/Sqrt[1-Norm[U]^2];
  58.  
  59. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  60. (* 3) Rotationsmatrix für die auf der Sichtebene eintreffenden Strahlen ||||||||||||||||| *)
  61. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  62.  
  63. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
  64. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  65. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  66.  
  67. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  68. (* 4) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  69. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  70.  
  71. raytracer[{Ф_, ϑ_}] :=
  72.  
  73. Quiet[Module[{V, W, vw, tMax, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
  74. DGL, sol, ε, Lz, pθi, pr0, Q, k, t10, r10, Θ10, Φ10, т0, R1, t, r, θ, φ, τ, plunge,
  75. X, Y, Z, rt, θt, φt, т, ξ, stepsize, laststep, mtl, mta, ft, fτ},
  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. (* Integrationsende *)
  96. mtl=If[a^2+℧^2>1,
  97. {"EventLocator", "Event"->r[τ]-R0-1.0},
  98. {"EventLocator", "Event"->If[(r[τ]==1.01rA || r[τ] == R0+1.0) == 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]==-((a (2 r0-℧^2) Sin[θi]^2 (vφ0i (-2 r0+r0^2+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+
  110. r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-
  111. 2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)]+
  112. (a vφ0i (2 r0-℧^2) Sin[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+
  113. a^2 Cos[θi]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2))))/((a^2-2 r0+r0^2+
  114. ℧^2) (r0^2+a^2 Cos[θi]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θi]^2)))+\[Sqrt](((vr0i^2 (a^2-2 r0+
  115. r0^2+℧^2)+vθ0i^2 (a^2-2 r0+r0^2+℧^2)) (r0^2+a^2 Cos[θi]^2) (-2 r0+r0^2+℧^2+a^2 Cos[θi]^2)+
  116. (a^2 (-2 r0+℧^2)^2 Sin[θi]^4 (vφ0i (-2 r0+r0^2+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+r0^2)^2-
  117. a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-2 r0+
  118. r0^2+℧^2) (r0^2+a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)]+
  119. (a vφ0i (2 r0-℧^2) Sin[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+
  120. a^2 Cos[θi]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)))^2)/((a^2-2 r0+r0^2+
  121. ℧^2) (r0^2+a^2 Cos[θi]^2)^2)-((2 r0-r0^2-℧^2-a^2 Cos[θi]^2) Sin[θi]^2 ((a^2+r0^2)^2-
  122. a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2) (vφ0i (-2 r0+r0^2+a^2 Cos[θi]^2) Csc[θi] Sqrt[((a^2+
  123. r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2-
  124. 2 r0+r0^2+℧^2) (r0^2+a^2 Cos[θi]^2))/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)]+
  125. (a vφ0i (2 r0-℧^2) Sin[θi] Sqrt[((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)/(r0^2+
  126. a^2 Cos[θi]^2)])/((a^2+r0^2)^2-a^2 (a^2-2 r0+r0^2+℧^2) Sin[θi]^2)))^2)/((a^2-2 r0+r0^2+
  127. ℧^2) (r0^2+a^2 Cos[θi]^2)^2))/((a^2-2 r0+r0^2+℧^2) (-2 r0+r0^2+℧^2+a^2 Cos[θi]^2)^2)),
  128. t[0]==0,
  129.  
  130. r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
  131. (a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
  132. r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
  133. r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
  134. 8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
  135. Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
  136. 8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
  137.  
  138. r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θi]^2)/(a^2+(-2+r0) r0+℧^2)],
  139. r[0]==r0,
  140.  
  141. θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
  142. r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
  143. r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
  144. r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
  145. 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+
  146. a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
  147. r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
  148.  
  149. θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2],
  150. θ[0]==θi,
  151.  
  152. φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
  153. 8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
  154. Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
  155. 16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
  156. θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
  157. 2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),
  158.  
  159. φ'[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+
  160. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+
  161. 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-
  162. ℧^2) Sin[θi])/((r0^2+a^2 Cos[θi]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
  163. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θi]^2)),
  164. φ[0]==φ0
  165.  
  166. };
  167. (* Integrator *)
  168. sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
  169. WorkingPrecision-> wp,
  170. MaxSteps-> Infinity,
  171. Method-> mta,
  172. InterpolationOrder-> All,
  173. StepMonitor :> (laststep=plunge; plunge=τ;
  174. stepsize=plunge-laststep;), Method->{"EventLocator",
  175. "Event" :> (If[stepsize<2*^-2, 0, 1])}];
  176. (* Integrationszeit *)
  177. tMax = Max[tmax, plunge+1/10];
  178. (* Raumkoordinaten *)
  179. rt[τ_] := Evaluate[r[τ]/.sol][[1]];
  180. θt[τ_] := Evaluate[θ[τ]/.sol][[1]]+π/2;
  181. φt[τ_] :=-Evaluate[φ[τ]/.sol][[1]]-π;
  182. (* Affiner Parameter bei Emission *)
  183. т[coord_, dist_] := ξ/.FindRoot[coord[ξ]-dist, {ξ, tMax 9/10, tMax, -1}];
  184. т0 = т[rt, R0];
  185. R1 = rt[т0]; (* Check ob die Photonen von der Hemisphäre kommen *)
  186. (* Berechung der Ursprungskoordinaten der Photonen *)
  187. If[т0>0, {0, -π/2}, If[R1<10r0, {0, -π/2}, If[R1>10R0, {0, -π/2}, {φt[т0], θt[т0]}]]]]]
  188.  
  189. mem : raytrace[{Ф_, ϑ_}] := mem = raytracer[{Ф, ϑ}]
  190.  
  191. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  192. (* 5) Testbild laden und transformieren ||||||||||||||||||||||||||||||||||||||||||||||||| *)
  193. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  194.  
  195. fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]}
  196.  
  197. pcr = ParallelTable[
  198. ImageTransformation[pic, fpt, DataRange->{{-1, 1}, {0, 1}},
  199. PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Padding->"Periodic"],
  200. {x, 1, 2 kernels}];
  201. pct = ImageAssemble[Table[{pcr[[x]]}, {x, 2 kernels, 1, -1}]];
  202.  
  203. width = ImageDimensions[pic][[1]]; height = ImageDimensions[pic][[2]];
  204. hzoom = If[breite>2 hoehe, 1/zoom, 1/zoom/2/hoehe*breite];
  205. vzoom = If[breite>2 hoehe, 1/zoom*2 hoehe/breite, 1/zoom];
  206.  
  207. FOV -> {360.0 hzoom "degree", 180.0 vzoom "degree"}
  208.  
  209. img = ParallelTable[
  210. ImageTransformation[pct, raytrace, {breite, Ceiling[hoehe/kernels/grain]},
  211. DataRange->{{-π, π-2π/width}, {-π/2, 3π/2}},
  212. PlotRange->{{-π, π-2π/width} hzoom, {-π/2+x, -π/2+x+π/kernels/grain} vzoom},
  213. Padding->"Periodic"],
  214. {x, 0, π-π/kernels/grain, π/kernels/grain}]
  215.  
  216. source = ImageResize[pic, {2 hoehe, hoehe}]
  217. image = ImageAssemble[Table[{img[[x]]}, {x, kernels grain, 1, -1}]]
RAW Paste Data