daily pastebin goal
52%
SHARE
TWEET

Schwarzschild Simulator

Yukterez Feb 13th, 2018 (edited) 74 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (*||| Schwarzschild Simulator |||||||| yukterez.net ||||||||||| Version 13.02.2018 |||||||| Syntax: Mathematica |||*)
  2. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  3.  
  4. ClearAll["Global`*"]                                       (* Variablen frei machen *)
  5. G = 1; M = 1; c = 1; rs = 2 G M/c^2;                       (* Einheiten *)
  6. wp = MachinePrecision;                                     (* numerische Genauigkeit *)
  7. set= {"GlobalAdaptive", "MaxErrorIncreases" -> 100, Method -> "GaussKronrodRule"};
  8. para = 20; pstep = 1/3;                                    (* Paraboloid Grid *)
  9.  
  10. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  11.  
  12. j[v_] := Sqrt[1 - v^2/c^2];    J = j[v0];                  (* Gammafaktor *)
  13. k[r_] := Sqrt[1 - rs/r];       κ = k[r0];                  (* Schwarzschildfaktor *)
  14.  
  15. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  16.  
  17. r0 = 3/2 rs;                                               (* Startradius *)
  18. v0 = 7/10 c;                                               (* lokale Startgeschwindigkeit *)
  19. φ  = Pi/4;                                                 (* Abschusswinkel *)
  20. θ0 = 0;                                                    (* Startwinkel *)
  21. τM = 20 G M/c^3;                                           (* Simulationsdauer *)
  22.  
  23. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  24.  
  25. τMax = Min[τM, plunge-0.0001]; TMax=Min[τM, и[plunge-0.0001]];
  26. (* г = Sqrt[χ^2 + γ^2]; Θ = ArcSin[γ/г]; *)                (* Kartesisch auf Polar *)
  27. vr0 = v0 Sin[φ] κ/J; vθ0 = v0/r0 Cos[φ]/J;                 (* Längenkontraktion und Tiefenexpansion *)
  28. d1 = τM/10; d2 = d1; f = 3;                                (* Schweifdauer und Frameanzahl *)
  29.  
  30. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  31.  
  32. sol = NDSolve[{                                            (* Differentialgleichung *)
  33. r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
  34. r'[0]  == vr0,
  35. r[0]   == r0,
  36. θ''[t] == -((2 r'[t] θ'[t])/r[t]),
  37. θ'[0]  == vθ0,
  38. θ[0]   == θ0,
  39. τ'[t]  == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
  40. τ[0]   == 0,
  41. cl'[t] == ((r'[t] / k[r[t]])^2 + (θ'[t] r[t])^2)/c^2,
  42. cl[0]  == 0
  43. }, {r, θ, τ, cl}, {t, 0, τM},
  44. MaxSteps -> Infinity,
  45. Method -> Automatic,
  46. WorkingPrecision -> wp,
  47. InterpolationOrder -> All,
  48. StepMonitor :> (laststep=plunge; plunge=t;
  49. stepsize=plunge-laststep;), Method->{"EventLocator",
  50. "Event" :> (If[stepsize<1*^-5, 0, 1])}];
  51.  
  52. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  53.  
  54. t[ξ_] := Quiet[                                            (* Eigenzeit nach Koordinatenzeit *)
  55. χ /.FindRoot[Evaluate[τ[χ] /. sol][[1]] - ξ, {χ, 0}, WorkingPrecision -> wp, Method -> Automatic]];
  56. Τ := Quiet[t[ι]];
  57.  
  58. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  59.  
  60. grr[r_]:=(1-2/r)^-1;
  61. Ȓ[я_]:=Quiet[NIntegrate[Sqrt[Abs[grr[r]]], {r, 0, Abs[я]}, Method -> set, MaxRecursion -> 100]];
  62. Ř[я_]:=Quiet[NIntegrate[Sqrt[Abs[grr[r]-1]], {r, Abs[я], para}, Method -> set, MaxRecursion -> 100]];
  63. ř[я_]:=Quiet[Ř[0]-Ř[я]];
  64. grid[я_]:=Quiet[x/.FindRoot[Ȓ[x]-я, {x, 1}]];
  65.  
  66. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  67.  
  68. x[t_] := (Sin[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
  69. y[t_] := (Cos[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
  70.  
  71. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  72.  
  73. R[t_] := Evaluate[r[t] /. sol][[1]];                       (* Radialkoordinate *)
  74. γ[t_] := Evaluate[τ'[t] /. sol][[1]];                      (* Zeitdilatation *)
  75. и[t_] := Evaluate[τ[t] /. sol][[1]];                       (* Koordinatenzeit *)
  76.  
  77. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  78.  
  79. crθ[t_] := Evaluate[cl'[t] /. sol][[1]];                   (* Celerität *)
  80. vrθ[t_] := crθ[t]/Sqrt[1 + crθ[t]^2];
  81. clr[t_] := Evaluate[r'[t] /. sol][[1]];
  82. clθ[t_] := R[t] Evaluate[θ'[t] /. sol][[1]];
  83.  
  84. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  85.  
  86. vr[t_] := clr[t]/γ[t]/k[R[t]]^2;                           (* lokale Geschwindigkeit, radial *)
  87. vt[t_] := clθ[t]/γ[t]/k[R[t]];                             (* lokale Geschwindigkeit, transversal *)
  88. vp[t_] := Sqrt[vr[t]^2 + vt[t]^2];                         (* lokale Geschwindigkeit, total *)
  89. vc[t_] := Sqrt[vr[t]^2 k[R[t]]^2 + vt[t]^2] k[R[t]];       (* Koordinatengeschwindigkeit, total *)
  90.  
  91. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  92.  
  93. s[text_] := Style[text, FontSize -> font];  font = 11;     (* Stil der Anzeigetafel *)
  94. PR = 2 r0;                                                 (* Plot Range *)
  95.  
  96. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  97.  
  98. Do[Print[                                                  (* Animation nach Eigenzeit *)
  99. Rasterize[Grid[{{Show[
  100.  
  101. Graphics[{Table[{LightGray, Circle[{0, 0}, grid[n]]}, {n, pstep, para, pstep}]},
  102. Frame -> True, ImageSize -> 360,
  103. PlotRange -> PR, ImagePadding ->  Automatic],
  104. Graphics[{Magenta, Opacity[0.1], Annulus[{0, 0}, {2, para}]}],
  105. Graphics[{White, Opacity[0.7], Disk[{0, 0}, 2]}],
  106. Graphics[{Cyan, Opacity[0.1], Disk[{0, 0}, 2]}],
  107. Graphics[{Black, Dashed, Circle[{0, 0}, r0]}],
  108.  
  109. Graphics[{PointSize[0.01], Red, Point[{x[т], y[т]}]}],
  110.  
  111. ParametricPlot[{x[η], y[η]}, {η, 0, Max[1*^-16, т - d1]},
  112. PlotStyle->{Thickness[0.004], LightGray}],
  113. ParametricPlot[{x[η], y[η]}, {η, Max[1*^-16, т - d2], т},
  114. ColorFunction -> Function[{x, y, η},
  115. Hue[0, 1, 0.5, Max[Min[(-т + (η + d2))/d2, 1], 0]]],
  116. ColorFunctionScaling -> False]]},
  117.  
  118.       {Grid[{
  119.       {s["  Eigenzeit"],           " = ",    s[N[т, 8]],                                      s[" GM/c³"]},
  120.       {s["  Koordinatenzeit"],     " = ",    s[N[Evaluate[τ[т] /. sol][[1]], 8]],             s[" GM/c³"]},
  121.       {s["  Zeitdilatation"],      " = ",    s[N[γ[т], 8]],                                   s[" dt/dτ"]},
  122.       {s["  Winkel"],              " = ",    s[N[Evaluate[(θ[т] /. sol) 180/Pi][[1]], 8]],    s[" Grad"]},
  123.       {s["  Radialkoordinate"],    " = ",    s[N[R[т] , 8]],                                  s[" GM/c²"]},
  124.       {s["  x-Achse"],             " = ",    s[N[x[т], 8]],                                   s[" GM/c²"]},
  125.       {s["  y-Achse"],             " = ",    s[N[y[т], 8]],                                   s[" GM/c²"]},
  126.       {s["  v lokal"],             " = ",    s[N[vp[т], 8]],                                  s[" c"]},
  127.       {s["  v extern"],            " = ",    s[N[vc[т], 8]],                                  s[" c"]},
  128.       {s["  kinetische Energie"],  " = ",    s[N[1/Sqrt[1 - vp[т]^2] - 1, 8]],                s[" mc²"]},
  129.       {s[" "],                     "   ",    s["                 "],                          s[" "]}
  130.       }, Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]
  131.       ], {т, τMax/f, τMax, τMax/f}]
  132.    
  133. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  134.  
  135. Do[Print[                                                  (* Animation nach Koordinatenzeit *)
  136. Rasterize[Grid[{{Show[
  137.  
  138. Graphics[{Table[{LightGray, Circle[{0, 0}, grid[n]]}, {n, pstep, para, pstep}]},
  139. Frame -> True, ImageSize -> 360,
  140. PlotRange -> PR, ImagePadding ->  Automatic],
  141. Graphics[{Magenta, Opacity[0.1], Annulus[{0, 0}, {2, para}]}],
  142. Graphics[{White, Opacity[0.7], Disk[{0, 0}, 2]}],
  143. Graphics[{Cyan, Opacity[0.1], Disk[{0, 0}, 2]}],
  144. Graphics[{Black, Dashed, Circle[{0, 0}, r0]}],
  145.  
  146. Graphics[{PointSize[0.01], Red, Point[{x[Τ], y[Τ]}]}],
  147.  
  148. ParametricPlot[{x[η], y[η]}, {η, 0, Max[1*^-16, Τ - d1]},
  149. PlotStyle->{Thickness[0.004], LightGray}],
  150. ParametricPlot[{x[η], y[η]}, {η, Max[1*^-16, Τ - d2], Τ},
  151. ColorFunction -> Function[{x, y, η},
  152. Hue[0, 1, 0.5, Max[Min[(-Τ + (η + d2))/d2, 1], 0]]],
  153. ColorFunctionScaling -> False]]},
  154.  
  155.       {Grid[{
  156.       {s["  Eigenzeit"],           " = ",    s[N[Τ, 8]],                                      s[" GM/c³"]},
  157.       {s["  Koordinatenzeit"],     " = ",    s[N[ι, 8]],                                      s[" GM/c³"]},
  158.       {s["  Zeitdilatation"],      " = ",    s[N[γ[Τ], 8]],                                   s[" dt/dτ"]},
  159.       {s["  Winkel"],              " = ",    s[N[Evaluate[(θ[Τ] /. sol) 180/Pi][[1]], 8]],    s[" Grad"]},
  160.       {s["  Radialkoordinate"],    " = ",    s[N[R[Τ], 8]],                                   s[" GM/c²"]},
  161.       {s["  x-Achse"],             " = ",    s[N[x[Τ], 8]],                                   s[" GM/c²"]},
  162.       {s["  y-Achse"],             " = ",    s[N[y[Τ], 8]],                                   s[" GM/c²"]},
  163.       {s["  v lokal"],             " = ",    s[N[vp[Τ], 8]],                                  s[" c"]},
  164.       {s["  v extern"],            " = ",    s[N[vc[Τ], 8]],                                  s[" c"]},
  165.       {s["  kinetische Energie"],  " = ",    s[N[1/Sqrt[1 - vp[Τ]^2] - 1, 8]],                s[" mc²"]},
  166.       {s[" "],                     "   ",    s["                 "],                          s[" "]}
  167.       }, Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]
  168.       ], {ι, TMax/f, TMax, TMax/f}]
  169.    
  170. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  171. (* Syntax: Mathematica || schwarzschild.yukterez.net || kerr.yukterez.net ||| Simon Tyran - [Симон Тыран] - Vienna *)
RAW Paste Data
Top