Yukterez

Schwarzschild Plunge

Feb 17th, 2018
157
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (*||| Schwarzschild Simulator |||||||| yukterez.net ||||||||||| Simplified Version |||||||| 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. τMax = Min[τM, plunge-0.000001]; (* Integrationsende *)
  9. para = 20; pstep = 1/3; (* Paraboloid Grid *)
  10.  
  11. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  12.  
  13. j[v_] := Sqrt[1 - v^2/c^2]; J = j[v0]; (* Gammafaktor *)
  14. k[r_] := Sqrt[1 - rs/r]; κ = k[r0]; (* Schwarzschildfaktor *)
  15.  
  16. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  17.  
  18. r0 = 3/2 rs; (* Startradius *)
  19. v0 = 7/10 c; (* lokale Startgeschwindigkeit *)
  20. φ = Pi/4; (* Abschusswinkel *)
  21. θ0 = 0; (* Startwinkel *)
  22. τM = 18 G M/c^3; (* Simulationsdauer *)
  23.  
  24. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  25.  
  26. vr0 = v0 Sin[φ] κ/J; vθ0 = v0/r0 Cos[φ]/J; (* Längenkontraktion und Tiefenexpansion *)
  27. d1 = τM; d2 = d1; f = 1; (* Schweifdauer und Frameanzahl *)
  28.  
  29. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  30.  
  31. sol = NDSolve[{ (* Differentialgleichung *)
  32. r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
  33. r'[0] == vr0,
  34. r[0] == r0,
  35. θ''[t] == -((2 r'[t] θ'[t])/r[t]),
  36. θ'[0] == vθ0,
  37. θ[0] == θ0
  38. }, {r, θ}, {t, 0, τM},
  39. MaxSteps -> Infinity,
  40. Method -> Automatic,
  41. WorkingPrecision -> wp,
  42. InterpolationOrder -> All,
  43. StepMonitor :> (laststep=plunge; plunge=t;
  44. stepsize=plunge-laststep;), Method->{"EventLocator",
  45. "Event" :> (If[stepsize<1*^-5, 0, 1])}];
  46.  
  47. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  48.  
  49. grr[r_]:=(1-2/Abs[r])^-1;
  50. R[t_] := Evaluate[r[t] /. sol][[1]];
  51. Ȓ[я_]:=Quiet[NIntegrate[Sqrt[Abs[grr[r]]], {r, 0, Abs[я]}, Method -> set, MaxRecursion -> 100]];
  52. Ř[я_]:=Quiet[NIntegrate[Sqrt[Abs[grr[r]-1]], {r, Abs[я], para}, Method -> set, MaxRecursion -> 100]];
  53. ř[я_]:=Quiet[Ř[0]-Ř[я]];
  54. grid[я_]:=Quiet[x/.FindRoot[Ȓ[x]-я, {x, 1}]];
  55.  
  56. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  57.  
  58. x[t_] := (Sin[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
  59. y[t_] := (Cos[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
  60.  
  61. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  62.  
  63. s[text_] := Style[text, FontSize -> font]; font = 11; (* Stil der Anzeigetafel *)
  64. PR = 2 r0; (* Plot Range *)
  65.  
  66. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  67.  
  68. Quiet[Do[Print[ (* Animation nach Eigenzeit *)
  69. Rasterize[Grid[{{Show[
  70. Graphics[{Table[{LightGray, Circle[{0, 0}, grid[n]]}, {n, pstep, para, pstep}]},
  71. Frame -> True, ImageSize -> 360,
  72. PlotRange -> PR, ImagePadding -> Automatic],
  73. Graphics[{Magenta, Opacity[0.1], Annulus[{0, 0}, {2, para}]}],
  74. Graphics[{White, Opacity[0.7], Disk[{0, 0}, 2]}],
  75. Graphics[{Cyan, Opacity[0.1], Disk[{0, 0}, 2]}],
  76. Graphics[{Black, Dashed, Circle[{0, 0}, r0]}],
  77. Graphics[{PointSize[0.01], Red, Point[{x[т], y[т]}]}],
  78. ParametricPlot[{x[η], y[η]}, {η, Max[1*^-16, т - d2], т},
  79. ColorFunction -> Function[{x, y, η},
  80. Hue[0, 1, 0.5, Max[Min[(-т + (η + d2))/d2, 1], 0]]],
  81. ColorFunctionScaling -> False]]},
  82. {Grid[{
  83. {s[" Eigenzeit"], " = ", s[N[т, 8]], s[" GM/c³"]},
  84. {s[" Radialkoordinate"], " = ", s[N[R[т] , 8]], s[" GM/c²"]},
  85. {s[" Winkelkoordinate"], " = ", s[N[Evaluate[(θ[т] /. sol) 180/Pi][[1]], 8]], s[" Grad°"]},
  86. {s[" "], " ", s[" "], s[" "]}
  87. }, Alignment -> Left, Spacings -> {0, 1/2}]}}, Alignment -> Left]]
  88. ], {т, τMax/f, τMax, τMax/f}]]
  89.  
  90. (*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*)
  91. (*||||| Syntax: Mathematica ||||| schwarzschild.yukterez.net |||| kerr.yukterez.net ||||| Simon Tyran, Vienna |||||*)
RAW Paste Data