Yukterez

Kerr-Schild, Photon Vs Ring Singularity

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