Yukterez

Gravity of a Planet with a Disk or Ring, Simulator

Apr 24th, 2020
62
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* ||||| Mathematica Syntax | yukterez.net | Planet, Ring & Disk Simulator | Version 2 ||||||| *)
  3. (* ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
  6. mt2 = {"ImplicitRungeKutta", "DifferenceOrder"-> 20};
  7. mt3 = {"EquationSimplification"-> "Residual"};
  8. mt0 = Automatic;
  9. mta = mt1;
  10. wp = MachinePrecision;
  11.  
  12. T = 2000; (* Simulationsdauer *)
  13.  
  14. G = 1; (* Gravitationskonstante *)
  15. M = 1/2; (* Masse der Scheibe *)
  16. Ḿ = 1; (* Masse der zentralen Kugel *)
  17.  
  18. я1 = 15; (* Scheibeninnenradius *)
  19. я2 = 20; (* Scheibenaußenradius *)
  20. я3 = 10; (* Kugelradius *)
  21.  
  22. x0 = 25; (* Startposition x *)
  23. y0 = 0; (* Startposition y *)
  24. z0 = 1/100000; (* Startposition z *)
  25.  
  26. v0 = Sqrt[vx^2+vy^2+vz^2]; (* Anfangsgeschwindigkeit *)
  27.  
  28. vx = 0; (* Anfangsgeschwindigkeit x *)
  29. vy = 0; (* Anfangsgeschwindigkeit y *)
  30. vz = 900/1019 Sqrt[G (M+Ḿ)/x0]; (* Anfangsgeschwindigkeit z *)
  31.  
  32. ρ = M/(π я2^2-π я1^2); (* Flächendichte der Scheibe *)
  33. m = If[я3==0, Ḿ, Ḿ Sqrt[x[t]^2+y[t]^2+z[t]^2]^3/я3]; (* innere Kugelrestmasse *)
  34.  
  35. r[t_] := Sqrt[x[t]^2+y[t]^2]; (* zylindrischer Radius *)
  36. R[t_] := Sqrt[x[t]^2+y[t]^2+z[t]^2]; (* sphärischer Radius *)
  37.  
  38. k[t_, я_] := Sqrt[4r[t] я/(R[t]^2+я^2+2r[t] я)]; (* Funktionen *)
  39. α[t_, я_] := Sqrt[4 r[t] я/(я+r[t])^2];
  40. φ[t_, я_] := ArcSin[z[t]/(R[t]^2+я^2-2r[t] я)];
  41.  
  42. ax[t_, я_] := (-4G ρ x[t] Sqrt[я]/(k[t, я] r[t]^(3/2)))((1-
  43. k[t, я]^2/2) EllipticK[k[t, я]^2]-EllipticE[k[t, я]^2]);
  44.  
  45. ay[t_, я_] := (-4G ρ y[t] Sqrt[я]/(k[t, я] r[t]^(3/2)))((1-
  46. k[t, я]^2/2) EllipticK[k[t, я]^2]-EllipticE[k[t, я]^2]);
  47.  
  48. az[t_, я_] := If[z0==0 && vz==0, 0,
  49. 2G ρ z[t]/Abs[z[t]](-π+(R[t]^2+я^2+2r[t]я)^(-1/2)((R[t]-я)Sqrt[(R[t]-r[t])/(R[t]+
  50. r[t])] EllipticPi[2r[t]/(R[t]+r[t]), k[t, я]^2]+(R[t]+я)Sqrt[(R[t]+r[t])/(R[t]-
  51. r[t])] EllipticPi[-2r[t]/(R[t]-r[t]), k[t, я]^2]))]; (* Scheibenfeld *)
  52.  
  53. g[χ_] := -G Min[m, Ḿ] χ[t]/Sqrt[(x[t]^2+y[t]^2+z[t]^2)^3]; (* Kugelfeld *)
  54.  
  55. V[t_]:=2G ρ(z[t](π/2+π/2 Sign[я-r[t]])-Sqrt[R[t]^2+я^2+2r[t]я] EllipticE[k[t]^2]-(я^2-
  56. r[t]^2)/Sqrt[R[t]^2+я^2+2r[t] я] EllipticK[k[t]^2]-(я-r[t])/(я+
  57. r[t]) (z[t]^2)/Sqrt[R[t]^2+я^2+2r[t] я]EllipticPi[α[t]^2,k[t]^2]); (* Potential *)
  58.  
  59. Z0=z0; If[N[Z0]==0.0&&N[vz]=!=0.0, (z0=1/1*^6;), (z0=Z0;)];
  60.  
  61. sol=NDSolve[{ (* Differentialgleichung *)
  62.  
  63. x''[t]==ax[t, я2]-If[я1==0, 0, ax[t, я1]]+g[x], (* Beschleunigung x *)
  64. y''[t]==ay[t, я2]-If[я1==0, 0, ay[t, я1]]+g[y], (* Beschleunigung y *)
  65. z''[t]==az[t, я2]-If[я1==0, 0, az[t, я1]]+g[z], (* Beschleunigung z *)
  66.  
  67. x'[0]==vx,
  68. y'[0]==vy,
  69. z'[0]==vz,
  70.  
  71. x[0]==x0,
  72. y[0]==y0,
  73. z[0]==z0},
  74.  
  75. {x, y, z}, {t, 0, T},
  76.  
  77. WorkingPrecision-> wp,
  78. MaxSteps-> Infinity,
  79. Method-> mta,
  80. InterpolationOrder-> All,
  81. StepMonitor :> (laststep=plunge; plunge=t;
  82. stepsize=plunge-laststep;), Method->{"EventLocator",
  83. "Event" :> (If[stepsize<1*^-4, 0, 1])}];
  84.  
  85. X[t_]:=x[t]/.sol[[1]];
  86. Y[t_]:=y[t]/.sol[[1]];
  87. Z[t_]:=z[t]/.sol[[1]];
  88.  
  89. XYZ[t_]:=Sqrt[X[t]^2+Y[t]^2+Z[t]^2];
  90.  
  91. Xyz[{x_, y_, z_}, α_] := {x Cos[α]-y Sin[α], x Sin[α]+y Cos[α], z}; (* Rotationsmatrix *)
  92. xYz[{x_, y_, z_}, β_] := {x Cos[β]+z Sin[β], y, z Cos[β]-x Sin[β]};
  93. xyZ[{x_, y_, z_}, ψ_] := {x, y Cos[ψ]-z Sin[ψ], y Sin[ψ]+z Cos[ψ]};
  94.  
  95. φ[t_] := ArcTan[Y[t], X[t]]; (* Breitengrad *)
  96. θ[t_] := ArcCos[Z[t]/XYZ[t]]; (* Längengrad *)
  97.  
  98. cd = If[Ḿ<M, 2/5, 4/5]; ck = If[Ḿ<M, 4/5, 2/5]; (* Farbfunktion *)
  99.  
  100. annulus3dF[color_: GrayLevel[cd], o : OptionsPattern[]]:= (* Scheibe *)
  101. Graphics3D[{Glow[color], Opacity[0.8], EdgeForm[None], Black,
  102. ChartElementData["CylindricalSector3D"][{##}, 1]}, o,
  103. Boxed->False] &;
  104.  
  105. dp= Style[\!\(\*SuperscriptBox[\(Y\),\(Y\)]\), White]; n0[z_] := Chop[Re[N[Simplify[z]]]];
  106. s[text_]:=Style[text, FontFamily->"Consolas", FontSize->11]; (* Anzeigestil *)
  107.  
  108. PR = 40; (* Plot Range *)
  109. d1 = 50; (* Schweiflänge *)
  110. Plp = Max[100, Round[tp/2]]; (* Kurven Details *)
  111.  
  112. Mrec = 5000; (* Rekursionslimit *)
  113. mrec = 15; (* Parametric Plot Subdivisionen *)
  114. imgsize = 380; (* Bildgröße *)
  115.  
  116. display[tp_] := Grid[{{ (* numerisches Display *)
  117. Grid[{
  118. {s[" "], " ", s[" "], s[dp]},
  119. {s[" t"], " = ", s[n0[tp]], s[dp]},
  120. {s[" R"], " = ", s[n0[XYZ[tp]]], s[dp]},
  121. {s[" θ"], " = ", s[n0[θ[tp]]], s[dp]},
  122. {s[" φ"], " = ", s[n0[φ[tp]]], s[dp]},
  123. {s[" x"], " = ", s[n0[X[tp]]], s[dp]},
  124. {s[" y"], " = ", s[n0[Y[tp]]], s[dp]},
  125. {s[" z"], " = ", s[n0[Z[tp]]], s[dp]}
  126. }, Alignment-> Left, Spacings-> {0, 0}],
  127. Grid[{
  128. {s[" "], " ", s[" "], s[dp]},
  129. {s[" vt"], " = ", s[n0[Sqrt[X'[tp]^2+Y'[tp]^2+Z'[tp]^2]]], s[dp]},
  130. {s[" vR"], " = ", s[n0[XYZ'[tp]]], s[dp]},
  131. {s[" vθ"], " = ", s[n0[XYZ[tp] θ'[tp]]], s[dp]},
  132. {s[" vφ"], " = ", s[n0[XYZ[tp] φ'[tp]]], s[dp]},
  133. {s[" vx"], " = ", s[n0[X'[tp]]], s[dp]},
  134. {s[" vy"], " = ", s[n0[Y'[tp]]], s[dp]},
  135. {s[" vz"], " = ", s[n0[Z'[tp]]], s[dp]}
  136. }, Alignment-> Left, Spacings-> {0, 0}],
  137. Grid[{
  138. {s[" "], " ", s[" "], s[dp]},
  139. {s[" at"], " = ", s[n0[Sqrt[X''[tp]^2+Y''[tp]^2+Z''[tp]^2]]], s[dp]},
  140. {s[" aR"], " = ", s[n0[XYZ''[tp]]], s[dp]},
  141. {s[" aθ"], " = ", s[n0[XYZ[tp] θ''[tp]]], s[dp]},
  142. {s[" aφ"], " = ", s[n0[XYZ[tp] φ''[tp]]], s[dp]},
  143. {s[" ax"], " = ", s[n0[X''[tp]]], s[dp]},
  144. {s[" ay"], " = ", s[n0[Y''[tp]]], s[dp]},
  145. {s[" az"], " = ", s[n0[Z''[tp]]], s[dp]}
  146. }, Alignment-> Left, Spacings-> {0, 0}],
  147. Grid[{
  148. {s[" "], " ", s[" "], s[dp]},
  149. {s[" M"], " = ", s[n0[M]], s[dp]},
  150. {s[" Ḿ"], " = ", s[n0[Ḿ]], s[dp]},
  151. {s[" я1"], " = ", s[n0[я1]], s[dp]},
  152. {s[" я2"], " = ", s[n0[я2]], s[dp]},
  153. {s[" я3"], " = ", s[n0[я3]], s[dp]},
  154. {},
  155. {s[" "], " ", s["yukterez.net"], s[dp]}
  156. }, Alignment-> Left, Spacings-> {0, 0}]}
  157. }, Alignment-> Left, Spacings-> {0, 0}];
  158.  
  159. plot1b[{xx_, yy_, zz_, tp_}] := (* Animation *)
  160. Show[
  161.  
  162. Graphics3D[{Glow[GrayLevel[ck]], Black, Opacity[1], Sphere[{0, 0, 0}, я3]},
  163. ImageSize-> imgsize,
  164. PlotRange-> {
  165. {-(2 Sign[Abs[xx]]+1) PR, +(2 Sign[Abs[xx]]+1) PR},
  166. {-(2 Sign[Abs[yy]]+1) PR, +(2 Sign[Abs[yy]]+1) PR},
  167. {-(2 Sign[Abs[zz]]+1) PR, +(2 Sign[Abs[zz]]+1) PR}
  168. },
  169. SphericalRegion->False,
  170. ImagePadding-> 1], (* Kugel *)
  171.  
  172. annulus3dF[][{0, 2 Pi}, {я1, я2}, {-PR/150, PR/150}], (* Scheibe *)
  173.  
  174. Graphics3D[{ (* Testpartikel *)
  175. {PointSize[0.015], Red, Point[
  176. {X[tp], Y[tp], Z[tp]}]}},
  177. ImageSize-> imgsize,
  178. PlotRange-> PR,
  179. SphericalRegion->False,
  180. ImagePadding-> 1],
  181.  
  182. If[tp==0, {}, (* Schweif *)
  183. Block[{$RecursionLimit = Mrec},
  184. ParametricPlot3D[
  185. {X[tt], Y[tt], Z[tt]}, {tt, If[tp<0, Min[0, tp+d1], Max[0, tp-d1]], tp},
  186. PlotStyle-> {Thickness[PR/6000]},
  187. ColorFunction-> Function[{X, Y, Z, t},
  188. Hue[0, 1, 0.5, If[tp<0, Max[Min[(+tp+(-t+d1))/d1, 1], 0],
  189. Max[Min[(-tp+(t+d1))/d1, 1], 0]]]],
  190. ColorFunctionScaling-> False,
  191. PlotPoints-> Automatic,
  192. MaxRecursion-> mrec]]],
  193.  
  194. If[tp==0, {}, (* Trajektorie *)
  195. Block[{$RecursionLimit = Mrec},
  196. ParametricPlot3D[
  197. {X[tt], Y[tt], Z[tt]}, {tt, 0, If[tp<0,
  198. Min[-1*^-16, tp+d1/3], Max[1*^-16, tp-d1/3]]},
  199. PlotStyle-> {Thickness[PR/10000], Opacity[1], Darker@Gray},
  200. PlotPoints-> Plp,
  201. MaxRecursion-> mrec]]],
  202.  
  203. ViewPoint-> {xx, yy, zz}];
  204.  
  205. Quiet[Do[
  206. Print[Rasterize[Grid[{{
  207. Grid[{{
  208. plot1b[{0, -Infinity, 0, tp}],
  209. plot1b[{0, 0, +Infinity, tp}]
  210. }}, Alignment->Left]},
  211. {display[tp]}},
  212. Alignment->Left]]],
  213. {tp, 0, plunge, plunge/2}]]
  214.  
  215. (* Export als HTML Dokument *)
  216. (* Export["Y:\\export\\dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
  217. (* Export direkt als Bildsequenz *)
  218. (* Quiet[Do[Export["Y:\\export\\dateiname" <> ToString[tp] <> ".png", Rasterize[... *)
RAW Paste Data