Yukterez

Raytracer (Accretion Disk)

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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×