SHARE
TWEET

Raytracer (Old Version)

Yukterez Apr 7th, 2018 (edited) 258 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. ClearAll["Global`*"]; ClearAll["Local`*"];
  6. Needs["DifferentialEquations`NDSolveProblems`"];
  7. Needs["DifferentialEquations`NDSolveUtilities`"];
  8.                                        
  9. kernels = 6;                                                          (* Parallelisierung *)
  10. grain   = 18;                           (* Subparallelisierung auf kernels*grain Streifen *)
  11. breite  = 216;                                               (* Zielabmessungen in Pixeln *)
  12. hoehe   = 108;          (* Höhe sollte ein ganzzahliges Vielfaches von kernels*grain sein *)
  13. zoom    = 1;                                  (* doppelter Zoom ergibt halben Sichtwinkel *)
  14.  
  15. LaunchKernels[kernels]
  16. wp = MachinePrecision;                                                     (* Genauigkeit *)
  17.  
  18. pic = Import["http://yukterez.net/mw/2/flip90.png"];         (* Hintergrundpanorama laden *)
  19.  
  20. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  21. (* 1) Startbedingungen und Position des Beobachters ||||||||||||||||||||||||||||||||||||| *)
  22. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  23.  
  24. r0   = 5;                                             (* Radialkoordinate des Beobachters *)
  25. θ0   = π/2;                                                                (* Breitengrad *)
  26. φ0   = 0;                                                                   (* Längengrad *)
  27.  
  28. R0   = 500;                              (* Radius des umspannenden Kugelschalenpanoramas *)
  29. tmax =-5/3 R0;                                          (* zeitlicher Integrationsbereich *)
  30.  
  31. a    = 0.7;                                                              (* Spinparameter *)
  32. ℧    = 0.7;                                     (* spezifische Ladung des schwarzen Lochs *)
  33. v0   = 1;                                                  (* Geschwindigkeit des Photons *)
  34.  
  35. vr   = 0;                                      (* Radiale Geschwindigkeit des Beobachters *)
  36. vθ   = 0;                                       (* Polare Geschwindigkeit des Beobachters *)
  37. vφ   = 0;     (* Azimutale Geschwindigkeit des Beobachters: 0 für ZAMO, -й0 für stationär *)
  38.  
  39. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  40. (* 2) Metrische Koeffizienten und Formeln ||||||||||||||||||||||||||||||||||||||||||||||| *)
  41. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  42.  
  43. Σ = r0^2+a^2 Cos[θ0]^2;
  44. Δ = r0^2-2 r0+a^2+℧^2;
  45. Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
  46. μ = 0; q = 0;
  47.  
  48. gtt = (2r0-℧^2)/Σ-1;
  49. grr = Σ/Δ;
  50. gθθ = Σ;
  51. gφφ = Χ/Σ Sin[θ0]^2;
  52. gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;
  53.  
  54. θi  =-θ0+π;
  55. rA  = 1+Sqrt[1-a^2-℧^2];
  56. й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-
  57. 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-
  58. 2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]);
  59.  
  60. U={+vr, +vθ, +vφ};
  61. γ=1/Sqrt[1-Norm[U]^2];
  62.  
  63. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  64. (* 3) Rotationsmatrix für die auf der Sichtebene eintreffenden Strahlen ||||||||||||||||| *)
  65. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  66.  
  67. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
  68. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  69. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  70.  
  71. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  72. (* 4) Raytracing Funktionscontainer ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  73. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  74.  
  75. raytracer[{Ф_, ϑ_}] :=
  76.  
  77. Quiet[Module[{V, W, vw, tMax, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
  78. DGL, sol, ε, Lz, pθi, pr0, Q, k, t10, r10, Θ10, Φ10, т0, R1, t, r, θ, φ, τ, plunge,
  79. X, Y, Z, rt, θt, φt, т, ξ, stepsize, laststep, mtl, mta, ft, fτ},
  80.  
  81. vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
  82.                                  (* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
  83. vr0a = vw[[3]] Sqrt[grr];
  84. vφ0a = vw[[2]] Sqrt[gφφ]/r0/Sin[θi];
  85. vθia = vw[[1]] Sqrt[gθθ]/r0;
  86.                                                                                 (* Betrag *)
  87. vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2];
  88.                                                                             (* Normierung *)
  89. vr0n = vr0a/vt0a;
  90. vφ0n = vφ0a/vt0a;
  91. vθ0n = vθia/vt0a;
  92.                                               (* Relativistische Geschwindigkeitsaddition *)
  93. V={vr0n, vθ0n, vφ0n};
  94. W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V);
  95.                                                                             (* Aberration *)
  96. vr0i = W[[1]];
  97. vθ0i = W[[2]];
  98. vφ0i = W[[3]];
  99.                                                                       (* Integrationsende *)                        
  100. mtl=If[a^2+℧^2>1,
  101. {"EventLocator", "Event"->r[τ]-R0-1.0},
  102. {"EventLocator", "Event"->If[(r[τ]==1.01rA || r[τ] == R0+1.0) == True, 0, 1]}];
  103. mta=mtl;
  104.  
  105. DGL = {                                               (* Kerr Newman Bewegungsgleichungen *)
  106.  
  107. t''[τ]==-(((r'[τ] ((a^2+r[τ]^2) (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+
  108. 2 (-℧^2+r[τ]) t'[τ]))+a (2 a^4 Cos[θ[τ]]^2+a^2 ℧^2 (3+Cos[2 θ[τ]]) r[τ]-
  109. a^2 (3+Cos[2 θ[τ]]) r[τ]^2+4 ℧^2 r[τ]^3-6 r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))/(a^2+℧^2+(-2+
  110. r[τ]) r[τ])+a^2 θ'[τ] (Sin[2 θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])-2 a Cos[θ[τ]] (℧^2-
  111. 2 r[τ]) Sin[θ[τ]]^3 φ'[τ]))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^2),
  112.  
  113. 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+
  114. 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-
  115. 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)]+
  116. (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+
  117. 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+
  118. ℧^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+
  119. 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)+
  120. (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-
  121. 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+
  122. 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)]+
  123. (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+
  124. 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+
  125. ℧^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-
  126. 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+
  127. 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-
  128. 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)]+
  129. (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+
  130. 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+
  131. ℧^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)),
  132. t[0]==0,
  133.  
  134. r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
  135. (a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
  136. r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
  137. r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
  138. 8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
  139. Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
  140. 8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
  141.  
  142. r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θi]^2)/(a^2+(-2+r0) r0+℧^2)],
  143. r[0]==r0,
  144.  
  145. θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
  146. r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
  147. r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
  148. r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
  149. 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+
  150. a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
  151. r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
  152.  
  153. θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2],
  154. θ[0]==θi,
  155.  
  156. φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
  157. 8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
  158. Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
  159. 16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
  160. θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
  161. 2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),
  162.  
  163. φ'[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+
  164. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+
  165. 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-
  166. ℧^2) Sin[θi])/((r0^2+a^2 Cos[θi]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
  167. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θi]^2)),
  168. φ[0]==φ0
  169.  
  170. };
  171.                                                                             (* Integrator *)
  172. sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
  173. WorkingPrecision-> wp,
  174. MaxSteps-> Infinity,
  175. Method-> mta,
  176. InterpolationOrder-> All,
  177. StepMonitor :> (laststep=plunge; plunge=τ;
  178. stepsize=plunge-laststep;), Method->{"EventLocator",
  179. "Event" :> (If[stepsize<2*^-2, 0, 1])}];
  180.                                                                       (* Integrationszeit *)
  181. tMax = Max[tmax, plunge+1/10];
  182.                                                                        (* Raumkoordinaten *)
  183. rt[τ_] := Evaluate[r[τ]/.sol][[1]];
  184. θt[τ_] := Evaluate[θ[τ]/.sol][[1]]+π/2;
  185. φt[τ_] :=-Evaluate[φ[τ]/.sol][[1]]-π;
  186.                                                         (* Affiner Parameter bei Emission *)
  187. т[coord_, dist_] := ξ/.FindRoot[coord[ξ]-dist, {ξ, tMax 9/10, tMax, -1}];
  188. т0 = т[rt, R0];
  189. R1 = rt[т0];                           (* Check ob die Photonen von der Hemisphäre kommen *)
  190.                                        (* Berechung der Ursprungskoordinaten der Photonen *)
  191. If[т0>0, {0, -π/2}, If[R1<10r0, {0, -π/2}, If[R1>10R0, {0, -π/2}, {φt[т0], θt[т0]}]]]]]
  192.  
  193. mem : raytrace[{Ф_, ϑ_}] := mem = raytracer[{Ф, ϑ}]
  194.  
  195. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  196. (* 5) Testbild laden und transformieren ||||||||||||||||||||||||||||||||||||||||||||||||| *)
  197. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  198.  
  199. fpt[{x_, y_}] := {If[y<0, x+1, x], If[y<0, -y, y]}
  200.  
  201. pcr = ParallelTable[
  202. ImageTransformation[pic, fpt, DataRange->{{-1, 1}, {0, 1}},
  203. PlotRange->{{-1, 1}, {-1+(x-1)/kernels, -1+x/kernels}}, Padding->"Periodic"],
  204. {x, 1, 2 kernels}];
  205. pct = ImageAssemble[{Table[{pcr[[x]]}, {x, 2 kernels, 1, -1}]}];
  206.  
  207. width = ImageDimensions[pic][[1]]; height = ImageDimensions[pic][[2]];
  208. hzoom = If[breite>2 hoehe, 1/zoom, 1/zoom/2/hoehe*breite];
  209. vzoom = If[breite>2 hoehe, 1/zoom*2 hoehe/breite, 1/zoom];
  210.  
  211. FOV -> {360.0 hzoom "degree", 180.0 vzoom "degree"}
  212.  
  213. img = ParallelTable[
  214. ImageTransformation[pct, raytrace, {breite, Ceiling[hoehe/kernels/grain]},
  215. DataRange->{{-π, π-2π/width}, {-π/2, 3π/2}},
  216. PlotRange->{{-π, π-2π/width} hzoom, {-π/2+x, -π/2+x+π/kernels/grain} vzoom},
  217. Padding->"Periodic"],
  218. {x, 0, π-π/kernels/grain, π/kernels/grain}]
  219.  
  220. source = ImageResize[pic, {2 hoehe, hoehe}]
  221. image  = ImageAssemble[{Table[{img[[x]]}, {x, kernels grain, 1, -1}]}]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top