Public Pastes
SHARE
TWEET

Raytracer (Accretion Disk)

Yukterez Jun 8th, 2019 (edited) 273 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
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