Yukterez

Raytracer (Redshift)

Jun 19th, 2019 (edited)
196
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* > raytracing.yukterez.net | 29.04.2018 - 20.06.2019 | Version 7B | 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. round = 0; (* 0: undurchsichtig, 1: Vorderseite, 2: Rückseite, 3 und 4: Lichtechos *)
  11.  
  12. LaunchKernels[kernels]
  13. wp = MachinePrecision; (* Genauigkeit *)
  14.  
  15. pic = Import["http://yukterez.net/mw/gradient4.png"]; (* Gradienten laden *)
  16. snc = 1; (* Scheibensynchronisierung: 1 = statisch, 2|3 = Empfang|Emission, 4 = vOrbit *)
  17.  
  18. r0 = 100; (* Radialkoordinate des Beobachters *)
  19. si = isco; (* Akkretionsscheibe Innenradius *)
  20. sr = 7; (* Akkretionsscheibe Außenradius *)
  21. θ0 = 70 π/180; (* Breitengrad *)
  22. φ0 = 0; (* Längengrad *)
  23.  
  24. tmax =-3 r0; (* zeitlicher Integrationsbereich *)
  25.  
  26. a = 0.7; (* Spinparameter *)
  27. ℧ = 0.7; (* spezifische Ladung des schwarzen Lochs *)
  28. v0 = 1; (* Geschwindigkeit des Photons *)
  29.  
  30. vr = 0; (* Radiale Geschwindigkeit des Beobachters *)
  31. vϑ = 0; (* Polare Geschwindigkeit des Beobachters *)
  32. vφ = 0; (* Azimutale Geschwindigkeit des Beobachters: 0 für ZAMO, -й0 für stationär *)
  33.  
  34. hvs = 0; (* horizontaler Versatz in Radianten *)
  35. vvs = 0; (* vertikaler Versatz in Radianten *)
  36.  
  37. Σ = r0^2+a^2 Cos[θ0]^2;
  38. Δ = r0^2-2 r0+a^2+℧^2;
  39. Χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 Δ;
  40. Σs[rs_] := rs^2;
  41. Δs[rs_] := rs^2-2 rs+a^2+℧^2;
  42. Χs[rs_] := (rs^2+a^2)^2-a^2 Δs[rs];
  43. κs[rs_] := a;
  44.  
  45. vθ =-vϑ;
  46. θs = π/2;
  47. θi =-θ0+π;
  48. q = 0;
  49. rA = 1+Sqrt[1-a^2-℧^2];
  50. rE = 1+Sqrt[1-℧^2];
  51. ς = Sqrt[Χ/Δ/Σ];
  52.  
  53. rp = rf/.Solve[4 a^2 (rf-℧^2)-(rf^2-3 rf+2 ℧^2)^2 == 0 && rf >= 1, rf];
  54. rP = 1.01 Min[rp]; Rp = 1.01 Max[rp];
  55. isco = Quiet[Min[RI/.NSolve[0 == RI (6 RI-RI^2-9 ℧^2+3 a^2)+4 ℧^2 (℧^2-a^2)-8 a (RI-℧^2)^(3/2) && RI>=If[Element[rA, Reals], rA, 0], RI]]];
  56. {"r horizon" -> N@rA, "r ergosphere" -> N@rE, "r isco" -> N@isco, "r photon pro" -> N@Min[rp], "r photon ret" -> N@Max[rp], "r disk" -> N@sr, "r observer" -> N@r0, "θ observer" -> N@θ0}
  57.  
  58. j[v_] := Sqrt[1-v^2];
  59. Ы[rs_] := Sqrt[Χs[rs]/Σs[rs]];
  60. ωs[rs_] := (a (2 rs - ℧^2))/Χs[rs];
  61.  
  62. ε[rs_] := Sqrt[Δs[rs] Σs[rs]/Χs[rs]]/j[vt]+Lz[rs] ωs[rs];
  63. Lz[rs_] := vt Ы[rs]/j[vt];
  64.  
  65. red[rs_] := Quiet[Reduce[
  66. dt == (Lz[rs] (-a (a^2+rs^2)+Δs[rs] κs[rs])+ε[rs] ((a^2+rs^2)^2-Δs[rs] κs[rs]^2))/(Δs[rs] Σs[rs]) &&
  67. 0 == ((a^2+(-2+rs) rs+℧^2) (16 a dt dΦ rs (rs-℧^2)+8 dt^2 rs (-rs+℧^2)+dΦ^2 rs (8 rs (-a^2+rs^3)+a^2 (4 a^2+4 ℧^2-4 (a-℧) (a+℧)))))/(8 rs^6) &&
  68. dΦ == (Lz[rs] (-a^2+Δs[rs])+ε[rs] (a (a^2+rs^2)-Δs[rs] κs[rs]))/(Δs[rs] Σs[rs]) &&
  69. vt > 0,
  70. {vt,dΦ,dt},
  71. Reals]];
  72.  
  73. vs = Interpolation[Table[{rr, If[Quiet@NumericQ[red[rr][[1, 2]]], red[rr][[1, 2]], 0]}, {rr, 0, sr, 0.02}]];
  74. φs = Interpolation[Table[{rr, If[Quiet@NumericQ[red[rr][[2, 2]]], red[rr][[2, 2]], 0]}, {rr, 0, sr, 0.02}]];
  75. ts = Interpolation[Table[{rr, If[Quiet@NumericQ[red[rr][[3, 2]]], red[rr][[3, 2]], 0]}, {rr, 0, sr, 0.02}]];
  76.  
  77. plot[func_, label_] := Plot[func, {rr, rP, sr}, GridLines -> {{Min[rp], Max[rp], rA, si, isco, rE, sr}, {}},
  78. Frame -> True, ImagePadding -> {{40,1}, {12,1}}, ImageSize -> 340, PlotLabel -> label, PlotRange->{{0,sr}, Automatic}]
  79.  
  80. plot[Sqrt[Χs[rr]/Δs[rr]/Σs[rr]], "Gravitational time dilation: y=dt/dт, x=r"]
  81. plot[ts[rr], "Total time dilation: y=dt/dτ, x=r"]
  82. plot[(a (2 rr-℧^2))/((a^2+rr^2)^2-a^2 (a^2-2 rr+rr^2+℧^2)), "Frame Dragging: y=dφ/dт, x=r"]
  83. plot[φs[rr]/ts[rr], "Shapirodelayed angular velocity: y=dφ/dt, x=r"]
  84. plot[φs[rr], "Coordinate speed: y=dφ/dτ, x=r"]
  85. plot[vs[rr], "Local velocity: y=v=dl/dτ, x=r"]
  86.  
  87. gtt = (2r0-℧^2)/Σ-1;
  88. grr = Σ/Δ;
  89. gθθ = Σ;
  90. gφφ = Χ/Σ Sin[θ0]^2;
  91. gtφ =-a (2r0-℧^2) Sin[θ0]^2/Σ;
  92.  
  93. й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-
  94. 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-
  95. 2 r0+r0^2+℧^2) Sin[θ0]^2)/(r0^2+a^2 Cos[θ0]^2)]);
  96.  
  97. U={+vr, +vθ, +vφ};
  98. γ=1/Sqrt[1-Norm[U]^2];
  99.  
  100. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z};
  101. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  102. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  103.  
  104. raytracer[{Ф_, ϑ_}] :=
  105.  
  106. Quiet[Module[{V, W, vw, vr0i, vθ0i, vφ0i, vr0n, vθ0n, vφ0n, vr0a, vθia, vφ0a, vt0a,
  107. DGL, sol, εj, pθi, pr0, Q, k, t10, r10, Θ10, Φ10, т0, т1, R1, t, r, θ, φ, τ, plunge, plunge2,
  108. X, Y, Z, tt, rt, θt, φt, т, ξ, stepsize, laststep, mtl, mta, ft, fτ, varb, Σj, Δj, Χj,
  109. ωj, dθ, dφ, dphi, dtheta, drad, dtau, vrj, vθj, vφj, vtj, vφt, shf, rnd},
  110.  
  111. vw=xyZ[Xyz[{0, 1, 0}, ϑ], Ф+π/2];
  112. (* Übersetzung des Einfallswinkels in den lokalen Tetrad *)
  113. vr0a = vw[[3]] Sqrt[grr];
  114. vφ0a = vw[[2]] Sqrt[gφφ]/r0/Sin[θi];
  115. vθia = vw[[1]] Sqrt[gθθ]/r0;
  116. (* Betrag *)
  117. vt0a = Sqrt[vr0a^2+vφ0a^2+vθia^2];
  118. (* Normierung *)
  119. vr0n = vr0a/vt0a;
  120. vφ0n = vφ0a/vt0a;
  121. vθ0n = vθia/vt0a;
  122. (* Relativistische Geschwindigkeitsaddition *)
  123. V={vr0n, vθ0n, vφ0n};
  124. W=(U+V+γ/(1+γ)(U\[Cross](U\[Cross]V)))/(1+U.V);
  125. (* Aberration *)
  126. vr0i = W[[1]];
  127. vθ0i = W[[2]];
  128. vφ0i = W[[3]];
  129.  
  130. rnd = If[round == 2,
  131. {WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0,
  132. (plunge=τ) && (tt=If[snc==1, 0, t[τ]]) && (rt=r[τ]) && (θt=θ[τ]) &&
  133. (φt=φ[τ]) && (dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) &&
  134. (dphi=φ'[τ]); "StopIntegration"]},
  135. If[round == 1,
  136. {WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0,
  137. (plunge=τ) && (tt=If[snc==1, 0, t[τ]]) && (rt=r[τ]) && (θt=θ[τ]) &&
  138. (φt=φ[τ]) && (dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) &&
  139. (dphi=φ'[τ]); "StopIntegration"]},
  140. If[round == 4,
  141. {WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0,
  142. (plunge2=τ)],
  143. WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]<0
  144. &&τ<plunge2-0.1,(plunge=τ)&&(tt=If[snc==1,0,t[τ]])&&(rt=r[τ])&&(θt=θ[τ])&&(φt=φ[τ]) &&
  145. (dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) && (dphi=φ'[τ]); "StopIntegration"]},
  146. If[round == 3,
  147. {WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0,
  148. (plunge2=τ)],
  149. WhenEvent[Mod[θ[τ],π]==π/2.0 && r[τ]>si && r[τ]<sr && θ'[τ]>0
  150. &&τ<plunge2-0.1,(plunge=τ)&&(tt=If[snc==1,0,t[τ]])&&(rt=r[τ])&&(θt=θ[τ])&&(φt=φ[τ]) &&
  151. (dtau=t'[τ]) && (drad=r'[τ]) && (dtheta=θ'[τ]) && (dphi=φ'[τ]); "StopIntegration"]},
  152. {WhenEvent[Mod[θ[τ],π]==Pi/2.0 && r[τ]>si && r[τ]<sr, (plunge=τ) &&
  153. (tt=If[snc==1, 0, t[τ]]) && (rt=r[τ]) && (θt=θ[τ]) && (φt=φ[τ]) && (dtau=t'[τ]) &&
  154. (drad=r'[τ]) && (dtheta=θ'[τ]) && (dphi=φ'[τ]); "StopIntegration"]}
  155. ]]]];
  156.  
  157. DGL = { (* Kerr Newman Bewegungsgleichungen *)
  158.  
  159. t''[τ]==-(((r'[τ] ((a^2+r[τ]^2) (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+
  160. 2 (-℧^2+r[τ]) t'[τ]))+a (2 a^4 Cos[θ[τ]]^2+a^2 ℧^2 (3+Cos[2 θ[τ]]) r[τ]-
  161. a^2 (3+Cos[2 θ[τ]]) r[τ]^2+4 ℧^2 r[τ]^3-6 r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))/(a^2+℧^2+(-2+
  162. r[τ]) r[τ])+a^2 θ'[τ] (Sin[2 θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])-2 a Cos[θ[τ]] (℧^2-
  163. 2 r[τ]) Sin[θ[τ]]^3 φ'[τ]))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^2),
  164.  
  165. t'[0]==ς,
  166. t[0]==0,
  167.  
  168. r''[τ]==((-1+r[τ])/(a^2+℧^2+(-2+r[τ]) r[τ])-r[τ]/(a^2 Cos[θ[τ]]^2+r[τ]^2)) r'[τ]^2+
  169. (a^2 Sin[2 θ[τ]] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(8 (a^2 Cos[θ[τ]]^2+
  170. r[τ]^2)^3))(a^2+℧^2+(-2+r[τ]) r[τ]) (8 t'[τ] (a^2 Cos[θ[τ]]^2 (-q ℧+t'[τ])+
  171. r[τ] (q ℧ r[τ]+(℧^2-r[τ]) t'[τ]))+8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+
  172. 8 a Sin[θ[τ]]^2 (a^2 Cos[θ[τ]]^2 (q ℧-2 t'[τ])+r[τ] (-q ℧ r[τ]+2 (-℧^2+r[τ]) t'[τ])) φ'[τ]+
  173. Sin[θ[τ]]^2 (r[τ] (a^2 (3 a^2+4 ℧^2+4 (a-℧) (a+℧) Cos[2 θ[τ]]+a^2 Cos[4 θ[τ]])+
  174. 8 r[τ] (2 a^2 Cos[θ[τ]]^2 r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]^2),
  175.  
  176. r'[0]==vr0i/Sqrt[(r0^2+a^2 Cos[θi]^2)/(a^2+(-2+r0) r0+℧^2)],
  177. r[0]==r0,
  178.  
  179. θ''[τ]==-((a^2 Cos[θ[τ]] Sin[θ[τ]] r'[τ]^2)/((a^2+℧^2+(-2+r[τ]) r[τ]) (a^2 Cos[θ[τ]]^2+
  180. r[τ]^2)))-(2 r[τ] r'[τ] θ'[τ])/(a^2 Cos[θ[τ]]^2+r[τ]^2)+(1/(16 (a^2 Cos[θ[τ]]^2+
  181. r[τ]^2)^3))Sin[2 θ[τ]] (a^2 (-8 t'[τ] (2 q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+8 (a^2 Cos[θ[τ]]^2+
  182. r[τ]^2)^2 θ'[τ]^2)+16 a (a^2+r[τ]^2) (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ]) φ'[τ]+(3 a^6-5 a^4 ℧^2+
  183. 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+
  184. a^4 Cos[4 θ[τ]] (a^2+℧^2+(-2+r[τ]) r[τ])+4 a^2 Cos[2 θ[τ]] (a^2+℧^2+(-2+
  185. r[τ]) r[τ]) (a^2+2 r[τ]^2)) φ'[τ]^2),
  186.  
  187. θ'[0]==vθ0i/Sqrt[r0^2+a^2 Cos[θi]^2],
  188. θ[0]==θi,
  189.  
  190. φ''[τ]==-(1/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2))((r'[τ] (4 a q ℧ (a^2 Cos[θ[τ]]^2-r[τ]^2)-
  191. 8 a (a^2 Cos[θ[τ]]^2+(℧^2-r[τ]) r[τ]) t'[τ]+(a^2 (3 a^2+8 ℧^2+a^2 (4 Cos[2 θ[τ]]+
  192. Cos[4 θ[τ]])) r[τ]-4 a^2 (3+Cos[2 θ[τ]]) r[τ]^2+8 (a^2+℧^2+a^2 Cos[2 θ[τ]]) r[τ]^3-
  193. 16 r[τ]^4+8 r[τ]^5+2 a^4 Sin[2 θ[τ]]^2) φ'[τ]))/(a^2+℧^2+(-2+r[τ]) r[τ])+
  194. θ'[τ] (8 a Cot[θ[τ]] (q ℧ r[τ]+(℧^2-2 r[τ]) t'[τ])+(8 Cot[θ[τ]] (a^2+r[τ]^2)^2-
  195. 2 a^2 (3 a^2+2 ℧^2+4 (-1+r[τ]) r[τ]) Sin[2 θ[τ]]-a^4 Sin[4 θ[τ]]) φ'[τ])),
  196.  
  197. φ'[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+
  198. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)]+a (2 r0-℧^2) (Sqrt[((a^2+(-2+r0) r0+℧^2) (r0^2+
  199. 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-
  200. ℧^2) Sin[θi])/((r0^2+a^2 Cos[θi]^2) Sqrt[((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0+
  201. ℧^2) Sin[θi]^2)/(r0^2+a^2 Cos[θi]^2)])))/((a^2+(-2+r0) r0+℧^2) (r0^2+a^2 Cos[θi]^2)),
  202. φ[0]==φ0,
  203.  
  204. rnd};
  205. (* Integrator *)
  206. sol = NDSolve[DGL, {t, r, θ, φ}, {τ, 0, tmax},
  207. WorkingPrecision-> wp,
  208. MaxSteps-> Infinity,
  209. InterpolationOrder-> All];
  210.  
  211. Σj = rt^2;
  212. Δj = rt^2-2 rt+a^2+℧^2;
  213. Χj = (rt^2+a^2)^2-a^2 Δj;
  214. ωj = (a(2 rt-℧^2))/Χj;
  215.  
  216. т0 = If[NumericQ[plunge], plunge, tmax];
  217. т1 = If[NumericQ[plunge], tt, 0];
  218.  
  219. εj = dtau (1-(2 rt-℧^2)/rt^2)+(a dphi (2 rt-℧^2))/rt^2;
  220.  
  221. vrj = drad/Sqrt[Δj] Sqrt[Σj];
  222. vθj = dtheta Sqrt[Σj];
  223. vφj = (-(((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+rt^2) Sin[(π/2)] (-dphi-
  224. (a q ℧ rt)/((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+rt^2))+
  225. (εj Csc[(π/2)]^2 (a (-a^2-℧^2+2 rt-rt^2) Sin[(π/2)]^2+a (a^2+
  226. rt^2) Sin[(π/2)]^2))/((a^2 Cos[(π/2)]^2+rt^2) (a^2+℧^2-2 rt+rt^2))+(a q ℧ rt (a^2+
  227. ℧^2-2 rt+rt^2-a^2 Sin[(π/2)]^2))/((a^2 Cos[(π/2)]^2+rt^2)^2 (a^2+℧^2-2 rt+
  228. rt^2))))/((a^2+℧^2-2 rt+rt^2-a^2 Sin[(π/2)]^2) Sqrt[((a^2+rt^2)^2-
  229. a^2 (a^2+℧^2-2 rt+rt^2) Sin[(π/2)]^2)/(a^2 Cos[(π/2)]^2+rt^2)])));
  230. vtj = Sqrt[vrj^2+vθj^2+vφj^2];
  231. vφt = vφj/vtj;
  232. shf = ς/(1-vs[rt] vφt)/dtau;
  233.  
  234. dφ = If[snc==1, 0, If[snc==2, т1 φs[rt]/ts[rt], If[snc==3, (т1+r0) φs[rt]/ts[rt], (т1-r0) φs[rt]/ts[rt]]]];
  235.  
  236. If[т0>tmax+1,If[rt>sr,{0,0},If[rt<si,{0,0},{0,Min[2, shf]}]],{0,0}]]];
  237.  
  238. width=ImageDimensions[pic][[1]];
  239. height=ImageDimensions[pic][[2]];
  240. hzoom=If[breite>2 hoehe,1/zoom,1/zoom/2/hoehe*breite];
  241. vzoom=If[breite>2 hoehe,1/zoom*2 hoehe/breite,1/zoom];
  242.  
  243. FOV->{360.0 hzoom "degree",180.0 vzoom "degree"}
  244.  
  245. img=ParallelTable[ImageTransformation[pic,raytracer,{breite,Ceiling[hoehe/kernels/grain]},
  246. DataRange->{{0,2π},{0,2}},
  247. PlotRange->{{-π+hvs/hzoom, π+hvs/hzoom} hzoom, {-π/2+vvs/vzoom+x, -π/2+vvs/vzoom+x+π/kernels/grain} vzoom},
  248. Padding->"Fixed"],
  249. {x, 0, π-π/kernels/grain, π/kernels/grain}]
  250.  
  251. image = ImageAssemble[Table[{img[[x]]}, {x, kernels grain, 1, -1}]]
RAW Paste Data