API tools faq
paste
Yukterez

Raytracer (Accretion Disk)

Yukterez
Jun 8th, 2019 (edited)
354
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 13.37 KB
raw download clone embed print report
  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
Public Pastes
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
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!