Yukterez

Disk with Density Function, Orbit Simulator

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