Yukterez

Gödel Metric Simulator

Jan 21st, 2020
118
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* ||| Gödel Metric | Mathematica Syntax | geodesics.yukterez.net | Version 22.1.2020 ||| *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. ClearAll["Local`*"];
  6. Needs["DifferentialEquations`NDSolveProblems`"];
  7. Needs["DifferentialEquations`NDSolveUtilities`"];
  8.  
  9. wp=MachinePrecision;
  10. mt1=Automatic;
  11. mt2={"EquationSimplification"-> "Residual"};
  12. mt3={"ImplicitRungeKutta", "DifferenceOrder"-> 20};
  13. mt4={"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
  14. mta=mt1; (* mt1: Speed, mt3: Accuracy *)
  15. dgl=1;
  16.  
  17. tmax = 50; (* Eigenzeit *)
  18. Tmax = 50; (* Koordinatenzeit *)
  19. TMax = Min[Tmax, т[plunge-1*^-3]]; tMax = Min[tmax, plunge]; (* Integrationsende *)
  20.  
  21. ω=1/2; (* Rotationsparameter *)
  22.  
  23. r0=Sqrt[2]/ω; (* Anfangsposition *)
  24. φ0=0;
  25. z0=-1;
  26.  
  27. vr=Sqrt[89/2]/10; (* Anfangsgeschwindigkeit *)
  28. vφ=Sqrt[89/2]/10;
  29. vz=Sqrt[11]/10;
  30.  
  31. v0=Sqrt[vr^2+vφ^2+vz^2]; "v"->N@v0
  32. μ=If[Abs[v0]==1,0,-1]; (* Partikel: μ=-1, Photon: μ=0 *)
  33.  
  34. dp= Style[\!\(\*SuperscriptBox[\(Y\),\(Y\)]\), White]; n0[z_]:=Chop[N[Simplify[z]]];
  35.  
  36. dt=Quiet[Sqrt[(2-ω^2 r0^2)/(2+ω^2 r0^2)]/Sqrt[1-v0^2 μ^2]];
  37. dr=Quiet[vr/(Sqrt[2] Sqrt[1-v0^2 μ^2] Sqrt[1/(2+r0^2 ω^2)])];
  38. dφ=Quiet[(Sqrt[4 r0^2-2 r0^4 ω^2] (-vφ-(2 dt r0^2 Sqrt[1-v0^2 μ^2] ω)/Sqrt[4 r0^2-
  39. 2 r0^4 ω^2]))/(r0^2 Sqrt[1-v0^2 μ^2] (-2+r0^2 ω^2))];
  40. dz=Quiet[vz/Sqrt[1-v0^2 μ^2]];
  41.  
  42. DGL=Quiet[{
  43.  
  44. t''[τ]==-((2 ω^2 r[τ] r'[τ] (2 t'[τ]+ω r[τ]^2 φ'[τ]))/(2+ω^2 r[τ]^2)),
  45. t'[0]==dt,
  46. t[0]==0,
  47.  
  48. r''[τ]==-((r[τ] (-2 ω^2 r'[τ]^2+(2+ω^2 r[τ]^2)^2 φ'[τ] (2 ω t'[τ]+(-1+
  49. ω^2 r[τ]^2) φ'[τ])))/(4+2 ω^2 r[τ]^2)),
  50. r'[0]==dr,
  51. r[0]==r0,
  52.  
  53. φ''[τ]==+(4 r'[τ] (ω t'[τ]-φ'[τ]))/(r[τ] (2+ω^2 r[τ]^2)),
  54. φ'[0]==dφ,
  55. φ[0]==φ0,
  56.  
  57. z''[τ]==0,
  58. z'[0]==dz,
  59. z[0]==z0
  60.  
  61. }];
  62.  
  63. sol=Quiet[NDSolve[DGL, {t, r, φ, z}, {τ, 0, tmax},
  64. WorkingPrecision-> wp,
  65. MaxSteps-> Infinity,
  66. Method-> mta,
  67. InterpolationOrder-> All,
  68. StepMonitor :> (laststep=plunge; plunge=τ;
  69. stepsize=plunge-laststep;), Method->{"EventLocator",
  70. "Event" :> (If[stepsize<1*^-4, 0, 1])}]];
  71.  
  72. T[τ_]:=Chop[Evaluate[t[τ]/.sol][[1]]];
  73. R[τ_]:=Chop[Evaluate[r[τ]/.sol][[1]]];
  74. Φ[τ_]:=Chop[Evaluate[φ[τ]/.sol][[1]]];
  75. Z[τ_]:=Chop[Evaluate[z[τ]/.sol][[1]]];
  76. Y[τ_]:=Chop[R[τ] Sin[Φ[τ]]];
  77. X[τ_]:=Chop[R[τ] Cos[Φ[τ]]];
  78.  
  79. XYZ[τ_]:=Sqrt[X[τ]^2+Y[τ]^2+Z[τ]^2]; XY[τ_]:=Sqrt[X[τ]^2+Y[τ]^2]; (* kartesischer Radius *)
  80.  
  81. Xyz[{x_, y_, z_}, α_]:={x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z}; (* Rotationsmatrix *)
  82. xYz[{x_, y_, z_}, β_]:={x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  83. xyZ[{x_, y_, z_}, ψ_]:={x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  84.  
  85. v[τ_]:=If[μ==0, 1, (* Lokalgeschwindigkeit *)
  86. Sqrt[-2+2 T'[τ]^2+ω^2 R[τ]^2+T'[τ]^2 ω^2 R[τ]^2]/(T'[τ] Sqrt[2+ω^2 R[τ]^2])];
  87.  
  88. PR=2 r0; (* Plot Range *)
  89. d1=tMax/10; (* Schweiflänge *)
  90. plp=50; (* Flächenplot Details *)
  91. Plp=Automatic; (* Kurven Details *)
  92.  
  93. w1l=0; w2l=0; w1r=0; w2r=0; (* Startperspektiven *)
  94. Mrec=100; mrec=10; (* Parametric Plot Subdivisionen *)
  95. imgsize=380; (* Bildgröße *)
  96.  
  97. s[text_]:=Style[text, FontFamily->"Consolas", FontSize->11]; (* Anzeigestil *)
  98.  
  99. display1[tp_]:=Quiet[Grid[{
  100. {},
  101. {If[μ==0, s[" affineP"], s[" τ propr"]], " = ", s[n0[tp]], s["GM/c³"], s[dp]},
  102. {s[" t coord"], " = ", s[n0[T[tp]]], s["GM/c³"], s[dp]},
  103. {s[" ṫ total"], " = ", s[n0[T'[tp]]], s["dt/dτ"], s[dp]},
  104. {},
  105. {s[" x carts"], " = ", s[n0[X[tp]]], s["GM/c²"], s[dp]},
  106. {s[" y carts"], " = ", s[n0[Y[tp]]], s["GM/c²"], s[dp]},
  107. {s[" z carts"], " = ", s[n0[Z[tp]]], s["GM/c²"], s[dp]},
  108.  
  109. {s[" "], s[" "], s[" "], s[" "]}},
  110. Alignment-> Left, Spacings-> {0, 0}]];
  111.  
  112. display2[tp_]:=Quiet[Grid[{
  113. {},
  114. {s[" d¹x/dτ¹"], " = ", s[n0[X'[tp]]], s["c"], s[dp]},
  115. {s[" d¹y/dτ¹"], " = ", s[n0[Y'[tp]]], s["c"], s[dp]},
  116. {s[" d¹z/dτ¹"], " = ", s[n0[Z'[tp]]], s["c"], s[dp]},
  117. {},
  118. {s[" r critc"], " = ", s[n0[Sqrt[2]/ω]], s["GM/c²"], s[dp]},
  119. {s[" r cylnd"], " = ", s[n0[R[tp]]], s["GM/c²"], s[dp]},
  120. {s[" φ cylnd"], " = ", s[n0[Φ[tp] 180/π]], s["degree"], s[dp]},
  121.  
  122. {s[" "], s[" "], s[" "], s[" "]}},
  123. Alignment-> Left, Spacings-> {0, 0}]];
  124.  
  125. display3[tp_]:=Quiet[Grid[{
  126. {},
  127. {s[" d²r/dτ²"], " = ", s[n0[R''[tp]]], s["c⁶/G/M"], s[dp]},
  128. {s[" d²φ/dτ²"], " = ", s[n0[Φ''[tp]]], s["c⁶/G²/M²"], s[dp]},
  129. {},
  130. {s[" d¹r/dτ¹"], " = ", s[n0[R'[tp]]], s["c"], s[dp]},
  131. {s[" d¹φ/dτ¹"], " = ", s[n0[Φ'[tp] 180/π]], s["c³/G/M"], s[dp]},
  132. {},
  133. {s[" v local"], " = ", s[n0[v[tp]]], s["c"], s[dp]},
  134.  
  135. {s[" "], s[" "], s[" "], s[" "]}},
  136. Alignment-> Left, Spacings-> {0, 0}]];
  137.  
  138. plot1b[{X_, Y_, Z_, xx_, yy_, zz_, tp_, w1_, w2_}]:= (* Animation *)
  139. Quiet[Show[
  140.  
  141. Graphics3D[{
  142. {PointSize[0.015], Red, Point[
  143. Xyz[xyZ[{X[tp], Y[tp], Z[tp]}, w1], w2]]}},
  144. ImageSize-> imgsize,
  145. PlotRange-> {
  146. {-(9 Sign[Abs[xx]]+1) PR, +(9 Sign[Abs[xx]]+1) PR},
  147. {-(9 Sign[Abs[yy]]+1) PR, +(9 Sign[Abs[yy]]+1) PR},
  148. {-(9 Sign[Abs[zz]]+1) PR, +(9 Sign[Abs[zz]]+1) PR}
  149. },
  150. SphericalRegion->False,
  151. ImagePadding-> 1],
  152.  
  153. Graphics3D[{
  154. {PointSize[0.01], Green, Point[
  155. Xyz[xyZ[{X[0], Y[0], Z[0]}, w1], w2]]}}],
  156.  
  157. Graphics3D[{
  158. {PointSize[0.01], Blue, Point[
  159. Xyz[xyZ[{0, 0, 0}, w1], w2]]}}],
  160.  
  161. If[tp==0, {},
  162. Block[{$RecursionLimit = Mrec},
  163. ParametricPlot3D[
  164. Xyz[xyZ[{X[tt], Y[tt], Z[tt]}, w1], w2], {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp},
  165. PlotStyle-> {Thickness[0.004]},
  166. ColorFunction-> Function[{x, y, z, t},
  167. Hue[0, 1, 0.5, If[tp<0,
  168. Max[Min[(+tp+(-t+d1))/d1, 1], 0], Max[Min[(-tp+(t+d1))/d1, 1], 0]]]],
  169. ColorFunctionScaling-> False,
  170. PlotPoints-> Automatic,
  171. MaxRecursion-> mrec]]],
  172.  
  173. If[tp==0, {},
  174. Block[{$RecursionLimit = Mrec},
  175. ParametricPlot3D[
  176. Xyz[xyZ[{X[tt], Y[tt], Z[tt]}, w1], w2],
  177. {tt, 0, If[tp<0, Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]},
  178. PlotStyle-> {Thickness[0.004], Opacity[0.6], Darker[Gray]},
  179. PlotPoints-> Plp,
  180. MaxRecursion-> mrec]]],
  181.  
  182. ViewPoint-> {xx, yy, zz}]];
  183.  
  184. Quiet[Do[
  185. Print[Rasterize[Grid[{{
  186. plot1b[{X, Y, T, 0, -Infinity, 0, tp, w1l, w2l}],
  187. plot1b[{X, Y, Z, 0, -Infinity, 0, tp, w1l, w2l}],
  188. plot1b[{X, Y, Z, 0, 0, +Infinity, tp, w1r, w2r}]},
  189. {"x,t", "x,z", "x,y"},
  190. {display1[tp], display2[tp], display3[tp]}
  191. }, Alignment->Left]]],
  192. {tp, 0, tMax, tMax/2}]]
  193.  
  194. "r of τ"->Quiet[Plot[R[τ], {τ, 0, tMax}, Frame->True]]
  195. "ṫ of τ"->Quiet[Plot[T'[τ], {τ, 0, tMax}, Frame->True]]
  196. "t of τ"->Quiet[Plot[T[τ], {τ, 0, tMax}, Frame->True]]
RAW Paste Data