Yukterez

Gravity & Charge, 3 Body Simulator

Feb 16th, 2019 (edited)
52
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* ||| Mathematica Syntax || yukterez.net || 3 Body Newtonian Mass & Charge Simulator ||| *)
  3. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  4.  
  5. ClearAll["Global`*"]; ClearAll["Local`*"];
  6. Needs["DifferentialEquations`NDSolveProblems`"];
  7. Needs["DifferentialEquations`NDSolveUtilities`"];
  8.  
  9. Amp = 1; kg = 1; m = 1; sek = 1; km = 1000 m; (* SI Einheiten *)
  10.  
  11. mt1 = {"StiffnessSwitching", Method-> {"ExplicitRungeKutta", Automatic}};
  12. mt2 = {"ImplicitRungeKutta", "DifferenceOrder"-> 20};
  13. mt3 = {"EquationSimplification"-> "Residual"};
  14. mt0 = Automatic;
  15. mta = mt2;
  16. wp = MachinePrecision;
  17.  
  18. (* Plot Optionen *)
  19.  
  20. Tmax = 10000 sek;
  21. tMax = Min[Tmax, plunge];
  22. trail = 3000 sek;
  23. point = 0.015;
  24. thk = 0.004;
  25. plotrange = 1 m {{-0.2, +1.2}, {-0.6, +0.6}, {-0.2, +1.2}};
  26. viewpoint = {0, Infinity, 0};
  27. imagesize = 430;
  28. startpos = 0;
  29.  
  30. (* Konstanten *)
  31.  
  32. G = 667384/10^16 m^3/kg/sek^2;
  33. Λ = 11056*^-56/m^2;
  34. ε0 = 8854187817*^-21 Amp^2 sek^4/kg/m^3;
  35. c = 299792458 m/sek;
  36. Au = 149597870700 m;
  37. dy = 24*3600 sek;
  38. yr = 36525*dy/100;
  39.  
  40. (* Körper 1 *)
  41.  
  42. m1 = 1000 kg;
  43. q1 = 0;
  44.  
  45. x1x = 1/2 m;
  46. y1y = 0 m;
  47. z1z = Sqrt[3]/2 m;
  48.  
  49. v1x = 0 m/sek;
  50. v1y = 0 m/sek;
  51. v1z = 0 m/sek;
  52.  
  53. (* Körper 2 *)
  54.  
  55. m2 = 2000/3 kg;
  56. q2 = 0 Amp sek;
  57.  
  58. x2x = 1/10 m;
  59. y2y = 0 m;
  60. z2z = 0 m;
  61.  
  62. v2x = 0 m/sek;
  63. v2y = 0 m/sek;
  64. v2z = 0 m/sek;
  65.  
  66. (* Körper 3 *)
  67.  
  68. m3 = 500 kg;
  69. q3 = 0 Amp sek;
  70.  
  71. x3x = 4/5 m;
  72. y3y = 0 m;
  73. z3z = 1/5 m;
  74.  
  75. v3x = 0 m/sek;
  76. v3y = 0 m/sek;
  77. v3z = 0 m/sek;
  78.  
  79. (* Differentialgleichung *)
  80.  
  81. nds=NDSolve[{
  82.  
  83. x1'[t] == vx1[t], y1'[t] == vy1[t], z1'[t] == vz1[t],
  84. x2'[t] == vx2[t], y2'[t] == vy2[t], z2'[t] == vz2[t],
  85. x3'[t] == vx3[t], y3'[t] == vy3[t], z3'[t] == vz3[t],
  86.  
  87. vx1'[t] ==
  88. (G m2 (x2[t]-x1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  89. (G m3 (x3[t]-x1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  90. If[q1 == 0, 0,
  91. (-q1*q2/(4Pi ε0 )/m1 (x2[t]-x1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  92. (-q1*q3/(4Pi ε0 )/m1 (x3[t]-x1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]]+
  93. Λ*c^2*x1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  94.  
  95. vy1'[t] ==
  96. (G m2 (y2[t]-y1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  97. (G m3 (y3[t]-y1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  98. If[q1 == 0, 0,
  99. (-q1*q2/(4Pi ε0 )/m1 (y2[t]-y1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  100. (-q1*q3/(4Pi ε0 )/m1 (y3[t]-y1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]]+
  101. Λ*c^2*y1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  102.  
  103. vz1'[t] ==
  104. (G m2 (z2[t]-z1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  105. (G m3 (z3[t]-z1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  106. If[q1 == 0, 0,
  107. (-q1*q2/(4Pi ε0 )/m1 (z2[t]-z1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  108. (-q1*q3/(4Pi ε0 )/m1 (z3[t]-z1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]]+
  109. Λ*c^2*z1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  110.  
  111. vx2'[t] ==
  112. (G m1 (x1[t]-x2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  113. (G m3 (x3[t]-x2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  114. If[q2 == 0, 0,
  115. (-q2*q1/(4Pi ε0 )/m2 (x1[t]-x2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  116. (-q2*q3/(4Pi ε0 )/m2 (x3[t]-x2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]]+
  117. Λ*c^2*x2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  118.  
  119. vy2'[t] ==
  120. (G m1 (y1[t]-y2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  121. (G m3 (y3[t]-y2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  122. If[q2 == 0, 0,
  123. (-q2*q1/(4Pi ε0 )/m2 (y1[t]-y2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  124. (-q2*q3/(4Pi ε0 )/m2 (y3[t]-y2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]]+
  125. Λ*c^2*y2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  126.  
  127. vz2'[t] ==
  128. (G m1 (z1[t]-z2[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  129. (G m3 (z3[t]-z2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  130. If[q2 == 0, 0,
  131. (-q2*q1/(4Pi ε0 )/m2 (z1[t]-z2[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  132. (-q2*q3/(4Pi ε0 )/m2 (z3[t]-z2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]]+
  133. Λ*c^2*z2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  134.  
  135. vx3'[t] ==
  136. (G m1 (x1[t]-x3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  137. (G m2 (x2[t]-x3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  138. If[q3 == 0, 0,
  139. (-q3*q1/(4Pi ε0 )/m3 (x1[t]-x3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  140. (-q3*q2/(4Pi ε0 )/m3 (x2[t]-x3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]]+
  141. Λ*c^2*x3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  142.  
  143. vy3'[t] ==
  144. (G m1 (y1[t]-y3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  145. (G m2 (y2[t]-y3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  146. If[q3 == 0, 0,
  147. (-q3*q1/(4Pi ε0 )/m3 (y1[t]-y3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  148. (-q3*q2/(4Pi ε0 )/m3 (y2[t]-y3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]]+
  149. Λ*c^2*y3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  150.  
  151. vz3'[t] ==
  152. (G m1 (z1[t]-z3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  153. (G m2 (z2[t]-z3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  154. If[q3 == 0, 0,
  155. (-q3*q1/(4Pi ε0 )/m3 (z1[t]-z3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  156. (-q3*q2/(4Pi ε0 )/m3 (z2[t]-z3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]]+
  157. Λ*c^2*z3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  158.  
  159. x1[0] == x1x, y1[0] == y1y, z1[0] == z1z,
  160. x2[0] == x2x, y2[0] == y2y, z2[0] == z2z,
  161. x3[0] == x3x, y3[0] == y3y, z3[0] == z3z,
  162.  
  163. vx1[0] == v1x, vy1[0] == v1y, vz1[0] == v1z,
  164. vx2[0] == v2x, vy2[0] == v2y, vz2[0] == v2z,
  165. vx3[0] == v3x, vy3[0] == v3y, vz3[0] == v3z},
  166.  
  167. {x1, x2, x3, y1, y2, y3, z1, z2, z3,
  168. vx1, vx2, vx3, vy1, vy2, vy3, vz1, vz2, vz3},
  169.  
  170. {t, 0, Tmax},
  171.  
  172. WorkingPrecision-> wp,
  173. MaxSteps-> Infinity,
  174. Method-> mta,
  175. InterpolationOrder-> All,
  176. StepMonitor :> (laststep=plunge; plunge=t;
  177. stepsize=plunge-laststep;), Method->{"EventLocator",
  178. "Event" :> (If[stepsize<1*^-4, 0, 1])}];
  179.  
  180. (* Position, Geschwindigkeit *)
  181.  
  182. f2p[t_]={{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}, {x3[t], y3[t], z3[t]}}/.nds[[1]];
  183. f2v[t_]={{vx1[t], vy1[t], vz1[t]}, {vx2[t], vy2[t], vz2[t]}, {vx3[t], vy3[t], vz3[t]}}/.nds[[1]];
  184. swp[t_]=(m1 Evaluate[f2p[t][[1]]]+m2 Evaluate[f2p[t][[2]]]+m3 Evaluate[f2p[t][[3]]])/(m1+m2+m3);
  185.  
  186. (* Formatierung *)
  187.  
  188. s[text_]=Style[text, FontSize->11];
  189. sw[text_]=Style[text, White, FontSize->11];
  190. colorfunc[n_]=Function[{x, y, z, t},
  191. Hue[0, n, 0.5,
  192. If[Tmax<0, Max[Min[(+T+(-t+trail))/trail, 1], 0],
  193. Max[Min[(-T+(t+trail))/trail, 1], 0]]]];
  194.  
  195. (* Animation *)
  196.  
  197. Do[Print[Rasterize[
  198. Grid[{{
  199. Show[
  200.  
  201. If[T == 0, {},
  202.  
  203. ParametricPlot3D[Evaluate[f2p[t]],
  204. {t, Max[0, T-trail], T},
  205.  
  206. PlotStyle->{
  207. {Thickness[thk], Red},
  208. {Thickness[thk], Blue},
  209. {Thickness[thk], Green}},
  210.  
  211. PlotRange->plotrange, AspectRatio->1, MaxRecursion->15, Axes->True, ImageSize->imagesize]],
  212.  
  213. Graphics3D[
  214. If[startpos==1, {
  215. {PointSize[2point/3], Lighter[Red], Point[{x1x, y1y, z1z}]},
  216. {PointSize[2point/3], Lighter[Blue], Point[{x2x, y2y, z2z}]},
  217. {PointSize[2point/3], Lighter[Green], Point[{x3x, y3y, z3z}]}
  218. }, {}],
  219.  
  220. PlotRange->plotrange, AspectRatio->1, Axes->True, ImageSize->imagesize],
  221.  
  222. Graphics3D[{PointSize[point], Red, Point[Evaluate[f2p[T]][[1]]]}],
  223. Graphics3D[{PointSize[point], Blue, Point[Evaluate[f2p[T]][[2]]]}],
  224. Graphics3D[{PointSize[point], Green, Point[Evaluate[f2p[T]][[3]]]}],
  225.  
  226. ViewPoint->viewpoint]},
  227.  
  228. { },
  229. {s["t"->N[T]], sw[1/2]},
  230. { },
  231. {s["p1{x,y,z}"-> Evaluate[f2p[T][[1]]]], sw[1/2]},
  232. {s["v1{x,y,z}"-> Evaluate[f2v[T][[1]]]], sw[1/2]},
  233. {s["v1{total}"->{Evaluate[Chop@Norm[f2v[T][[1]]]]}], sw[1/2]},
  234. { },
  235. {s["p2{x,y,z}"-> Evaluate[f2p[T][[2]]]], sw[1/2]},
  236. {s["v2{x,y,z}"-> Evaluate[f2v[T][[2]]]], sw[1/2]},
  237. {s["v2{total}"->{Evaluate[Chop@Norm[f2v[T][[2]]]]}], sw[1/2]},
  238. { },
  239. {s["p3{x,y,z}"-> Evaluate[f2p[T][[3]]]], sw[1/2]},
  240. {s["v3{x,y,z}"-> Evaluate[f2v[T][[3]]]], sw[1/2]},
  241. {s["v3{total}"->{Evaluate[Chop@Norm[f2v[T][[3]]]]}], sw[1/2]},
  242. { },
  243. {s["ps{x,y,z}"-> swp[T]], sw[1/2]},
  244. {s["vs{x,y,z}"-> swp'[T]], sw[1/2]},
  245. {s["vs{total}"->{Chop@Norm[swp'[T]]}], sw[1/2]}
  246. }, Alignment->Left]]],
  247.  
  248. (* Zeitregler *)
  249.  
  250. {T, 0, tMax, tMax/5}]
  251.  
  252. (* Export als HTML Dokument *)
  253. (* Export["dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
  254. (* Export direkt als Bildsequenz *)
  255. (* ParallelDo[Export["dateiname" <> ToString[T] <> ".png", Rasterize[...] ], {T, 0, 10, 5}] *)
  256.  
  257.  
  258.  
RAW Paste Data