Public Pastes
SHARE
TWEET

Raytracer (Accretion Disk)

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