# 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,
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