Yukterez

Kerr Photon Orbit

Jul 21st, 2017
154
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* |||||||| Mathematica Syntax | http://kerr.yukterez.net | Version: 20.07.2017 |||||||| *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. mt1={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
  6. mt2={"EventLocator", "Event"-> (r[t]-1001/1000 rA)};
  7. mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
  8. mt4={"EquationSimplification"-> "Residual"};
  9. mt0=Automatic;
  10. mta=mt0;
  11. wp=MachinePrecision;
  12.  
  13. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  14. (* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *)
  15. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  16.  
  17. A=a; (* pseudosphärisch [BL]: A=0, kartesisch [KS]: A=a *)
  18.  
  19. tmax=500; (* Eigenzeit *)
  20. Tmax=400; (* Koordinatenzeit *)
  21.  
  22. r0=3; (* Startradius *)
  23. θ0=π/2; (* Breitengrad *)
  24. φ0=0; (* Längengrad *)
  25. q=0; (* elektrische Ladung *)
  26. a=1; (* Spinparameter *)
  27. μ=0; (* Baryon: μ=-1, Photon: μ=0 *)
  28.  
  29. v0=1; (* Anfangsgeschwindigkeit *)
  30. α0=0; (* vertikaler Abschusswinkel *)
  31. ψ0=π/2+ArcSin[1/3]; (* Bahninklinationswinkel *)
  32.  
  33. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  34. (* |||||||| 2) GESCHWINDIGKEITS-, ENERGIE UND DREHIMPULSKOMPONENTEN ||||||||||||||||||||| *)
  35. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  36.  
  37. vr0=v0 Sin[α0]; (* radiale Geschwindigkeitskomponente *)
  38. vθ0=v0 Cos[α0] Sin[ψ0]; (* longitudinale Geschwindigkeitskomponente *)
  39. vφ0=v0 Cos[α0] Cos[ψ0]; (* latitudinale Geschwindigkeitskomponente *)
  40.  
  41. x0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0]; (* Anfangskoordinaten *)
  42. y0[A_]:=Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
  43. z0[A_]:=r0 Cos[θ0];
  44.  
  45. ε=Sqrt[δ Ξ/χ]/j[v0]+Lz щ; (* Energie und Drehimpulskomponenten *)
  46. Lz=vφ0 Ы/j[v0];
  47. pθ0=vθ0 Sqrt[Ξ]/j[v0];
  48. pr0=vr0 Sqrt[(Ξ/δ)/j[v0]^2];
  49. Q=pθ0^2+(Lz^2 Csc[θ0]^2-a^2 (ε^2+μ)) Cos[θ0]^2; (* Carter Konstante *)
  50. k=Q+Lz^2+a^2 (ε^2+μ); (* Carter k *)
  51.  
  52. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  53. (* |||||||| 3) RADIUS NACH GESCHWINDIGKEIT UND VICE VERSA ||||||||||||||||||||||||||||||| *)
  54. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  55.  
  56. rPro=2 (1+Cos[2/3 ArcCos[-a]]); (* prograder Photonenorbitradius *)
  57. rRet=2 (1+Cos[2/3 ArcCos[+a]]); (* retrograder Photonenorbitradius *)
  58. rTeo=1+2 Sqrt[1-a^3/3] Cos[ArcCos[(1-a^2)/(1-a^2/3)^(3/2)]/3];
  59.  
  60. vPro=(a^2-2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a+r0^(3/2))); (* Kreisgeschwindigkeit + *)
  61. vRet=(a^2+2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a-r0^(3/2))); (* Kreisgeschwindigkeit - *)
  62.  
  63. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  64. (* |||||||| 4) HORIZONTE UND ERGOSPHÄREN RADIEN ||||||||||||||||||||||||||||||||||||||||| *)
  65. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  66.  
  67. rE=1+Sqrt[1-a^2 Cos[θ]^2-q^2]; (* äußere Ergosphäre *)
  68. RE[A_, w1_, w2_]:=Xyz[xyZ[
  69. {Sqrt[rE^2+A^2] Sin[θ] Cos[φ], Sqrt[rE^2+A^2] Sin[θ] Sin[φ], rE Cos[θ]}, w1], w2];
  70. rG=1-Sqrt[1-a^2 Cos[θ]^2-q^2]; (* innere Ergosphäre *)
  71. RG[A_, w1_, w2_]:=Xyz[xyZ[
  72. {Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2];
  73. rA=1+Sqrt[1-a^2-q^2]; (* äußerer Horizont *)
  74. RA[A_, w1_, w2_]:=Xyz[xyZ[
  75. {Sqrt[rA^2+A^2] Sin[θ] Cos[φ], Sqrt[rA^2+A^2] Sin[θ] Sin[φ], rA Cos[θ]}, w1], w2];
  76. rI=1-Sqrt[1-a^2-q^2]; (* innerer Horizont *)
  77. RI[A_, w1_, w2_]:=Xyz[xyZ[
  78. {Sqrt[rI^2+A^2] Sin[θ] Cos[φ], Sqrt[rI^2+A^2] Sin[θ] Sin[φ], rI Cos[θ]}, w1], w2];
  79.  
  80. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  81. (* |||||||| 5) HORIZONTE UND ERGOSPHÄREN PLOT ||||||||||||||||||||||||||||||||||||||||||| *)
  82. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  83.  
  84. horizons[A_, mesh_, w1_, w2_]:=Show[
  85. ParametricPlot3D[RE[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  86. Mesh -> mesh, PlotPoints -> plp, PlotStyle -> Directive[Blue, Opacity[0.10]]],
  87. ParametricPlot3D[RA[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  88. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
  89. ParametricPlot3D[RI[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  90. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.25]]],
  91. ParametricPlot3D[RG[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  92. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.35]]]];
  93. BLKS:=Grid[{{horizons[a, 35, 0, 0], horizons[0, 35, 0, 0]}}];
  94.  
  95. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  96. (* |||||||| 6) FUNKTIONEN ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  97. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  98.  
  99. j[v_]:=Sqrt[1+μ v^2]; (* Lorentzfaktor *)
  100. щ=2r0 a/χ; (* Frame Drag *)
  101. я=Sqrt[Χ/Σ]Sin[θ[τ]]; (* axialer Umfangsradius *)
  102. яi[τ_]:=Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
  103. Ы=Sqrt[χ/Ξ]Sin[θ0];
  104. Σ=r[τ]^2+a^2 Cos[θ[τ]]^2; (* poloidialer Umfangsradius *)
  105. Σi[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2;
  106. Ξ=r0^2+a^2 Cos[θ0]^2;
  107. Δ=r[τ]^2-2r[τ]+a^2+q^2;
  108. Δi[τ_]:=R[τ]^2-2R[τ]+a^2+q^2;
  109. δ=r0^2-2r0+a^2+q^2;
  110. Щ=Lz^2 Cot[θ[τ]]^2;
  111. Χ=(r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 Δ;
  112. Χi[τ_]:=(R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 Δi[τ];
  113. χ=(r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
  114.  
  115. т[τ_]:=Evaluate[t[τ]/.sol][[1]]; (* Koordinatenzeit nach Eigenzeit *)
  116. д[ξ_] :=Quiet[Ξ /.FindRoot[т[Ξ]-ξ, {Ξ, 0}]]; (* Eigenzeit nach Koordinatenzeit *)
  117. T :=Quiet[д[tk]];
  118.  
  119. ю[τ_]:=Evaluate[t'[τ]/.sol][[1]];
  120. γ[τ_]:=If[μ==0, "Infinity", ю[τ]]; (* totale ZD *)
  121. R[τ_]:=Evaluate[r[τ]/.sol][[1]]; (* Boyer-Lindquist Radius *)
  122. Φ[τ_]:=Evaluate[φ[τ]/.sol][[1]];
  123. Θ[τ_]:=Evaluate[θ[τ]/.sol][[1]];
  124. ß[τ_]:=Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ];
  125.  
  126. ς[τ_]:=Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0=Sqrt[χ/δ/Ξ]; (* gravitative ZD *)
  127. Λ[τ_]:=R[τ]^2+a^2-2 R[τ]; Λ0=r0^2+a^2-2 r0;
  128. Υ[τ_]:=(R[τ]^2+a^2)^2-a^2 Λ[τ] Sin[Θ[τ]]^2; Υ0=(r0^2+a^2)^2-a^2 Λ0 Sin[θ0]^2;
  129. ρ[τ_]:=R[τ]^2+a^2 Cos[Θ[τ]]^2; ρ0=r0^2+a^2 Cos[θ0]^2;
  130. ω[τ_]:=2R[τ] a/Υ[τ]; ω0=2r0 a/Υ0; (* Frame Dragging Winkelgeschwindigkeit *)
  131. Ω[τ_]:=ω[τ] Sqrt[X[τ]^2+Y[τ]^2]; (* Frame Dragging beobachtete Geschwindigkeit *)
  132. й[τ_]:=ω[τ] яi[τ] ς[τ]; й0=ω0 Ы ς0; (* Frame Dragging lokale Geschwindigkeit *)
  133.  
  134. ж[τ_]:=Sqrt[ς[τ]^2-1]/ς[τ]; ж0=Sqrt[ς0^2-1]/ς0; (* Fluchtgeschwindigkeit *)
  135. v[τ_]:=If[μ==0, 1, Abs[Re[-((\[Sqrt](-a^4(ε-Lz ω[τ])^2-2 a^2R[τ]^2 (ε-Lz ω[τ])^2-
  136. R[τ]^4(ε-Lz ω[τ])^2+Δi[τ](Σi[τ]+a^2 Sin[Θ[τ]]^2 (ε-
  137. Lz ω[τ])^2)))/(Sqrt[-(a^2+R[τ]^2)^2+
  138. a^2 Sin[Θ[τ]]^2 Δi[τ]](ε - Lz ω[τ])))]]]; (* lokale Dreiergeschwindigkeit *)
  139. pΘ[τ_]:=Evaluate[pθ[τ] /. sol][[1]];
  140. pR[τ_]:=Evaluate[pr[τ] /. sol][[1]];
  141. sh[τ_]:=Sqrt[ß[τ]^2-Ω[τ]^2];
  142. epot[τ_]:=ε-1-ekin[τ]; (* potentielle Energie *)
  143. ekin[τ_]:=If[μ==0, ς[τ], 1/Sqrt[1-v[τ]^2]-1]; (* kinetische Energie *)
  144.  
  145. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  146. (* |||||||| 7) DIFFERENTIALGLEICHUNG |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  147. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  148.  
  149. dp= \!\(\*SuperscriptBox[\(Y\),\(Y\)]\); n0[z_]:=Chop[N[z]];
  150.  
  151. DGL={
  152. t'[τ]==ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ),
  153. t[0]==0,
  154. r'[τ]==(pr[τ] Δ)/Σ,
  155. r[0]==r0,
  156. θ'[τ]==pθ[τ]/Σ,
  157. θ[0]==θ0,
  158. φ'[τ]==(2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ),
  159. φ[0]==φ0,
  160. pr'[τ]==1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+
  161. 2r[τ](r[τ]^2+a^2) ε^2-2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ,
  162. pr[0]==pr0,
  163. pθ'[τ]==(Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ)),
  164. pθ[0]==pθ0
  165. };
  166.  
  167. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  168. (* |||||||| 8) INTEGRATION ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  169. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  170.  
  171. sol=NDSolve[DGL, {t, r, θ, φ, pr, pθ}, {τ, 0, tmax},
  172. WorkingPrecision-> wp,
  173. MaxSteps-> Infinity,
  174. Method-> mta,
  175. InterpolationOrder-> All];
  176.  
  177. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  178. (* |||||||| 9) KOORDINATEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  179. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  180.  
  181. X[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]]; (* kartesisch *)
  182. Y[τ_]:=Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
  183. Z[τ_]:=Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];
  184.  
  185. x[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]]; (* Plotkoordinaten *)
  186. y[τ_]:=Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
  187. z[τ_]:=Z[τ];
  188.  
  189. XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2]; (* kartesischer Radius *)
  190.  
  191. Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z}; (* Rotationsmatrix *)
  192. xYz[{x_, y_, z_}, β_]:={x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  193. xyZ[{x_, y_, z_}, ψ_]:={x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  194.  
  195. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  196. (* |||||||| 10) PLOT EINSTELLUNGEN |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  197. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  198.  
  199. PR=1.2r0; (* Plot Range *)
  200. VP={r0, r0, r0}; (* Perspektive x,y,z*)
  201. d1=10; (* Schweiflänge *)
  202. plp=50; (* Flächenplot Details *)
  203. w1l=0; w2l=0; w1r=0; w2r=0; (* Startperspektiven *)
  204. Mrec=100; mrec=10; (* Parametric Plot Subdivisionen *)
  205. imgsize=380; (* Bildgröße *)
  206.  
  207. s[text_]:=Style[text, FontSize->font]; font=11; (* Anzeigestil *)
  208.  
  209. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  210. (* |||||||| 11) PLOT NACH KOORDINATENZEIT ||||||||||||||||||||||||||||||||||||||||||||||| *)
  211. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  212.  
  213. display[T_]:=Grid[{
  214. {s[" t coord"], " = ", s[n0[tk]], s["GM/c³"], s[dp]},
  215. {If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c³"], s[dp]},
  216. {s[" γ total"], " = ", s[n0[γ[T]]], s["dt/dτ"], s[dp]},
  217. {s[" ς gravt"], " = ", s[n0[ς[T]]], s["dt/dτ"], s[dp]},
  218. {s[" r coord"], " = ", s[n0[R[T]]], s["GM/c²"], s[dp]},
  219. {s[" φ longd"], " = ", s[n0[Φ[T]]], s["rad"], s[dp]},
  220. {s[" θ lattd"], " = ", s[n0[Θ[T]]], s["rad"], s[dp]},
  221. {s[" Ř crc.φ"], " = ", s[n0[яi[T]]], s["GM/c²"], s[dp]},
  222. {s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[T]]]], s["GM/c²"], s[dp]},
  223. {s[" Δ rad.f"], " = ", s[n0[Sqrt[Δi[T]]]], s["GM/c²"], s[dp]},
  224. {s[" E kinet"], " = ", s[n0[ekin[T]]], s["mc²"], s[dp]},
  225. {s[" E poten"], " = ", s[n0[epot[T]]], s["mc²"], s[dp]},
  226. {s[" E total"], " = ", s[n0[ε]], s["mc²"], s[dp]},
  227. {s[" CarterQ"], " = ", s[N[Q]], s["GMm/c"], s[dp]},
  228. {s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
  229. {s[" L polar"], " = ", s[n0[pΘ[T]]], s["GMm/c"], s[dp]},
  230. {s[" p r.mom"], " = ", s[n0[pR[T]]], s["mc"], s[dp]},
  231. {s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c²"], s[dp]},
  232. {s[" x carts"], " = ", s[n0[X[T]]], s["GM/c²"], s[dp]},
  233. {s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c²"], s[dp]},
  234. {s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c²"], s[dp]},
  235. {s[" ω fdrag"], " = ", s[n0[ω[T]]], s["c³/G/M"], s[dp]},
  236. {s[" v fdrag"], " = ", s[n0[й[T]]], s["c"], s[dp]},
  237. {s[" Ω fdrag"], " = ", s[n0[Ω[T]]], s["c"], s[dp]},
  238. {s[" v obsvd"], " = ", s[n0[ß[T]]], s["c"], s[dp]},
  239. {s[" v escpe"], " = ", s[n0[ж[T]]], s["c"], s[dp]},
  240. {s[" v delay"], " = ", s[n0[sh[T]]], s["c"], s[dp]},
  241. {s[" v local"], " = ", s[n0[v[T]]], s["c"], s[dp]},
  242. {s[" "], s[" "], s[" "], s[" "]}},
  243. Alignment-> Left, Spacings-> {0, 0}];
  244.  
  245. plot1a[{xx_, yy_, zz_, tk_, w1_, w2_}]:= (* Animation *)
  246. Rasterize[
  247. Show[Graphics3D[{
  248. {PointSize[0.009], Red, Point[
  249. Xyz[xyZ[{x[T], y[T], z[T]}, w1], w2]]}},
  250. ImageSize-> imgsize,
  251. PlotRange-> PR,
  252. SphericalRegion->False,
  253. ImagePadding-> 1],
  254. horizons[A, None, w1, w2],
  255. If[a==0, {},
  256. Graphics3D[{{PointSize[0.009], Purple, Point[
  257. Xyz[xyZ[{
  258. Sin[-φ0-щ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
  259. Cos[-φ0-щ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
  260. z0[A]}, w1], w2]]}}]],
  261. If[tk==0, {}, If[a==0, {},
  262. ParametricPlot3D[
  263. Xyz[xyZ[{
  264. Sin[-φ0-щ tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
  265. Cos[-φ0-щ tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
  266. z0[A]}, w1], w2],
  267. {tt, Max[0, tk-199/100 π/щ], tk},
  268. PlotStyle -> {Thickness[0.001], Dashed, Purple},
  269. PlotPoints-> Automatic,
  270. MaxRecursion-> mrec]]],
  271. If[tk==0, {},
  272. Block[{$RecursionLimit = Mrec},
  273. ParametricPlot3D[
  274. Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, 0, Max[1*^-16, T-d1/3]},
  275. PlotStyle-> {Thickness[0.003], Gray},
  276. PlotPoints-> Automatic,
  277. MaxRecursion-> mrec]]],
  278. Block[{$RecursionLimit = Mrec},
  279. If[tk==0, {},
  280. ParametricPlot3D[
  281. Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, Max[0, T-d1], T},
  282. PlotStyle-> {Thickness[0.004]},
  283. ColorFunction-> Function[{x, y, z, t},
  284. Hue[0, 1, 0.5, Max[Min[(-T+(t+d1))/d1, 1], 0]]],
  285. ColorFunctionScaling-> False,
  286. PlotPoints-> Automatic,
  287. MaxRecursion-> mrec]]],
  288. ViewPoint-> {xx, yy, zz}]];
  289.  
  290. Do[
  291. Print[Rasterize[Grid[{{
  292. plot1a[{0, -Infinity, 0, tk, w1l, w2l}],
  293. plot1a[{0, 0, Infinity, tk, w1r, w2r}],
  294. display[Quiet[д[tk]]]
  295. }}, Alignment->Left]]],
  296. {tk, 0, Tmax, Tmax/10}]
  297.  
  298. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  299. (* |||||||||||||||| http://kerr.yukerez.net ||| Simon Tyran, Vienna ||||||||||||||||| *)
  300. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
RAW Paste Data