Yukterez

Gravity & Charge, 5 Body Simulator

Feb 14th, 2019 (edited)
139
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* ||| Mathematica Syntax || yukterez.net || 5 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 = 2π sek;
  21. tMax = Min[Tmax, plunge];
  22. trail = π/3 sek;
  23. point = 0.015;
  24. thk = 0.004;
  25. plotrange = 1.2 m {{-1, +1}, {-1, +1}, {-1, +1}};
  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 = 1 m^3/sek^2/G;
  43. q1 = 0 Amp sek;
  44.  
  45. x1x = 0.3673027525587564 m;
  46. y1y = 0 m;
  47. z1z = 0 m;
  48.  
  49. v1x = 0 m/sek;
  50. v1y = 0 m/sek;
  51. v1z = 2.134808238913507 m/sek;
  52.  
  53. (* Körper 2 *)
  54.  
  55. m2 = 1 m^3/sek^2/G;
  56. q2 = 0 Amp sek;
  57.  
  58. x2x = 0.6760231721222919 m;
  59. y2y = 0 m;
  60. z2z = -0.1852023639618821 m;
  61.  
  62. v2x = -2.138115663581648 m/sek;
  63. v2y = 0 m/sek;
  64. v2z = -0.7311482987140561 m/sek;
  65.  
  66. (* Körper 3 *)
  67.  
  68. m3 = 1 m^3/sek^2/G;
  69. q3 = 0 Amp sek;
  70.  
  71. x3x = -0.859674548416703 m;
  72. y3y = 0 m;
  73. z3z = -0.3897351750882059 m;
  74.  
  75. v3x = 0.5365854530502069 m/sek;
  76. v3y = 0 m/sek;
  77. v3z = -0.3362558207426968 m/sek;
  78.  
  79. (* Körper 4 *)
  80.  
  81. m4 = 1 m^3/sek^2/G;
  82. q4 = 0 Amp sek;
  83.  
  84. x4x = -0.8596745484016705 m;
  85. y4y = 0 m;
  86. z4z = 0.3897351750882059 m;
  87.  
  88. v4x = -0.5365854530502073 m/sek;
  89. v4y = 0 m/sek;
  90. v4z = -0.3362558207426963 m/sek;
  91.  
  92. (* Körper 5 *)
  93.  
  94. m5 = 1 m^3/sek^2/G;
  95. q5 = 0 Amp sek;
  96.  
  97. x5x = 0.6760231721222916 m;
  98. y5y = 0 m;
  99. z5z = 0.1852023639618821 m;
  100.  
  101. v5x = 2.138115663581649 m/sek;
  102. v5y = 0 m/sek;
  103. v5z = -0.7311482987140556 m/sek;
  104.  
  105. (* Differentialgleichung *)
  106.  
  107. nds=NDSolve[{
  108.  
  109. x1'[t] == vx1[t], y1'[t] == vy1[t], z1'[t] == vz1[t],
  110. x2'[t] == vx2[t], y2'[t] == vy2[t], z2'[t] == vz2[t],
  111. x3'[t] == vx3[t], y3'[t] == vy3[t], z3'[t] == vz3[t],
  112. x4'[t] == vx4[t], y4'[t] == vy4[t], z4'[t] == vz4[t],
  113. x5'[t] == vx5[t], y5'[t] == vy5[t], z5'[t] == vz5[t],
  114.  
  115. vx1'[t] ==
  116. (G m2 (x2[t]-x1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  117. (G m3 (x3[t]-x1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  118. (G m4 (x4[t]-x1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  119. (G m5 (x5[t]-x1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]+
  120. If[q1 == 0, 0,
  121. (-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]+
  122. (-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]+
  123. (-q1*q4/(4Pi ε0 )/m1 (x4[t]-x1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  124. (-q1*q5/(4Pi ε0 )/m1 (x5[t]-x1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]]+
  125. Λ*c^2*x1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  126.  
  127. vy1'[t] ==
  128. (G m2 (y2[t]-y1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  129. (G m3 (y3[t]-y1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  130. (G m4 (y4[t]-y1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  131. (G m5 (y5[t]-y1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]+
  132. If[q1 == 0, 0,
  133. (-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]+
  134. (-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]+
  135. (-q1*q4/(4Pi ε0 )/m1 (y4[t]-y1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  136. (-q1*q5/(4Pi ε0 )/m1 (y5[t]-y1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]]+
  137. Λ*c^2*y1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  138.  
  139. vz1'[t] ==
  140. (G m2 (z2[t]-z1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  141. (G m3 (z3[t]-z1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  142. (G m4 (z4[t]-z1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  143. (G m5 (z5[t]-z1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]+
  144. If[q1 == 0, 0,
  145. (-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]+
  146. (-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]+
  147. (-q1*q4/(4Pi ε0 )/m1 (z4[t]-z1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  148. (-q1*q5/(4Pi ε0 )/m1 (z5[t]-z1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]]+
  149. Λ*c^2*z1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  150.  
  151. vx2'[t] ==
  152. (G m1 (x1[t]-x2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  153. (G m3 (x3[t]-x2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  154. (G m4 (x4[t]-x2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  155. (G m5 (x5[t]-x2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]+
  156. If[q2 == 0, 0,
  157. (-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]+
  158. (-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]+
  159. (-q2*q4/(4Pi ε0 )/m2 (x4[t]-x2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  160. (-q2*q5/(4Pi ε0 )/m2 (x5[t]-x2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]]+
  161. Λ*c^2*x2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  162.  
  163. vy2'[t] ==
  164. (G m1 (y1[t]-y2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  165. (G m3 (y3[t]-y2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  166. (G m4 (y4[t]-y2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  167. (G m5 (y5[t]-y2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]+
  168. If[q2 == 0, 0,
  169. (-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]+
  170. (-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]+
  171. (-q2*q4/(4Pi ε0 )/m2 (y4[t]-y2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  172. (-q2*q5/(4Pi ε0 )/m2 (y5[t]-y2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]]+
  173. Λ*c^2*y2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  174.  
  175. vz2'[t] ==
  176. (G m1 (z1[t]-z2[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  177. (G m3 (z3[t]-z2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  178. (G m4 (z4[t]-z2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  179. (G m5 (z5[t]-z2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]+
  180. If[q2 == 0, 0,
  181. (-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]+
  182. (-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]+
  183. (-q2*q4/(4Pi ε0 )/m2 (z4[t]-z2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  184. (-q2*q5/(4Pi ε0 )/m2 (z5[t]-z2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]]+
  185. Λ*c^2*z2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  186.  
  187. vx3'[t] ==
  188. (G m1 (x1[t]-x3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  189. (G m2 (x2[t]-x3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  190. (G m4 (x4[t]-x3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  191. (G m5 (x5[t]-x3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]+
  192. If[q3 == 0, 0,
  193. (-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]+
  194. (-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]+
  195. (-q3*q4/(4Pi ε0 )/m3 (x4[t]-x3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  196. (-q3*q5/(4Pi ε0 )/m3 (x5[t]-x3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]]+
  197. Λ*c^2*x3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  198.  
  199. vy3'[t] ==
  200. (G m1 (y1[t]-y3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  201. (G m2 (y2[t]-y3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  202. (G m4 (y4[t]-y3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  203. (G m5 (y5[t]-y3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]+
  204. If[q3 == 0, 0,
  205. (-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]+
  206. (-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]+
  207. (-q3*q4/(4Pi ε0 )/m3 (y4[t]-y3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  208. (-q3*q5/(4Pi ε0 )/m3 (y5[t]-y3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]]+
  209. Λ*c^2*y3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  210.  
  211. vz3'[t] ==
  212. (G m1 (z1[t]-z3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  213. (G m2 (z2[t]-z3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  214. (G m4 (z4[t]-z3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  215. (G m5 (z5[t]-z3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]+
  216. If[q3 == 0, 0,
  217. (-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]+
  218. (-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]+
  219. (-q3*q4/(4Pi ε0 )/m3 (z4[t]-z3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  220. (-q3*q5/(4Pi ε0 )/m3 (z5[t]-z3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]]+
  221. Λ*c^2*z3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  222.  
  223. vx4'[t] ==
  224. (G m1 (x1[t]-x4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  225. (G m2 (x2[t]-x4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  226. (G m3 (x3[t]-x4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  227. (G m5 (x5[t]-x4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]+
  228. If[q4 == 0, 0,
  229. (-q4*q1/(4Pi ε0 )/m4 (x1[t]-x4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  230. (-q4*q2/(4Pi ε0 )/m4 (x2[t]-x4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  231. (-q4*q3/(4Pi ε0 )/m4 (x3[t]-x4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  232. (-q4*q5/(4Pi ε0 )/m4 (x5[t]-x4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]]+
  233. Λ*c^2*x4[t]^2/Sqrt[x4[t]^2+y4[t]^2+z4[t]^2],
  234.  
  235. vy4'[t] ==
  236. (G m1 (y1[t]-y4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  237. (G m2 (y2[t]-y4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  238. (G m3 (y3[t]-y4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  239. (G m5 (y5[t]-y4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]+
  240. If[q4 == 0, 0,
  241. (-q4*q1/(4Pi ε0 )/m4 (y1[t]-y4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  242. (-q4*q2/(4Pi ε0 )/m4 (y2[t]-y4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  243. (-q4*q3/(4Pi ε0 )/m4 (y3[t]-y4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  244. (-q4*q5/(4Pi ε0 )/m4 (y5[t]-y4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]]+
  245. Λ*c^2*y4[t]^2/Sqrt[x4[t]^2+y4[t]^2+z4[t]^2],
  246.  
  247. vz4'[t] ==
  248. (G m1 (z1[t]-z4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  249. (G m2 (z2[t]-z4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  250. (G m3 (z3[t]-z4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  251. (G m5 (z5[t]-z4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]+
  252. If[q4 == 0, 0,
  253. (-q4*q1/(4Pi ε0 )/m4 (z1[t]-z4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  254. (-q4*q2/(4Pi ε0 )/m4 (z2[t]-z4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  255. (-q4*q3/(4Pi ε0 )/m4 (z3[t]-z4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  256. (-q4*q5/(4Pi ε0 )/m4 (z5[t]-z4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]]+
  257. Λ*c^2*z4[t]^2/Sqrt[x4[t]^2+y4[t]^2+z4[t]^2],
  258.  
  259. vx5'[t] ==
  260. (G m1 (x1[t]-x5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  261. (G m2 (x2[t]-x5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  262. (G m3 (x3[t]-x5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  263. (G m4 (x4[t]-x5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]+
  264. If[q5 == 0, 0,
  265. (-q5*q1/(4Pi ε0 )/m5 (x1[t]-x5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  266. (-q5*q2/(4Pi ε0 )/m5 (x2[t]-x5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  267. (-q5*q3/(4Pi ε0 )/m5 (x3[t]-x5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  268. (-q5*q4/(4Pi ε0 )/m5 (x4[t]-x5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]]+
  269. Λ*c^2*x5[t]^2/Sqrt[x5[t]^2+y5[t]^2+z5[t]^2],
  270.  
  271. vy5'[t] ==
  272. (G m1 (y1[t]-y5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  273. (G m2 (y2[t]-y5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  274. (G m3 (y3[t]-y5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  275. (G m4 (y4[t]-y5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]+
  276. If[q5 == 0, 0,
  277. (-q5*q1/(4Pi ε0 )/m5 (y1[t]-y5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  278. (-q5*q2/(4Pi ε0 )/m5 (y2[t]-y5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  279. (-q5*q3/(4Pi ε0 )/m5 (y3[t]-y5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  280. (-q5*q4/(4Pi ε0 )/m5 (y4[t]-y5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]]+
  281. Λ*c^2*y5[t]^2/Sqrt[x5[t]^2+y5[t]^2+z5[t]^2],
  282.  
  283. vz5'[t] ==
  284. (G m1 (z1[t]-z5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  285. (G m2 (z2[t]-z5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  286. (G m3 (z3[t]-z5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  287. (G m4 (z4[t]-z5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]+
  288. If[q5 == 0, 0,
  289. (-q5*q1/(4Pi ε0 )/m5 (z1[t]-z5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  290. (-q5*q2/(4Pi ε0 )/m5 (z2[t]-z5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  291. (-q5*q3/(4Pi ε0 )/m5 (z3[t]-z5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  292. (-q5*q4/(4Pi ε0 )/m5 (z4[t]-z5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]]+
  293. Λ*c^2*z5[t]^2/Sqrt[x5[t]^2+y5[t]^2+z5[t]^2],
  294.  
  295. x1[0] == x1x, y1[0] == y1y, z1[0] == z1z,
  296. x2[0] == x2x, y2[0] == y2y, z2[0] == z2z,
  297. x3[0] == x3x, y3[0] == y3y, z3[0] == z3z,
  298. x4[0] == x4x, y4[0] == y4y, z4[0] == z4z,
  299. x5[0] == x5x, y5[0] == y5y, z5[0] == z5z,
  300.  
  301. vx1[0] == v1x, vy1[0] == v1y, vz1[0] == v1z,
  302. vx2[0] == v2x, vy2[0] == v2y, vz2[0] == v2z,
  303. vx3[0] == v3x, vy3[0] == v3y, vz3[0] == v3z,
  304. vx4[0] == v4x, vy4[0] == v4y, vz4[0] == v4z,
  305. vx5[0] == v5x, vy5[0] == v5y, vz5[0] == v5z},
  306.  
  307. {x1, x2, x3, x4, x5, y1, y2, y3, y4, y5, z1, z2, z3, z4, z5,
  308. vx1, vx2, vx3, vx4, vx5, vy1, vy2, vy3, vy4, vy5, vz1, vz2, vz3, vz4, vz5},
  309.  
  310. {t, 0, Tmax},
  311.  
  312. WorkingPrecision-> wp,
  313. MaxSteps-> Infinity,
  314. Method-> mta,
  315. InterpolationOrder-> All,
  316. StepMonitor :> (laststep=plunge; plunge=t;
  317. stepsize=plunge-laststep;), Method->{"EventLocator",
  318. "Event" :> (If[stepsize<1*^-4, 0, 1])}];
  319.  
  320. (* Position, Geschwindigkeit *)
  321.  
  322. f2p[t_]={{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}, {x3[t], y3[t], z3[t]}, {x4[t], y4[t], z4[t]}, {x5[t], y5[t], z5[t]}}/.nds[[1]];
  323. f2v[t_]={{vx1[t], vy1[t], vz1[t]}, {vx2[t], vy2[t], vz2[t]}, {vx3[t], vy3[t], vz3[t]}, {vx4[t], vy4[t], vz4[t]}, {vx5[t], vy5[t], vz5[t]}}/.nds[[1]];
  324. swp[t_]=(m1 Evaluate[f2p[t][[1]]]+m2 Evaluate[f2p[t][[2]]]+m3 Evaluate[f2p[t][[3]]]+m4 Evaluate[f2p[t][[4]]]+m5 Evaluate[f2p[t][[5]]])/(m1+m2+m3+m4+m5);
  325.  
  326. (* Formatierung *)
  327.  
  328. s[text_]=Style[text, FontSize->11];
  329. sw[text_]=Style[text, White, FontSize->11];
  330. colorfunc[n_]=Function[{x, y, z, t},
  331. Hue[0, n, 0.5,
  332. If[Tmax<0, Max[Min[(+T+(-t+trail))/trail, 1], 0],
  333. Max[Min[(-T+(t+trail))/trail, 1], 0]]]];
  334.  
  335. (* Animation *)
  336.  
  337. Do[Print[Rasterize[
  338. Grid[{{
  339. Show[
  340.  
  341. If[T == 0, {},
  342.  
  343. ParametricPlot3D[Evaluate[f2p[t]],
  344. {t, Max[0, T-trail], T},
  345.  
  346. PlotStyle->{
  347. {Thickness[thk], Red},
  348. {Thickness[thk], Blue},
  349. {Thickness[thk], Green},
  350. {Thickness[thk], Magenta},
  351. {Thickness[thk], Cyan}},
  352.  
  353. PlotRange->plotrange, AspectRatio->1, MaxRecursion->15, Axes->True, ImageSize->imagesize]],
  354.  
  355. Graphics3D[
  356. If[startpos==1, {
  357. {PointSize[2point/3], Lighter[Red], Point[{x1x, y1y, z1z}]},
  358. {PointSize[2point/3], Lighter[Blue], Point[{x2x, y2y, z2z}]},
  359. {PointSize[2point/3], Lighter[Green], Point[{x3x, y3y, z3z}]},
  360. {PointSize[2point/3], Lighter[Magenta], Point[{x4x, y4y, z4z}]},
  361. {PointSize[2point/3], Lighter[Cyan], Point[{x5x, y5y, z5z}]}
  362. }, {}],
  363.  
  364. PlotRange->plotrange, AspectRatio->1, Axes->True, ImageSize->imagesize],
  365.  
  366. Graphics3D[{PointSize[point], Red, Point[Evaluate[f2p[T]][[1]]]}],
  367. Graphics3D[{PointSize[point], Blue, Point[Evaluate[f2p[T]][[2]]]}],
  368. Graphics3D[{PointSize[point], Green, Point[Evaluate[f2p[T]][[3]]]}],
  369. Graphics3D[{PointSize[point], Magenta, Point[Evaluate[f2p[T]][[4]]]}],
  370. Graphics3D[{PointSize[point], Cyan, Point[Evaluate[f2p[T]][[5]]]}],
  371.  
  372. ViewPoint->viewpoint]},
  373.  
  374. { },
  375. {s["t"->N[T]], sw[1/2]},
  376. { },
  377. {s["p1{x,y,z}"-> Evaluate[f2p[T][[1]]]], sw[1/2]},
  378. {s["v1{x,y,z}"-> Evaluate[f2v[T][[1]]]], sw[1/2]},
  379. {s["v1{total}"->{Evaluate[Chop@Norm[f2v[T][[1]]]]}], sw[1/2]},
  380. { },
  381. {s["p2{x,y,z}"-> Evaluate[f2p[T][[2]]]], sw[1/2]},
  382. {s["v2{x,y,z}"-> Evaluate[f2v[T][[2]]]], sw[1/2]},
  383. {s["v2{total}"->{Evaluate[Chop@Norm[f2v[T][[2]]]]}], sw[1/2]},
  384. { },
  385. {s["p3{x,y,z}"-> Evaluate[f2p[T][[3]]]], sw[1/2]},
  386. {s["v3{x,y,z}"-> Evaluate[f2v[T][[3]]]], sw[1/2]},
  387. {s["v3{total}"->{Evaluate[Chop@Norm[f2v[T][[3]]]]}], sw[1/2]},
  388. { },
  389. {s["p4{x,y,z}"-> Evaluate[f2p[T][[4]]]], sw[1/2]},
  390. {s["v4{x,y,z}"-> Evaluate[f2v[T][[4]]]], sw[1/2]},
  391. {s["v4{total}"->{Evaluate[Chop@Norm[f2v[T][[4]]]]}], sw[1/2]},
  392. { },
  393. {s["p5{x,y,z}"-> Evaluate[f2p[T][[5]]]], sw[1/2]},
  394. {s["v5{x,y,z}"-> Evaluate[f2v[T][[5]]]], sw[1/2]},
  395. {s["v5{total}"->{Evaluate[Chop@Norm[f2v[T][[5]]]]}], sw[1/2]},
  396. { },
  397. {s["ps{x,y,z}"-> swp[T]], sw[1/2]},
  398. {s["vs{x,y,z}"-> swp'[T]], sw[1/2]},
  399. {s["vs{total}"->{Chop@Norm[swp'[T]]}], sw[1/2]}
  400. }, Alignment->Left]]],
  401.  
  402. (* Zeitregler *)
  403.  
  404. {T, 0, tMax, tMax/5}]
  405.  
  406. (* Export als HTML Dokument *)
  407. (* Export["dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
  408. (* Export direkt als Bildsequenz *)
  409. (* ParallelDo[Export["dateiname" <> ToString[T] <> ".png", Rasterize[...] ], {T, 0, 10, 5}] *)
  410.  
  411.  
RAW Paste Data