Yukterez

Pure Kerr 3D Simulator

Aug 21st, 2017
207
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* |||| Mathematica Syntax | kerr.yukterez.net | 22.06.2016 - 04.04.2019, Version 17 |||| *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
  6. mt2 = {"EventLocator", "Event"-> (r[τ]-1001/1000 rA)};
  7. mt3 = {"ImplicitRungeKutta", "DifferenceOrder"-> 20};
  8. mt4 = {"EquationSimplification"-> "Residual"};
  9. mt0 = Automatic;
  10. mta = mt0; (* mt0: Speed, mt3: Accuracy *)
  11. wp = MachinePrecision;
  12.  
  13. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  14. (* |||||||| 1) STARTBEDINGUNGEN EINGEBEN |||||||||||||||||||||||||||||||||||||||||||||||| *)
  15. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  16.  
  17. A = a; (* pseudosphärisch: A=0, kartesisch: A=a *)
  18. crd = 1; (* Boyer-Lindquist: crd=1, Kerr-Schild: crd=2 *)
  19. dsp = 1; (* Display Modus *)
  20.  
  21. tmax = 120; (* Eigenzeit *)
  22. Tmax = 120; (* Koordinatenzeit *)
  23. TMax = Min[Tmax, т[plunge-1*^-3]]; tMax = Min[tmax, plunge]; (* Integrationsende *)
  24.  
  25. r0 = 3; (* Startradius *)
  26. r1 = 5; (* Endradius wenn v0=vr0=vr1 *)
  27. θ0 = π/2; (* Breitengrad *)
  28. φ0 = 0; (* Längengrad *)
  29. a = 1; (* Spinparameter *)
  30.  
  31. v0 = 1; (* Anfangsgeschwindigkeit *)
  32. α0 = 0; (* vertikaler Abschusswinkel *)
  33. ψ0 = δp[r0, a]; (* Bahninklinationswinkel *)
  34.  
  35. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  36. (* |||||||| 2) GESCHWINDIGKEITS-, ENERGIE UND DREHIMPULSKOMPONENTEN ||||||||||||||||||||| *)
  37. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  38.  
  39. vr0 = v0 Sin[α0]; (* radiale Geschwindigkeitskomponente *)
  40. vθ0 = v0 Cos[α0] Sin[ψ0]; (* longitudinale Geschwindigkeitskomponente *)
  41. vφ0 = v0 Cos[α0] Cos[ψ0]; (* latitudinale Geschwindigkeitskomponente *)
  42.  
  43. vrj[τ_]:=R'[τ]/Sqrt[Δi[τ]] Sqrt[Σi[τ] (1+μ v[τ]^2)];
  44. vθj[τ_]:=Θ'[τ] Sqrt[Σi[τ] (1+μ v[τ]^2)];
  45. vφBL[τ_]:=-((Sin[Θ[τ]] Sqrt[1+μ v[τ]^2] (-a^5 ε Cos[Θ[τ]]^2-a ε R[τ]^4+
  46. a^2 Δi[τ] (xJ[τ] ε Cot[Θ[τ]]^2+Φ'[τ] Cos[Θ[τ]]^2 Σi[τ])+
  47. R[τ]^2 (-a^3 ε (1+Cos[Θ[τ]]^2)+Δi[τ] (xJ[τ] ε Csc[Θ[τ]]^2+Φ'[τ] Σi[τ]))))/((a^2 Cos[Θ[τ]]^2+
  48. R[τ]^2) (a^2 Sin[Θ[τ]]^2-Δi[τ]) Sqrt[Χi[τ]/Σi[τ]]));
  49. vφKS[τ_]:=(j[v[τ]] Sin[Θ[τ]]^2 (2 a ε R[τ]-Φ'[τ] Δi[τ] Σi[τ]+
  50. a Σi[τ] R'[τ]))/(Ыi[τ] (2 R[τ]-Σi[τ]))
  51. vφj[τ_]:=If[crd==2, vφKS[τ], vφBL[τ]];
  52.  
  53. vtj[τ_]:=Sqrt[vrj[τ]^2+vθj[τ]^2+vφj[τ]^2];
  54. vr[τ_]:=vrj[τ]/vtj[τ]*v[τ];
  55. vθ[τ_]:=vθj[τ]/vtj[τ]*v[τ];
  56. vφ[τ_]:=vφj[τ]/vtj[τ]*v[τ];
  57.  
  58. x0BL[A_] := Sqrt[r0^2+A^2] Sin[θ0] Cos[φ0]; (* Anfangskoordinaten *)
  59. y0BL[A_] := Sqrt[r0^2+A^2] Sin[θ0] Sin[φ0];
  60. z0[A_] := r0 Cos[θ0];
  61.  
  62. x0KS[A_] := ((r0 Cos[φ0]+A Sin[φ0]) Sin[θ0]);
  63. y0KS[A_] := ((r0 Sin[φ0]-A Cos[φ0]) Sin[θ0]);
  64.  
  65. x0[A_] := If[crd==1, x0BL[A], x0KS[A]];
  66. y0[A_] := If[crd==1, y0BL[A], y0KS[A]];
  67.  
  68. ε = Sqrt[δ Ξ/χ]/j[v0]+Lz ω0; (* Energie und Drehimpulskomponenten *)
  69. Lz = vφ0 Ы/j[v0];
  70. pθ0 = vθ0 Sqrt[Ξ]/j[v0];
  71. pr0 = vr0 Sqrt[(Ξ/δ)/j[v0]^2];
  72. Q = Simplify[Limit[pθ0^2+(Lz^2 Csc[ϑ]^2-a^2 (ε^2+μ)) Cos[ϑ]^2, ϑ->θ0]]; (* Carter Q *)
  73. k = Q+Lz^2+a^2 (ε^2+μ); (* Carter k *)
  74. μ = If[Abs[v0]==1, 0, -1]; (* Baryon: μ=-1, Photon: μ=0 *)
  75.  
  76.  
  77. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  78. (* |||||||| 3) RADIUS NACH GESCHWINDIGKEIT UND VICE VERSA ||||||||||||||||||||||||||||||| *)
  79. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  80.  
  81. rPro = 2 (1+Cos[2/3 ArcCos[-a]]); (* prograder Photonenorbitradius *)
  82. rRet = 2 (1+Cos[2/3 ArcCos[+a]]); (* retrograder Photonenorbitradius *)
  83. rTeo = 1+2 Sqrt[1-a^3/3] Cos[ArcCos[(1-a^2)/(1-a^2/3)^(3/2)]/3];
  84. (* ISCO *)
  85. isco = rISCO/.Solve[0 == rISCO (6 rISCO-rISCO^2+3 a^2)-8 a rISCO^(3/2) && rISCO>=rA, rISCO][[1]];
  86.  
  87. δp[r_,a_] := Quiet[δi/.NSolve[(a^4(-1+r)+2(-3+r)r^4+a^2r(6+r(-5+3 r))+ (* Eq. Ink. Winkel *)
  88. 4a Sqrt[a^2+(-2+r)r](a^2+3 r^2)Cos[δi]-a^2(3+r)(a^2+(-2+r)r)Cos[2δi])/(2r(a^2+
  89. (-2+r)r)(r^3+a^2(2+r)))==0&&δi<=π&&δi>=0,δi][[1]]];
  90.  
  91. vPro = (a^2-2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a+r0^(3/2)));(* Kreisgeschwindigkeit + *)
  92. vRet = (a^2+2a Sqrt[r0]+r0^2)/(Sqrt[a^2+(-2+r0)r0](a-r0^(3/2)));(* Kreisgeschwindigkeit - *)
  93. vr1 = \[Sqrt](((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2) ((a^2+r1^2)^2-a^2 (a^2+(-2+
  94. r1) r1) Sin[θ0]^2) (-1+((a^2+(-2+r1) r1) (r1^2+a^2 Cos[θ0]^2) (-(a^2+r0^2)^2+
  95. a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))/((a^2+(-2+r0) r0) (r0^2+a^2 Cos[θ0]^2) (-(a^2+
  96. r1^2)^2+a^2 (a^2+(-2+r1) r1) Sin[θ0]^2))))/((a^2+(-2+r1) r1) (r1^2+
  97. a^2 Cos[θ0]^2) ((a^2+r0^2)^2-a^2 (a^2+(-2+r0) r0) Sin[θ0]^2))); (* v Flucht von r0 bis r1 *)
  98.  
  99. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  100. (* |||||||| 4) HORIZONTE UND ERGOSPHÄREN RADIEN ||||||||||||||||||||||||||||||||||||||||| *)
  101. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  102.  
  103. rE = 1+Sqrt[1-a^2 Cos[θ]^2]; (* äußere Ergosphäre *)
  104. rG = 1-Sqrt[1-a^2 Cos[θ]^2]; (* innere Ergosphäre *)
  105. rA = 1+Sqrt[1-a^2]; (* äußerer Horizont *)
  106. rI = 1-Sqrt[1-a^2]; (* innerer Horizont *)
  107.  
  108. RE[A_, w1_, w2_] := Xyz[xyZ[
  109. {Sqrt[rE^2+A^2] Sin[θ] Cos[φ], Sqrt[rE^2+A^2] Sin[θ] Sin[φ], rE Cos[θ]}, w1], w2];
  110. RG[A_, w1_, w2_] := Xyz[xyZ[
  111. {Sqrt[rG^2+A^2] Sin[θ] Cos[φ], Sqrt[rG^2+A^2] Sin[θ] Sin[φ], rG Cos[θ]}, w1], w2];
  112. RA[A_, w1_, w2_] := Xyz[xyZ[
  113. {Sqrt[rA^2+A^2] Sin[θ] Cos[φ], Sqrt[rA^2+A^2] Sin[θ] Sin[φ], rA Cos[θ]}, w1], w2];
  114. RI[A_, w1_, w2_] := Xyz[xyZ[
  115. {Sqrt[rI^2+A^2] Sin[θ] Cos[φ], Sqrt[rI^2+A^2] Sin[θ] Sin[φ], rI Cos[θ]}, w1], w2];
  116.  
  117. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  118. (* |||||||| 5) HORIZONTE UND ERGOSPHÄREN PLOT ||||||||||||||||||||||||||||||||||||||||||| *)
  119. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  120.  
  121. horizons[A_, mesh_, w1_, w2_] := Show[
  122. ParametricPlot3D[RE[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  123. Mesh -> mesh, PlotPoints -> plp, PlotStyle -> Directive[Blue, Opacity[0.10]]],
  124. ParametricPlot3D[RA[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  125. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Cyan, Opacity[0.15]]],
  126. ParametricPlot3D[RI[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  127. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.25]]],
  128. ParametricPlot3D[RG[A, w1, w2], {φ, 0, 2 π}, {θ, 0, π},
  129. Mesh -> None, PlotPoints -> plp, PlotStyle -> Directive[Red, Opacity[0.35]]]];
  130. BLKS := Grid[{{horizons[a, 35, 0, 0], horizons[0, 35, 0, 0]}}];
  131.  
  132. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  133. (* |||||||| 6) FUNKTIONEN ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  134. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  135.  
  136. j[v_] := Sqrt[1+μ v^2]; (* Lorentzfaktor *)
  137. mirr = Sqrt[(Sqrt[1-a^2]+1)/2]; (* irreduzible Masse *)
  138. я = Sqrt[Χ/Σ]Sin[θ[τ]]; (* axialer Umfangsradius *)
  139. яi[τ_] := Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
  140.  
  141. Ы = Sqrt[χ/Ξ]Sin[θ0];
  142. Ыi[τ_] := Sqrt[Χi[τ]/Σi[τ]]Sin[Θ[τ]];
  143. ц = 2r[τ]/Σ; ц0=2r0/Ξ;
  144.  
  145. Σ = r[τ]^2+a^2 Cos[θ[τ]]^2; (* poloidialer Umfangsradius *)
  146. Σi[τ_] := R[τ]^2+a^2 Cos[Θ[τ]]^2;
  147. Ξ = r0^2+a^2 Cos[θ0]^2;
  148.  
  149. Δ = r[τ]^2-2r[τ]+a^2;
  150. Δi[τ_] := R[τ]^2-2R[τ]+a^2;
  151. δ = r0^2-2r0+a^2;
  152.  
  153. Χ = (r[τ]^2+a^2)^2-a^2 Sin[θ[τ]]^2 Δ;
  154. Χi[τ_] := (R[τ]^2+a^2)^2-a^2 Sin[Θ[τ]]^2 Δi[τ];
  155. χ = (r0^2+a^2)^2-a^2 Sin[θ0]^2 δ;
  156.  
  157. XJ = a Sin[θ[τ]]^2;
  158. xJ[τ_] := a Sin[Θ[τ]]^2;
  159. Xj = a Sin[θ0]^2;
  160.  
  161. т[τ_] := Evaluate[t[τ]/.sol][[1]]; (* Koordinatenzeit nach Eigenzeit *)
  162. д[ξ_] := Quiet[Ξ /.FindRoot[т[Ξ]-ξ, {Ξ, 0}]]; (* Eigenzeit nach Koordinatenzeit *)
  163. T := Quiet[д[tk]];
  164.  
  165. ю[τ_] := Evaluate[t'[τ]/.sol][[1]];
  166. γ[τ_] := If[μ==0, "Infinity", ю[τ]]; (* totale ZD *)
  167. R[τ_] := Evaluate[r[τ]/.sol][[1]]; (* Boyer-Lindquist Radius *)
  168. Φ[τ_] := Evaluate[φ[τ]/.sol][[1]];
  169. Θ[τ_] := Evaluate[θ[τ]/.sol][[1]];
  170. ß[τ_] := Re[Sqrt[X'[τ]^2+Y'[τ]^2+Z'[τ]^2 ]/ю[τ]];
  171.  
  172. gs[τ_] := (2 (R[τ]^2-a^2 Cos[Θ[τ]]^2) Sqrt[((a^2+R[τ]^2)^2-a^2 (a^2+(R[τ]- (* Gravitation *)
  173. 2) R[τ]) Sin[Θ[τ]]^2)/(a^2+2 R[τ]^2+a^2 Cos[2 Θ[τ]])^2])/(R[τ]^2+a^2 Cos[Θ[τ]]^2)^2;
  174.  
  175. ς[τ_] := Sqrt[Χi[τ]/Δi[τ]/Σi[τ]]; ς0 = Sqrt[χ/δ/Ξ]; (* gravitative ZD *)
  176. ω[τ_] := 2R[τ] a/Χi[τ]; ω0 = 2r0 a/χ; ωd=2r[τ] a/Χ; (* Frame Dragging *)
  177. Ω[τ_] := ω[τ] Sqrt[X[τ]^2+Y[τ]^2]; (* Frame Dragging beobachtete Geschwindigkeit *)
  178. й[τ_] := ω[τ] яi[τ] ς[τ]; й0 = ω0 Ы ς0; (* Frame Dragging lokale Geschwindigkeit *)
  179.  
  180. ж[τ_] := Sqrt[ς[τ]^2-1]/ς[τ]; ж0 = Sqrt[ς0^2-1]/ς0; (* Fluchtgeschwindigkeit *)
  181. vEsc = ж0;
  182.  
  183. vd[τ_] := Abs[-((\[Sqrt](-a^4(ε-Lz ωd)^2-2 a^2r[τ]^2 (ε-Lz ωd)^2-
  184. r[τ]^4(ε-Lz ωd)^2+Δ(Σ+a^2 Sin[θ[τ]]^2 (ε-
  185. Lz ωd)^2)))/(Sqrt[-(a^2+r[τ]^2)^2+
  186. a^2 Sin[θ[τ]]^2 Δ](ε - Lz ωd)))];
  187.  
  188. v[τ_] := If[μ==0, 1, Evaluate[vlt'[τ]/.sol][[1]]]; (* lokale Dreiergeschwindigkeit *)
  189. dst[τ_] := Evaluate[str[τ]/.sol][[1]]; (* Strecke *)
  190.  
  191. pΘ[τ_] := Evaluate[pθ[τ] /. sol][[1]]; pΘks[τ_] := Σi[τ] Θ'[τ]; (* Impuls *)
  192. pR[τ_] := Evaluate[pr[τ] /. sol][[1]]; pRks[τ_] := Σi[τ]/Δi[τ] R'[τ];
  193. sh[τ_] := Re[Sqrt[ß[τ]^2-Ω[τ]^2]];
  194. epot[τ_] := ε+μ-ekin[τ]; (* potentielle Energie *)
  195. ekin[τ_] := If[μ==0, ς[τ], 1/Sqrt[1-v[τ]^2]-1]; (* kinetische Energie *)
  196.  
  197. (* beobachtete Inklination *)
  198. ink0 := б/. Solve[Z'[0]/ю[0] Cos[б]==-Y'[0]/ю[0] Sin[б]&&б>0&&б<2π&&б<δp[r0, a], б][[1]];
  199.  
  200. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  201. (* |||||||| 7) DIFFERENTIALGLEICHUNG |||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  202. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  203.  
  204. dp= Style[\!\(\*SuperscriptBox[\(Y\),\(Y\)]\), White]; n0[z_] := Chop[Re[N[Simplify[z]]]];
  205.  
  206. (* Boyer-Lindquist-Koordinaten *)
  207.  
  208. pr2τ[τ_] := 1/(Σ Δ) (((r[τ]^2+a^2)μ-k)(r[τ]-1)+r[τ] Δ μ+2r[τ](r[τ]^2+a^2) ε^2-
  209. 2 a ε Lz)-(2pr[τ]^2 (r[τ]-1))/Σ;
  210. pθ2τ[τ_] := (Sin[θ[τ]]Cos[θ[τ]])/Σ (Lz^2/Sin[θ[τ]]^4-a^2 (ε^2+μ));
  211.  
  212. DG1={
  213. t'[τ] == ε+(2r[τ](r[τ]^2+a^2)ε-2 a r[τ] Lz)/(Σ Δ),
  214. t[0] == 0,
  215.  
  216. r'[τ] == (pr[τ] Δ)/Σ,
  217. r[0] == r0,
  218.  
  219. θ'[τ] == pθ[τ]/Σ,
  220. θ[0] == θ0,
  221.  
  222. φ'[τ] == (2 a r[τ] ε+(Σ-2r[τ])Lz Csc[θ[τ]]^2)/(Σ Δ),
  223. φ[0] == φ0,
  224.  
  225. pr'[τ] == pr2τ[τ],
  226. pr[0] == pr0,
  227.  
  228. pθ'[τ] == pθ2τ[τ],
  229. pθ[0] == pθ0,
  230.  
  231. str'[τ]== vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
  232. str[0] == 0,
  233.  
  234. vlt'[τ]== vd[τ],
  235. vlt[0] == 0
  236. };
  237.  
  238. (* Kerr-Schild-Koordinaten *)
  239.  
  240. dr = (pr0 δ)/Ξ; dθ=pθ0/Ξ;
  241. dφ = (2a r0 ε+(Ξ-2r0)Lz Csc[θ0]^2)/(Ξ δ)+dr a/δ;
  242. dΦ = If[θ0==0, 0, If[θ0==π, 0, dφ]];
  243. φk = φ0+cns/.FindRoot[Sqrt[r0^2+a^2] Cos[φ0+cns]-((r0 Cos[φ0]+a Sin[φ0])),{cns,1}];
  244.  
  245. DG2={
  246. t''[τ] == (2 ((a^4 Cos[θ[τ]]^4+a^2 Cos[θ[τ]]^2 r[τ]-r[τ]^3-r[τ]^4) r'[τ]^2+r[τ] ((a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-2 a^3 Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^3 θ'[τ] φ'[τ]+Sin[θ[τ]]^2 (r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]^2+a t'[τ] (a (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+2 (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ]))+r'[τ] ((a^4 Cos[θ[τ]]^4+2 a^2 Cos[θ[τ]]^2 r[τ]-2 r[τ]^3-r[τ]^4) t'[τ]+a (a r[τ] (2 a^2 Cos[θ[τ]]^3 Sin[θ[τ]]+r[τ]^2 Sin[2 θ[τ]]) θ'[τ]+(-a^4 Cos[θ[τ]]^4-2 a^2 Cos[θ[τ]]^2 r[τ]+2 r[τ]^3+r[τ]^4) Sin[θ[τ]]^2 φ'[τ]))))/(a^2 Cos[θ[τ]]^2+r[τ]^2)^3,
  247. t'[0] == Limit[(ц0 (-dr+a Sin[θ1]^2 dΦ))/(-1+ц0)+\[Sqrt]((1/((-1+ц0)^2 Ξ))(Ξ (dr^2+(-1+ц0) (μ-Ξ dθ^2))+Sin[θ1]^2 dΦ (-2a Ξ dr-(-1+ц0) χ dΦ+ц0^2 a^2 Ξ Sin[θ1]^2 dΦ))), θ1->θ0],
  248. t[0] == 0,
  249.  
  250. r''[τ] == (-8 (a^2 Cos[θ[τ]]^2-r[τ]^2) (a^2 Cos[2 θ[τ]]+r[τ] (2+r[τ])) r'[τ]^2+4 r'[τ] (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) (-2 r[τ]+a^2 Sin[θ[τ]]^2) t'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]-a Sin[θ[τ]]^2 (2 r[τ] (a^2 Cos[θ[τ]]^2 (-4+a^2+a^2 Cos[2 θ[τ]])+2 r[τ] ((2+a^2+a^2 Cos[2 θ[τ]]) r[τ]+r[τ]^3-a^2 Sin[θ[τ]]^2))+a^4 Sin[2 θ[τ]]^2) φ'[τ])+2 (a^2+(-2+r[τ]) r[τ]) (4 (a^2 Cos[θ[τ]]^2-r[τ]^2) t'[τ]^2+4 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2+8 a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 t'[τ] φ'[τ]+Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2))/(8 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3),
  251. r'[0] == dr,
  252. r[0] == r0,
  253.  
  254. θ''[τ] == (4 a^2 r[τ] Sin[2 θ[τ]] (r'[τ]+t'[τ])^2-8 r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 r'[τ] θ'[τ]+2 a^2 (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[2 θ[τ]] θ'[τ]^2-8 a ((Cos[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 Sin[θ[τ]]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]]) r'[τ]+r[τ] (a^2+r[τ]^2) Sin[2 θ[τ]] t'[τ]) φ'[τ]+(2 a^6 Cos[θ[τ]]^4+r[τ] (a^4 Cos[θ[τ]]^2 (5+Cos[2 θ[τ]]) r[τ]+2 a^2 (2+Cos[2 θ[τ]]) r[τ]^3+2 r[τ]^5+2 a^2 (a^2 (3+Cos[2 θ[τ]])+4 r[τ]^2) Sin[θ[τ]]^2)) Sin[2 θ[τ]] φ'[τ]^2)/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3),
  255. θ'[0] == dθ,
  256. θ[0] == θ0,
  257.  
  258. φ''[τ] == If[θ[τ]==0, 0, (4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) r'[τ]^2+4 (a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]^2+4 a r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2 θ'[τ]^2-8 (a^2 Cos[θ[τ]]^2+r[τ]^2) (Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 r[τ] Sin[2 θ[τ]]) θ'[τ] φ'[τ]+a Sin[θ[τ]]^2 (4 r[τ] ((a^2 Cos[θ[τ]]^2+r[τ]^2)^2-a^2 r[τ] Sin[θ[τ]]^2)+a^4 Sin[2 θ[τ]]^2) φ'[τ]^2+8 a t'[τ] (2 Cot[θ[τ]] r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2) θ'[τ]+a (-a^2 Cos[θ[τ]]^2+r[τ]^2) Sin[θ[τ]]^2 φ'[τ])+8 r'[τ] ((a^3 Cos[θ[τ]]^2-a r[τ]^2) t'[τ]+a Cot[θ[τ]] (a^2 Cos[θ[τ]]^2+r[τ]^2) (a^2 Cos[θ[τ]]^2+r[τ] (2+r[τ])) θ'[τ]-(r[τ] (a^2 Cos[θ[τ]]^2+r[τ]^2)^2+a^2 (a^2 Cos[θ[τ]]^2-r[τ]^2) Sin[θ[τ]]^2) φ'[τ]))/(4 (a^2 Cos[θ[τ]]^2+r[τ]^2)^3)],
  259. φ'[0] == dΦ,
  260. φ[0] == φk,
  261.  
  262. str'[τ]== vd[τ]/Max[1*^-16, Abs[Sqrt[1-vd[τ]^2]]],
  263. str[0] == 0,
  264.  
  265. vlt'[τ]== vd[τ],
  266. vlt[0] == 0
  267. };
  268.  
  269. DGL = If[crd==1, DG1, DG2];
  270.  
  271. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  272. (* |||||||| 8) INTEGRATION |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  273. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  274.  
  275. sol=
  276.  
  277. NDSolve[DGL, {t, r, θ, φ, str, vlt, pr, pθ}, {τ, 0, tmax},
  278. WorkingPrecision-> wp,
  279. MaxSteps-> Infinity,
  280. Method-> mta,
  281. InterpolationOrder-> All,
  282. StepMonitor :> (laststep=plunge; plunge=τ;
  283. stepsize=plunge-laststep;), Method->{"EventLocator",
  284. "Event" :> (If[stepsize<1*^-4, 0, 1])}];
  285.  
  286. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  287. (* |||||||| 9) KOORDINATEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  288. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  289.  
  290. XBL[τ_] := Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
  291. YBL[τ_] := Evaluate[Sqrt[r[τ]^2+a^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
  292. Z[τ_] := Evaluate[r[τ] Cos[θ[τ]]/.sol][[1]];
  293.  
  294. XKS[τ_] := Evaluate[((r[τ] Cos[φ[τ]]+a Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  295. YKS[τ_] := Evaluate[((r[τ] Sin[φ[τ]]-a Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  296.  
  297. X[τ_] := If[crd==1, XBL[τ], XKS[τ]];
  298. Y[τ_] := If[crd==1, YBL[τ], YKS[τ]];
  299.  
  300. xBL[τ_] := Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Cos[φ[τ]]/.sol][[1]];
  301. yBL[τ_] := Evaluate[Sqrt[r[τ]^2+A^2] Sin[θ[τ]] Sin[φ[τ]]/.sol][[1]];
  302. z[τ_] := Z[τ];
  303.  
  304. xKS[τ_] := Evaluate[((r[τ] Cos[φ[τ]]+A Sin[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  305. yKS[τ_] := Evaluate[((r[τ] Sin[φ[τ]]-A Cos[φ[τ]]) Sin[θ[τ]])/.sol][[1]];
  306.  
  307. x[τ_] := If[crd==1, xBL[τ], xKS[τ]];
  308. y[τ_] := If[crd==1, yBL[τ], yKS[τ]];
  309.  
  310. XYZ[τ_] := Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_] := Sqrt[X[τ]^2+Y[τ]^2]; (* kartesisches R *)
  311.  
  312. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z}; (* Rotationsmatrix *)
  313. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  314. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  315.  
  316. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  317. (* |||||||| 10) PLOT EINSTELLUNGEN |||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  318. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  319.  
  320. PR = 1.2 r1; (* Plot Range *)
  321. d1 = 2; (* Schweiflänge *)
  322. plp = 50; (* Flächenplot Details *)
  323. Plp = Automatic; (* Kurven Details *)
  324.  
  325. Mrec = 100; (* Parametric Plot Subdivisionen *)
  326. mrec = 10;
  327.  
  328. imgsize = 380; (* Bildgröße *)
  329. w1l = 0; (* Startperspektiven, Winkel *)
  330. w2l = 0;
  331. w1r = 0;
  332. w2r = 0;
  333.  
  334. s[text_]:=Style[text, FontFamily->"Consolas", FontSize->11]; (* Anzeigestil *)
  335.  
  336. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  337. (* |||||||| 11) PLOT NACH EIGENZEIT ||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  338. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  339.  
  340. display[T_] := Grid[{
  341. {If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[tp]], s["GM/c³"], s[dp]},
  342. {s[" t coord"], " = ", s[n0[т[tp]]], s["GM/c³"], s[dp]},
  343. (* {s[" т coord"], " = ", s[n0[т[tp]+Sign[1.5-dsp] 2 NIntegrate[R[Ti]/(R[Ti]^2-2 R[Ti]+a^2),{Ti,0,T}]]], s["GM/c³"], s[dp]}, *)
  344. {s[" ṫ total"], " = ", s[n0[If[dsp==1, ю[tp]-2 If[crd==1, 0, R'[tp] R[tp]/(R[tp]^2-2 R[tp]+a^2)], ю[T]]]], s[If[dsp==1, "dt/dτ", "dт/dτ"]], s[dp]},
  345. {s[" ς gravt"], " = ", s[n0[ς[tp]]], s["dt/dτ"], s[dp]},
  346. {s[" γ kinet"], " = ", If[μ==0, s[n0[ς[tp]]], s[n0[1/Sqrt[1-v[tp]^2]]]], s["dt/dτ"], s[dp]},
  347. {s[" R carts"], " = ", s[n0[XYZ[tp]]], s["GM/c²"], s[dp]},
  348. {s[" x carts"], " = ", s[n0[X[tp]]], s["GM/c²"], s[dp]},
  349. {s[" y carts"], " = ", s[n0[Y[tp]]], s["GM/c²"], s[dp]},
  350. {s[" z carts"], " = ", s[n0[Z[tp]]], s["GM/c²"], s[dp]},
  351. {s[" s dstnc"], " = ", s[n0[dst[tp]]], s["GM/c²"], s[dp]},
  352.  
  353. {s[" r coord"], " = ", s[n0[R[tp]]], s["GM/c²"], s[dp]},
  354. {s[" φ longd"], " = ", s[n0[Φ[tp] 180/π]], s["deg"], s[dp]},
  355. {s[" θ lattd"], " = ", s[n0[Θ[tp] 180/π]], s["deg"], s[dp]},
  356. {s[" d¹r/dτ¹"], " = ", s[n0[R'[tp]]], s["c"], s[dp]},
  357. {s[" d¹φ/dτ¹"], " = ", s[n0[Φ'[tp]]], s["c\.b3/G/M"], s[dp]},
  358. {s[" d¹θ/dτ¹"], " = ", s[n0[Θ'[tp]]], s["c\.b3/G/M"], s[dp]},
  359. {s[" d\.b2r/dτ\.b2"], " = ", s[n0[R''[tp]]], s["c⁴/G/M"], s[dp]},
  360. {s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Φ''[tp]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
  361. {s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Θ''[tp]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
  362. {s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
  363.  
  364. {s[" Ř crc.φ"], " = ", s[n0[яi[tp]]], s["GM/c²"], s[dp]},
  365. {s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[tp]]]], s["GM/c²"], s[dp]},
  366. {s[" M irred"], " = ", s[N[mirr]], s["M"], s[dp]},
  367. {s[" E kinet"], " = ", s[n0[ekin[tp]]], s["mc²"], s[dp]},
  368. {s[" E poten"], " = ", s[n0[epot[tp]]], s["mc²"], s[dp]},
  369. {s[" E total"], " = ", s[n0[ε]], s["mc²"], s[dp]},
  370. {s[" CarterQ"], " = ", s[n0[Q]], s["GMm/c"], s[dp]},
  371. {s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
  372. {s[" L polar"], " = ", s[n0[If[crd==1, pΘ[T], pΘks[T]]]], s["GMm/c"], s[dp]},
  373. {s[" p r.mom"], " = ", s[n0[If[crd==1, pR[tp], pRks[tp]]]], s["mc"], s[dp]},
  374.  
  375. {s[" ω fdrag"], " = ", s[n0[ω[tp]]], s["c³/G/M"], s[dp]},
  376. {s[" v fdrag"], " = ", s[n0[й[tp]]], s["c"], s[dp]},
  377. {s[" Ω fdrag"], " = ", s[n0[Ω[tp]]], s["c"], s[dp]},
  378. {s[" v propr"], " = ", s[n0[v[tp]/Sqrt[1-v[tp]^2]]], s["c"], s[dp]},
  379. {s[" v escpe"], " = ", s[n0[ж[tp]]], s["c"], s[dp]},
  380. {s[" v obsvd"], " = ", s[n0[ß[tp]]], s["c"], s[dp]},
  381. {s[" v r,loc"], " = ", s[n0[vr[tp]]], s["c"], s[dp]},
  382. {s[" v θ,loc"], " = ", s[n0[vθ[tp]]], s["c"], s[dp]},
  383. {s[" v φ,loc"], " = ", s[n0[vφ[tp]]], s["c"], s[dp]},
  384. {s[" v local"], " = ", s[n0[v[tp]]], s["c"], s[dp]},
  385. {s[" "], s[" "], s[" "], s[" "]}},
  386. Alignment-> Left, Spacings-> {0, 0}];
  387.  
  388. plot1b[{xx_, yy_, zz_, tk_, w1_, w2_}] := (* Animation *)
  389. Show[
  390.  
  391. Graphics3D[{
  392. {PointSize[0.011], Red, Point[
  393. Xyz[xyZ[{x[tp], y[tp], z[tp]}, w1], w2]]}},
  394. ImageSize-> imgsize,
  395. PlotRange-> {
  396. {-(2 Sign[Abs[xx]]+1) PR, +(2 Sign[Abs[xx]]+1) PR},
  397. {-(2 Sign[Abs[yy]]+1) PR, +(2 Sign[Abs[yy]]+1) PR},
  398. {-(2 Sign[Abs[zz]]+1) PR, +(2 Sign[Abs[zz]]+1) PR}
  399. },
  400. SphericalRegion->False,
  401. ImagePadding-> 1],
  402.  
  403. horizons[A, None, w1, w2],
  404.  
  405. If[A==0, {}, If[a==0, {}, ParametricPlot3D[
  406. Xyz[xyZ[{
  407. Sin[prm] a,
  408. Cos[prm] a,
  409. 0}, w1], w2],
  410. {prm, 0, 2π},
  411. PlotStyle -> {Thickness[0.005], Orange}]]],
  412.  
  413. If[crd==1, If[a==0, {},
  414. Graphics3D[{{PointSize[0.009], Purple, Point[
  415. Xyz[xyZ[{
  416. Sin[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  417. Cos[-φ0-ω0 т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  418. z0[A]}, w1], w2]]}}]],
  419. If[a==0, {},
  420. Graphics3D[{{PointSize[0.009], Purple, Point[
  421. Xyz[xyZ[{
  422. Sin[-φ0-ц0 a Ξ/χ т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  423. Cos[-φ0-ц0 a Ξ/χ т[tp]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  424. z0[A]}, w1], w2]]}}]]],
  425.  
  426. If[crd==1, If[tk==0, {}, If[a==0, {},
  427. ParametricPlot3D[
  428. Xyz[xyZ[{
  429. Sin[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  430. Cos[-φ0-ω0 т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  431. z0[A]}, w1], w2],
  432. {tt, Max[0, д[т[tp]-1/2 π/ω0]], tp},
  433. PlotStyle -> {Thickness[0.001], Dashed, Purple},
  434. PlotPoints-> Automatic,
  435. MaxRecursion-> 12]]],
  436. If[tk==0, {}, If[a==0, {},
  437. ParametricPlot3D[
  438. Xyz[xyZ[{
  439. Sin[-φ0-ц0 a Ξ/χ т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  440. Cos[-φ0-ц0 a Ξ/χ т[tt]+π/2] Sqrt[x0[A]^2+y0[A]^2],
  441. z0[A]}, w1], w2],
  442. {tt, Max[0, д[т[tp]-1/2 π/ω0]], tp},
  443. PlotStyle -> {Thickness[0.001], Dashed, Purple},
  444. PlotPoints-> Automatic,
  445. MaxRecursion-> 12]]]],
  446.  
  447. If[tk==0, {},
  448. Block[{$RecursionLimit = Mrec},
  449. ParametricPlot3D[
  450. Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp},
  451. PlotStyle-> {Thickness[0.004]},
  452. ColorFunction-> Function[{x, y, z, t},
  453. Hue[0, 1, 0.5, If[tp<0, Max[Min[(+tp+(-t+d1))/d1, 1], 0], Max[Min[(-tp+(t+d1))/d1, 1], 0]]]],
  454. ColorFunctionScaling-> False,
  455. PlotPoints-> Automatic,
  456. MaxRecursion-> mrec]]],
  457.  
  458. If[tk==0, {},
  459. Block[{$RecursionLimit = Mrec},
  460. ParametricPlot3D[
  461. Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]},
  462. PlotStyle-> {Thickness[0.004], Opacity[0.6], Darker@Gray},
  463. PlotPoints-> Plp,
  464. MaxRecursion-> mrec]]],
  465.  
  466. ViewPoint-> {xx, yy, zz}];
  467.  
  468. Quiet[Do[
  469. Print[Rasterize[Grid[{{
  470. plot1b[{0, -Infinity, 0, tp, w1l, w2l}],
  471. plot1b[{0, 0, +Infinity, tp, w1r, w2r}],
  472. display[tp]
  473. }, {" ", " ", " "}
  474. }, Alignment->Left]]],
  475. {tp, 0, tMax, tMax/1}]]
  476.  
  477. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  478. (* |||||||| 12) PLOT NACH KOORDINATENZEIT ||||||||||||||||||||||||||||||||||||||||||||||| *)
  479. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  480.  
  481. display[T_] := Grid[{
  482. {s[" t coord"], " = ", s[n0[tk]], s["GM/c³"], s[dp]},
  483. (* {s[" т coord"], " = ", s[n0[tk+Sign[1.5-dsp] 2 NIntegrate[R[Ti]/(R[Ti]^2-2 R[Ti]+a^2),{Ti,0,T}]]], s["GM/c³"], s[dp]}, *)
  484. {If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[T]], s["GM/c³"], s[dp]},
  485. {s[" ṫ total"], " = ", s[n0[If[dsp==1, ю[T]-2 If[crd==1, 0, R'[T] R[T]/(R[T]^2-2 R[T]+a^2)], ю[T]]]], s[If[dsp==1, "dt/dτ", "dт/dτ"]], s[dp]},
  486. {s[" ς gravt"], " = ", s[n0[ς[T]]], s["dt/dτ"], s[dp]},
  487. {s[" γ kinet"], " = ", If[μ==0, s[n0[ς[T]]], s[n0[1/Sqrt[1-v[T]^2]]]], s["dt/dτ"], s[dp]},
  488. {s[" R carts"], " = ", s[n0[XYZ[T]]], s["GM/c²"], s[dp]},
  489. {s[" x carts"], " = ", s[n0[X[T]]], s["GM/c²"], s[dp]},
  490. {s[" y carts"], " = ", s[n0[Y[T]]], s["GM/c²"], s[dp]},
  491. {s[" z carts"], " = ", s[n0[Z[T]]], s["GM/c²"], s[dp]},
  492. {s[" s dstnc"], " = ", s[n0[dst[T]]], s["GM/c²"], s[dp]},
  493.  
  494. {s[" r coord"], " = ", s[n0[R[T]]], s["GM/c²"], s[dp]},
  495. {s[" φ longd"], " = ", s[n0[Φ[T] 180/π]], s["deg"], s[dp]},
  496. {s[" θ lattd"], " = ", s[n0[Θ[T] 180/π]], s["deg"], s[dp]},
  497. {s[" d¹r/dτ¹"], " = ", s[n0[R'[T]]], s["c"], s[dp]},
  498. {s[" d¹φ/dτ¹"], " = ", s[n0[Φ'[T]]], s["c\.b3/G/M"], s[dp]},
  499. {s[" d¹θ/dτ¹"], " = ", s[n0[Θ'[T]]], s["c\.b3/G/M"], s[dp]},
  500. {s[" d\.b2r/dτ\.b2"], " = ", s[n0[R''[T]]], s["c⁴/G/M"], s[dp]},
  501. {s[" d\.b2φ/dτ\.b2"], " = ", s[n0[Φ''[T]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
  502. {s[" d\.b2θ/dτ\.b2"], " = ", s[n0[Θ''[T]]], s["c⁶/G\.b2/M\.b2"], s[dp]},
  503. {s[" a SpinP"], " = ", s[n0[a]], s["GM²/c"], s[dp]},
  504.  
  505. {s[" Ř crc.φ"], " = ", s[n0[яi[T]]], s["GM/c²"], s[dp]},
  506. {s[" Σ crc.θ"], " = ", s[n0[Sqrt[Σi[T]]]], s["GM/c²"], s[dp]},
  507. {s[" M irred"], " = ", s[N[mirr]], s["M"], s[dp]},
  508. {s[" E kinet"], " = ", s[n0[ekin[T]]], s["mc²"], s[dp]},
  509. {s[" E poten"], " = ", s[n0[epot[T]]], s["mc²"], s[dp]},
  510. {s[" E total"], " = ", s[n0[ε]], s["mc²"], s[dp]},
  511. {s[" CarterQ"], " = ", s[n0[Q]], s["GMm/c"], s[dp]},
  512. {s[" L axial"], " = ", s[n0[Lz]], s["GMm/c"], s[dp]},
  513. {s[" L polar"], " = ", s[n0[If[crd==1, pΘ[T], pΘks[T]]]], s["GMm/c"], s[dp]},
  514. {s[" p r.mom"], " = ", s[n0[If[crd==1, pR[T], pRks[T]]]], s["mc"], s[dp]},
  515.  
  516. {s[" ω fdrag"], " = ", s[n0[Abs[ω[T]]]], s["c³/G/M"], s[dp]},
  517. {s[" v fdrag"], " = ", s[n0[Abs[й[T]]]], s["c"], s[dp]},
  518. {s[" Ω fdrag"], " = ", s[n0[Abs[Ω[T]]]], s["c"], s[dp]},
  519. {s[" v propr"], " = ", s[n0[v[T]/Sqrt[1-v[T]^2]]], s["c"], s[dp]},
  520. {s[" v escpe"], " = ", s[n0[ж[T]]], s["c"], s[dp]},
  521. {s[" v obsvd"], " = ", s[n0[ß[T]]], s["c"], s[dp]},
  522. {s[" v r,loc"], " = ", s[n0[vr[T]]], s["c"], s[dp]},
  523. {s[" v θ,loc"], " = ", s[n0[vθ[T]]], s["c"], s[dp]},
  524. {s[" v φ,loc"], " = ", s[n0[vφ[T]]], s["c"], s[dp]},
  525. {s[" v local"], " = ", s[n0[v[T]]], s["c"], s[dp]},
  526. {s[" "], s[" "], s[" "], s[" "]}},
  527. Alignment-> Left, Spacings-> {0, 0}];
  528.  
  529. plot1a[{xx_, yy_, zz_, tk_, w1_, w2_}]:= (* Animation *)
  530. Show[
  531.  
  532. Graphics3D[{
  533. {PointSize[0.011], Red, Point[
  534. Xyz[xyZ[{x[T], y[T], z[T]}, w1], w2]]}},
  535. ImageSize-> imgsize,
  536. PlotRange-> {
  537. {-(2 Sign[Abs[xx]]+1) PR, +(2 Sign[Abs[xx]]+1) PR},
  538. {-(2 Sign[Abs[yy]]+1) PR, +(2 Sign[Abs[yy]]+1) PR},
  539. {-(2 Sign[Abs[zz]]+1) PR, +(2 Sign[Abs[zz]]+1) PR}
  540. },
  541. SphericalRegion->False,
  542. ImagePadding-> 1],
  543.  
  544. horizons[A, None, w1, w2],
  545.  
  546. If[A==0, {}, If[a==0, {}, ParametricPlot3D[
  547. Xyz[xyZ[{
  548. Sin[prm] a,
  549. Cos[prm] a,
  550. 0}, w1], w2],
  551. {prm, 0, 2π},
  552. PlotStyle -> {Thickness[0.005], Orange}]]],
  553.  
  554. If[crd==1, If[a==0, {},
  555. Graphics3D[{{PointSize[0.009], Purple, Point[
  556. Xyz[xyZ[{
  557. Sin[-φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
  558. Cos[-φ0-ω0 tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
  559. z0[A]}, w1], w2]]}}]],
  560. If[a==0, {},
  561. Graphics3D[{{PointSize[0.009], Purple, Point[
  562. Xyz[xyZ[{
  563. Sin[-φ0-ц0 a Ξ/χ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
  564. Cos[-φ0-ц0 a Ξ/χ tk+π/2] Sqrt[x0[A]^2+y0[A]^2],
  565. z0[A]}, w1], w2]]}}]]],
  566.  
  567. If[crd==1, If[tk==0, {}, If[a==0, {},
  568. ParametricPlot3D[
  569. Xyz[xyZ[{
  570. Sin[-φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
  571. Cos[-φ0-ω0 tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
  572. z0[A]}, w1], w2],
  573. {tt, Max[0, tk-1/2 π/ω0], tk},
  574. PlotStyle -> {Thickness[0.001], Dashed, Purple},
  575. PlotPoints-> Automatic,
  576. MaxRecursion-> mrec]]],
  577. If[tk==0, {}, If[a==0, {},
  578. ParametricPlot3D[
  579. Xyz[xyZ[{
  580. Sin[-φ0-ц0 a Ξ/χ tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
  581. Cos[-φ0-ц0 a Ξ/χ tt+π/2] Sqrt[x0[A]^2+y0[A]^2],
  582. z0[A]}, w1], w2],
  583. {tt, Max[0, tk-1/2 π/ω0], tk},
  584. PlotStyle -> {Thickness[0.001], Dashed, Purple},
  585. PlotPoints-> Automatic,
  586. MaxRecursion-> mrec]]]],
  587.  
  588. Block[{$RecursionLimit = Mrec},
  589. If[tk==0, {},
  590. ParametricPlot3D[
  591. Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, If[TMax<0, Min[0, T+d1], Max[0, T-d1]], T},
  592. PlotStyle-> {Thickness[0.004]},
  593. ColorFunction-> Function[{x, y, z, t},
  594. Hue[0, 1, 0.5, If[TMax<0, Max[Min[(+T+(-t+d1))/d1, 1], 0], Max[Min[(-T+(t+d1))/d1, 1], 0]]]],
  595. ColorFunctionScaling-> False,
  596. PlotPoints-> Automatic,
  597. MaxRecursion-> mrec]]],
  598.  
  599. If[tk==0, {},
  600. Block[{$RecursionLimit = Mrec},
  601. ParametricPlot3D[
  602. Xyz[xyZ[{x[tt], y[tt], z[tt]}, w1], w2], {tt, 0, If[Tmax<0, Min[-1*^-16, T+d1/3], Max[1*^-16, T-d1/3]]},
  603. PlotStyle-> {Thickness[0.004], Opacity[0.6], Darker@Gray},
  604. PlotPoints-> Plp,
  605. MaxRecursion-> mrec]]],
  606.  
  607. ViewPoint-> {xx, yy, zz}];
  608.  
  609. Quiet[Do[
  610. Print[Rasterize[Grid[{{
  611. plot1a[{0, -Infinity, 0, tk, w1l, w2l}],
  612. plot1a[{0, 0, Infinity, tk, w1r, w2r}],
  613. display[Quiet[д[tk]]]
  614. }, {" ", " ", " "}
  615. }, Alignment->Left]]],
  616. {tk, 0, TMax, TMax/1}]]
  617.  
  618. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  619. (* |||||||| 13) EXPORTOPTIONEN |||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  620. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  621.  
  622. (* Export als HTML Dokument *)
  623. (* Export["Y:\\export\\dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
  624. (* Export direkt als Bildsequenz *)
  625. (* Do[Export["Y:\\export\\dateiname" <> ToString[tk] <> ".png", Rasterize[...] *)
  626.  
  627. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  628. (* |||||||||||||||| http://kerr.yukerez.net ||||| Simon Tyran, Vienna ||||||||||||||||||| *)
  629. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
RAW Paste Data