Yukterez

Kerr Newman 3D Simulator, Doran, Raindrop

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