Yukterez

Gravity & Charge, 9 Body Simulator

Feb 19th, 2019 (edited)
106
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
  2. (* ||| Mathematica Syntax || yukterez.net || 9 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 = 40 yr;
  21. tMax = Min[Tmax, plunge];
  22. trail = 10 yr;
  23. point = 0.015;
  24. thk = 0.004;
  25. plotrange = 50 Au {{-1, +1}, {-1, +1}, {-1, +1}};
  26. viewpoint = {40, 30, 20};
  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. (* Ephemeriden vom 19.02.2019, 0:00:00 TDB *)
  40. (* Sonne *)
  41.  
  42. m1 = +1.988435*^30 kg;
  43. q1 = +77 Amp sek;
  44.  
  45. x1x = -1.147196570503204*^-03 Au;
  46. y1y = +7.515074431920434*^-03 Au;
  47. z1z = -4.730273651193038*^-05 Au;
  48.  
  49. v1x = -8.107931162902937*^-06 Au/dy;
  50. v1y = +1.520849732928662*^-06 Au/dy;
  51. v1z = +2.095554598567427*^-07 Au/dy;
  52.  
  53. (* Merkur *)
  54.  
  55. m2 = +3.30104*^23 kg;
  56. q2 = +0 Amp sek;
  57.  
  58. x2x = +2.493682187528474*^-01 Au;
  59. y2y = +2.060848667278006*^-01 Au;
  60. z2z = -6.803162776737710*^-03 Au;
  61.  
  62. v2x = -2.301828852252654*^-02 Au/dy;
  63. v2y = +2.326003199133993*^-02 Au/dy;
  64. v2z = +4.011640539083395*^-03 Au/dy;
  65.  
  66. (* Venus *)
  67.  
  68. m3 = +4.86732*^24 kg;
  69. q3 = +0 Amp sek;
  70.  
  71. x3x = -5.604572600267276*^-01 Au;
  72. y3y = -4.500554270408416*^-01 Au;
  73. z3z = +2.595073246894732*^-02 Au;
  74.  
  75. v3x = +1.265689462094818*^-02 Au/dy;
  76. v3y = -1.574829638876520*^-02 Au/dy;
  77. v3z = -9.467652690844731*^-04 Au/dy;
  78.  
  79. (* Erde + Mond *)
  80.  
  81. m4 = +5.9721986*^24 kg+7.3459*^22 kg;
  82. q4 = +0 Amp sek;
  83.  
  84. x4x = -8.552072163834489*^-01 Au;
  85. y4y = +5.049715021822364*^-01 Au;
  86. z4z = -6.849877545851131*^-05 Au;
  87.  
  88. v4x = -8.942912568116291*^-03 Au/dy;
  89. v4y = -1.492365678503182*^-02 Au/dy;
  90. v4z = +2.741178622694643*^-07 Au/dy;
  91.  
  92. (* Mars *)
  93.  
  94. m5 = +6.41693*^23 kg;
  95. q5 = +0 Amp sek;
  96.  
  97. x5x = +5.580724605736193*^-01 Au;
  98. y5y = +1.416261572201534*^+00 Au;
  99. z5z = +1.574925082740965*^-02 Au;
  100.  
  101. v5x = -1.248544019487808*^-02 Au/dy;
  102. v5y = +6.355083417008326*^-03 Au/dy;
  103. v5z = +4.394992947386628*^-04 Au/dy;
  104.  
  105. (* Jupiter *)
  106.  
  107. m6 = +1.89813*^27 kg;
  108. q6 = +0 Amp sek;
  109.  
  110. x6x = -1.795821860926694*^+00 Au;
  111. y6y = -5.016469167174772*^+00 Au;
  112. z6z = +6.097587180308248*^-02 Au;
  113.  
  114. v6x = +7.014525824256318*^-03 Au/dy;
  115. v6y = -2.183010990796764*^-03 Au/dy;
  116. v6z = -1.478090774743338*^-04 Au/dy;
  117.  
  118. (* Saturn *)
  119.  
  120. m7 = +5.68319*^26 kg;
  121. q7 = +0 Amp sek;
  122.  
  123. x7x = +2.211165351380597*^+00 Au;
  124. y7y = -9.803846216723874*^+00 Au;
  125. z7z = +8.244475037063657*^-02 Au;
  126.  
  127. v7x = +5.133965065556525*^-03 Au/dy;
  128. v7y = +1.210333590471664*^-03 Au/dy;
  129. v7z = -2.255855621236429*^-04 Au/dy;
  130.  
  131. (* Uranus *)
  132.  
  133. m8 = +8.68103*^25 kg;
  134. q8 = +0 Amp sek;
  135.  
  136. x8x = +1.691367572961052*^+01 Au;
  137. y8y = +1.040615964042521*^+01 Au;
  138. z8z = -1.804702052122950*^-01 Au;
  139.  
  140. v8x = -2.089933372733080*^-03 Au/dy;
  141. v8y = +3.166549064213605*^-03 Au/dy;
  142. v8z = +3.884093561739733*^-05 Au/dy;
  143.  
  144. (* Neptun *)
  145.  
  146. m9 = +1.02413*^26 kg;
  147. q9 = +0 Amp sek;
  148.  
  149. x9x = +2.901867480863295*^+01 Au;
  150. y9y = -7.331260396521146*^+00 Au;
  151. z9z = -5.177914737734761*^-01 Au;
  152.  
  153. v9x = +7.476131405747911*^-04 Au/dy;
  154. v9y = +3.062101642790218*^-03 Au/dy;
  155. v9z = -8.000840096853115*^-05 Au/dy;
  156.  
  157. (* Differentialgleichung *)
  158.  
  159. nds=NDSolve[{
  160.  
  161. x1'[t] == vx1[t], y1'[t] == vy1[t], z1'[t] == vz1[t],
  162. x2'[t] == vx2[t], y2'[t] == vy2[t], z2'[t] == vz2[t],
  163. x3'[t] == vx3[t], y3'[t] == vy3[t], z3'[t] == vz3[t],
  164. x4'[t] == vx4[t], y4'[t] == vy4[t], z4'[t] == vz4[t],
  165. x5'[t] == vx5[t], y5'[t] == vy5[t], z5'[t] == vz5[t],
  166. x6'[t] == vx6[t], y6'[t] == vy6[t], z6'[t] == vz6[t],
  167. x7'[t] == vx7[t], y7'[t] == vy7[t], z7'[t] == vz7[t],
  168. x8'[t] == vx8[t], y8'[t] == vy8[t], z8'[t] == vz8[t],
  169. x9'[t] == vx9[t], y9'[t] == vy9[t], z9'[t] == vz9[t],
  170.  
  171. vx1'[t] ==
  172. (G m2 (x2[t]-x1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  173. (G m3 (x3[t]-x1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  174. (G m4 (x4[t]-x1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  175. (G m5 (x5[t]-x1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]+
  176. (G m6 (x6[t]-x1[t]))/Sqrt[((x6[t]-x1[t])^2+(y6[t]-y1[t])^2+(z6[t]-z1[t])^2)^3]+
  177. (G m7 (x7[t]-x1[t]))/Sqrt[((x7[t]-x1[t])^2+(y7[t]-y1[t])^2+(z7[t]-z1[t])^2)^3]+
  178. (G m8 (x8[t]-x1[t]))/Sqrt[((x8[t]-x1[t])^2+(y8[t]-y1[t])^2+(z8[t]-z1[t])^2)^3]+
  179. (G m9 (x9[t]-x1[t]))/Sqrt[((x9[t]-x1[t])^2+(y9[t]-y1[t])^2+(z9[t]-z1[t])^2)^3]+
  180. If[q1 == 0, 0,
  181. (-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]+
  182. (-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]+
  183. (-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]+
  184. (-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]+
  185. (-q1*q6/(4Pi ε0 )/m1 (x6[t]-x1[t]))/Sqrt[((x6[t]-x1[t])^2+(y6[t]-y1[t])^2+(z6[t]-z1[t])^2)^3]+
  186. (-q1*q7/(4Pi ε0 )/m1 (x7[t]-x1[t]))/Sqrt[((x7[t]-x1[t])^2+(y7[t]-y1[t])^2+(z7[t]-z1[t])^2)^3]+
  187. (-q1*q8/(4Pi ε0 )/m1 (x8[t]-x1[t]))/Sqrt[((x8[t]-x1[t])^2+(y8[t]-y1[t])^2+(z8[t]-z1[t])^2)^3]+
  188. (-q1*q9/(4Pi ε0 )/m1 (x9[t]-x1[t]))/Sqrt[((x9[t]-x1[t])^2+(y9[t]-y1[t])^2+(z9[t]-z1[t])^2)^3]]+
  189. Λ*c^2*x1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  190.  
  191. vy1'[t] ==
  192. (G m2 (y2[t]-y1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  193. (G m3 (y3[t]-y1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  194. (G m4 (y4[t]-y1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  195. (G m5 (y5[t]-y1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]+
  196. (G m6 (y6[t]-y1[t]))/Sqrt[((x6[t]-x1[t])^2+(y6[t]-y1[t])^2+(z6[t]-z1[t])^2)^3]+
  197. (G m7 (y7[t]-y1[t]))/Sqrt[((x7[t]-x1[t])^2+(y7[t]-y1[t])^2+(z7[t]-z1[t])^2)^3]+
  198. (G m8 (y8[t]-y1[t]))/Sqrt[((x8[t]-x1[t])^2+(y8[t]-y1[t])^2+(z8[t]-z1[t])^2)^3]+
  199. (G m9 (y9[t]-y1[t]))/Sqrt[((x9[t]-x1[t])^2+(y9[t]-y1[t])^2+(z9[t]-z1[t])^2)^3]+
  200. If[q1 == 0, 0,
  201. (-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]+
  202. (-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]+
  203. (-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]+
  204. (-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]+
  205. (-q1*q6/(4Pi ε0 )/m1 (y6[t]-y1[t]))/Sqrt[((x6[t]-x1[t])^2+(y6[t]-y1[t])^2+(z6[t]-z1[t])^2)^3]+
  206. (-q1*q7/(4Pi ε0 )/m1 (y7[t]-y1[t]))/Sqrt[((x7[t]-x1[t])^2+(y7[t]-y1[t])^2+(z7[t]-z1[t])^2)^3]+
  207. (-q1*q8/(4Pi ε0 )/m1 (y8[t]-y1[t]))/Sqrt[((x8[t]-x1[t])^2+(y8[t]-y1[t])^2+(z8[t]-z1[t])^2)^3]+
  208. (-q1*q9/(4Pi ε0 )/m1 (y9[t]-y1[t]))/Sqrt[((x9[t]-x1[t])^2+(y9[t]-y1[t])^2+(z9[t]-z1[t])^2)^3]]+
  209. Λ*c^2*y1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  210.  
  211. vz1'[t] ==
  212. (G m2 (z2[t]-z1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  213. (G m3 (z3[t]-z1[t]))/Sqrt[((x3[t]-x1[t])^2+(y3[t]-y1[t])^2+(z3[t]-z1[t])^2)^3]+
  214. (G m4 (z4[t]-z1[t]))/Sqrt[((x4[t]-x1[t])^2+(y4[t]-y1[t])^2+(z4[t]-z1[t])^2)^3]+
  215. (G m5 (z5[t]-z1[t]))/Sqrt[((x5[t]-x1[t])^2+(y5[t]-y1[t])^2+(z5[t]-z1[t])^2)^3]+
  216. (G m6 (z6[t]-z1[t]))/Sqrt[((x6[t]-x1[t])^2+(y6[t]-y1[t])^2+(z6[t]-z1[t])^2)^3]+
  217. (G m7 (z7[t]-z1[t]))/Sqrt[((x7[t]-x1[t])^2+(y7[t]-y1[t])^2+(z7[t]-z1[t])^2)^3]+
  218. (G m8 (z8[t]-z1[t]))/Sqrt[((x8[t]-x1[t])^2+(y8[t]-y1[t])^2+(z8[t]-z1[t])^2)^3]+
  219. (G m9 (z9[t]-z1[t]))/Sqrt[((x9[t]-x1[t])^2+(y9[t]-y1[t])^2+(z9[t]-z1[t])^2)^3]+
  220. If[q1 == 0, 0,
  221. (-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]+
  222. (-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]+
  223. (-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]+
  224. (-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]+
  225. (-q1*q6/(4Pi ε0 )/m1 (z6[t]-z1[t]))/Sqrt[((x6[t]-x1[t])^2+(y6[t]-y1[t])^2+(z6[t]-z1[t])^2)^3]+
  226. (-q1*q7/(4Pi ε0 )/m1 (z7[t]-z1[t]))/Sqrt[((x7[t]-x1[t])^2+(y7[t]-y1[t])^2+(z7[t]-z1[t])^2)^3]+
  227. (-q1*q8/(4Pi ε0 )/m1 (z8[t]-z1[t]))/Sqrt[((x8[t]-x1[t])^2+(y8[t]-y1[t])^2+(z8[t]-z1[t])^2)^3]+
  228. (-q1*q9/(4Pi ε0 )/m1 (z9[t]-z1[t]))/Sqrt[((x9[t]-x1[t])^2+(y9[t]-y1[t])^2+(z9[t]-z1[t])^2)^3]]+
  229. Λ*c^2*z1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
  230.  
  231. vx2'[t] ==
  232. (G m1 (x1[t]-x2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  233. (G m3 (x3[t]-x2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  234. (G m4 (x4[t]-x2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  235. (G m5 (x5[t]-x2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]+
  236. (G m6 (x6[t]-x2[t]))/Sqrt[((x6[t]-x2[t])^2+(y6[t]-y2[t])^2+(z6[t]-z2[t])^2)^3]+
  237. (G m7 (x7[t]-x2[t]))/Sqrt[((x7[t]-x2[t])^2+(y7[t]-y2[t])^2+(z7[t]-z2[t])^2)^3]+
  238. (G m8 (x8[t]-x2[t]))/Sqrt[((x8[t]-x2[t])^2+(y8[t]-y2[t])^2+(z8[t]-z2[t])^2)^3]+
  239. (G m9 (x9[t]-x2[t]))/Sqrt[((x9[t]-x2[t])^2+(y9[t]-y2[t])^2+(z9[t]-z2[t])^2)^3]+
  240. If[q2 == 0, 0,
  241. (-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]+
  242. (-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]+
  243. (-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]+
  244. (-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]+
  245. (-q2*q6/(4Pi ε0 )/m2 (x6[t]-x2[t]))/Sqrt[((x6[t]-x2[t])^2+(y6[t]-y2[t])^2+(z6[t]-z2[t])^2)^3]+
  246. (-q2*q7/(4Pi ε0 )/m2 (x7[t]-x2[t]))/Sqrt[((x7[t]-x2[t])^2+(y7[t]-y2[t])^2+(z7[t]-z2[t])^2)^3]+
  247. (-q2*q8/(4Pi ε0 )/m2 (x8[t]-x2[t]))/Sqrt[((x8[t]-x2[t])^2+(y8[t]-y2[t])^2+(z8[t]-z2[t])^2)^3]+
  248. (-q2*q9/(4Pi ε0 )/m2 (x9[t]-x2[t]))/Sqrt[((x9[t]-x2[t])^2+(y9[t]-y2[t])^2+(z9[t]-z2[t])^2)^3]]+
  249. Λ*c^2*x2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  250.  
  251. vy2'[t] ==
  252. (G m1 (y1[t]-y2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
  253. (G m3 (y3[t]-y2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  254. (G m4 (y4[t]-y2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  255. (G m5 (y5[t]-y2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]+
  256. (G m6 (y6[t]-y2[t]))/Sqrt[((x6[t]-x2[t])^2+(y6[t]-y2[t])^2+(z6[t]-z2[t])^2)^3]+
  257. (G m7 (y7[t]-y2[t]))/Sqrt[((x7[t]-x2[t])^2+(y7[t]-y2[t])^2+(z7[t]-z2[t])^2)^3]+
  258. (G m8 (y8[t]-y2[t]))/Sqrt[((x8[t]-x2[t])^2+(y8[t]-y2[t])^2+(z8[t]-z2[t])^2)^3]+
  259. (G m9 (y9[t]-y2[t]))/Sqrt[((x9[t]-x2[t])^2+(y9[t]-y2[t])^2+(z9[t]-z2[t])^2)^3]+
  260. If[q2 == 0, 0,
  261. (-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]+
  262. (-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]+
  263. (-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]+
  264. (-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]+
  265. (-q2*q6/(4Pi ε0 )/m2 (y6[t]-y2[t]))/Sqrt[((x6[t]-x2[t])^2+(y6[t]-y2[t])^2+(z6[t]-z2[t])^2)^3]+
  266. (-q2*q7/(4Pi ε0 )/m2 (y7[t]-y2[t]))/Sqrt[((x7[t]-x2[t])^2+(y7[t]-y2[t])^2+(z7[t]-z2[t])^2)^3]+
  267. (-q2*q8/(4Pi ε0 )/m2 (y8[t]-y2[t]))/Sqrt[((x8[t]-x2[t])^2+(y8[t]-y2[t])^2+(z8[t]-z2[t])^2)^3]+
  268. (-q2*q9/(4Pi ε0 )/m2 (y9[t]-y2[t]))/Sqrt[((x9[t]-x2[t])^2+(y9[t]-y2[t])^2+(z9[t]-z2[t])^2)^3]]+
  269. Λ*c^2*y2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  270.  
  271. vz2'[t] ==
  272. (G m1 (z1[t]-z2[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
  273. (G m3 (z3[t]-z2[t]))/Sqrt[((x3[t]-x2[t])^2+(y3[t]-y2[t])^2+(z3[t]-z2[t])^2)^3]+
  274. (G m4 (z4[t]-z2[t]))/Sqrt[((x4[t]-x2[t])^2+(y4[t]-y2[t])^2+(z4[t]-z2[t])^2)^3]+
  275. (G m5 (z5[t]-z2[t]))/Sqrt[((x5[t]-x2[t])^2+(y5[t]-y2[t])^2+(z5[t]-z2[t])^2)^3]+
  276. (G m6 (z6[t]-z2[t]))/Sqrt[((x6[t]-x2[t])^2+(y6[t]-y2[t])^2+(z6[t]-z2[t])^2)^3]+
  277. (G m7 (z7[t]-z2[t]))/Sqrt[((x7[t]-x2[t])^2+(y7[t]-y2[t])^2+(z7[t]-z2[t])^2)^3]+
  278. (G m8 (z8[t]-z2[t]))/Sqrt[((x8[t]-x2[t])^2+(y8[t]-y2[t])^2+(z8[t]-z2[t])^2)^3]+
  279. (G m9 (z9[t]-z2[t]))/Sqrt[((x9[t]-x2[t])^2+(y9[t]-y2[t])^2+(z9[t]-z2[t])^2)^3]+
  280. If[q2 == 0, 0,
  281. (-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]+
  282. (-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]+
  283. (-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]+
  284. (-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]+
  285. (-q2*q6/(4Pi ε0 )/m2 (z6[t]-z2[t]))/Sqrt[((x6[t]-x2[t])^2+(y6[t]-y2[t])^2+(z6[t]-z2[t])^2)^3]+
  286. (-q2*q7/(4Pi ε0 )/m2 (z7[t]-z2[t]))/Sqrt[((x7[t]-x2[t])^2+(y7[t]-y2[t])^2+(z7[t]-z2[t])^2)^3]+
  287. (-q2*q8/(4Pi ε0 )/m2 (z8[t]-z2[t]))/Sqrt[((x8[t]-x2[t])^2+(y8[t]-y2[t])^2+(z8[t]-z2[t])^2)^3]+
  288. (-q2*q9/(4Pi ε0 )/m2 (z9[t]-z2[t]))/Sqrt[((x9[t]-x2[t])^2+(y9[t]-y2[t])^2+(z9[t]-z2[t])^2)^3]]+
  289. Λ*c^2*z2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
  290.  
  291. vx3'[t] ==
  292. (G m1 (x1[t]-x3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  293. (G m2 (x2[t]-x3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  294. (G m4 (x4[t]-x3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  295. (G m5 (x5[t]-x3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]+
  296. (G m6 (x6[t]-x3[t]))/Sqrt[((x6[t]-x3[t])^2+(y6[t]-y3[t])^2+(z6[t]-z3[t])^2)^3]+
  297. (G m7 (x7[t]-x3[t]))/Sqrt[((x7[t]-x3[t])^2+(y7[t]-y3[t])^2+(z7[t]-z3[t])^2)^3]+
  298. (G m8 (x8[t]-x3[t]))/Sqrt[((x8[t]-x3[t])^2+(y8[t]-y3[t])^2+(z8[t]-z3[t])^2)^3]+
  299. (G m9 (x9[t]-x3[t]))/Sqrt[((x9[t]-x3[t])^2+(y9[t]-y3[t])^2+(z9[t]-z3[t])^2)^3]+
  300. If[q3 == 0, 0,
  301. (-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]+
  302. (-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]+
  303. (-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]+
  304. (-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]+
  305. (-q3*q6/(4Pi ε0 )/m3 (x6[t]-x3[t]))/Sqrt[((x6[t]-x3[t])^2+(y6[t]-y3[t])^2+(z6[t]-z3[t])^2)^3]+
  306. (-q3*q7/(4Pi ε0 )/m3 (x7[t]-x3[t]))/Sqrt[((x7[t]-x3[t])^2+(y7[t]-y3[t])^2+(z7[t]-z3[t])^2)^3]+
  307. (-q3*q8/(4Pi ε0 )/m3 (x8[t]-x3[t]))/Sqrt[((x8[t]-x3[t])^2+(y8[t]-y3[t])^2+(z8[t]-z3[t])^2)^3]+
  308. (-q3*q9/(4Pi ε0 )/m3 (x9[t]-x3[t]))/Sqrt[((x9[t]-x3[t])^2+(y9[t]-y3[t])^2+(z9[t]-z3[t])^2)^3]]+
  309. Λ*c^2*x3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  310.  
  311. vy3'[t] ==
  312. (G m1 (y1[t]-y3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  313. (G m2 (y2[t]-y3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  314. (G m4 (y4[t]-y3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  315. (G m5 (y5[t]-y3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]+
  316. (G m6 (y6[t]-y3[t]))/Sqrt[((x6[t]-x3[t])^2+(y6[t]-y3[t])^2+(z6[t]-z3[t])^2)^3]+
  317. (G m7 (y7[t]-y3[t]))/Sqrt[((x7[t]-x3[t])^2+(y7[t]-y3[t])^2+(z7[t]-z3[t])^2)^3]+
  318. (G m8 (y8[t]-y3[t]))/Sqrt[((x8[t]-x3[t])^2+(y8[t]-y3[t])^2+(z8[t]-z3[t])^2)^3]+
  319. (G m9 (y9[t]-y3[t]))/Sqrt[((x9[t]-x3[t])^2+(y9[t]-y3[t])^2+(z9[t]-z3[t])^2)^3]+
  320. If[q3 == 0, 0,
  321. (-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]+
  322. (-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]+
  323. (-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]+
  324. (-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]+
  325. (-q3*q6/(4Pi ε0 )/m3 (y6[t]-y3[t]))/Sqrt[((x6[t]-x3[t])^2+(y6[t]-y3[t])^2+(z6[t]-z3[t])^2)^3]+
  326. (-q3*q7/(4Pi ε0 )/m3 (y7[t]-y3[t]))/Sqrt[((x7[t]-x3[t])^2+(y7[t]-y3[t])^2+(z7[t]-z3[t])^2)^3]+
  327. (-q3*q8/(4Pi ε0 )/m3 (y8[t]-y3[t]))/Sqrt[((x8[t]-x3[t])^2+(y8[t]-y3[t])^2+(z8[t]-z3[t])^2)^3]+
  328. (-q3*q9/(4Pi ε0 )/m3 (y9[t]-y3[t]))/Sqrt[((x9[t]-x3[t])^2+(y9[t]-y3[t])^2+(z9[t]-z3[t])^2)^3]]+
  329. Λ*c^2*y3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  330.  
  331. vz3'[t] ==
  332. (G m1 (z1[t]-z3[t]))/Sqrt[((x1[t]-x3[t])^2+(y1[t]-y3[t])^2+(z1[t]-z3[t])^2)^3]+
  333. (G m2 (z2[t]-z3[t]))/Sqrt[((x2[t]-x3[t])^2+(y2[t]-y3[t])^2+(z2[t]-z3[t])^2)^3]+
  334. (G m4 (z4[t]-z3[t]))/Sqrt[((x4[t]-x3[t])^2+(y4[t]-y3[t])^2+(z4[t]-z3[t])^2)^3]+
  335. (G m5 (z5[t]-z3[t]))/Sqrt[((x5[t]-x3[t])^2+(y5[t]-y3[t])^2+(z5[t]-z3[t])^2)^3]+
  336. (G m6 (z6[t]-z3[t]))/Sqrt[((x6[t]-x3[t])^2+(y6[t]-y3[t])^2+(z6[t]-z3[t])^2)^3]+
  337. (G m7 (z7[t]-z3[t]))/Sqrt[((x7[t]-x3[t])^2+(y7[t]-y3[t])^2+(z7[t]-z3[t])^2)^3]+
  338. (G m8 (z8[t]-z3[t]))/Sqrt[((x8[t]-x3[t])^2+(y8[t]-y3[t])^2+(z8[t]-z3[t])^2)^3]+
  339. (G m9 (z9[t]-z3[t]))/Sqrt[((x9[t]-x3[t])^2+(y9[t]-y3[t])^2+(z9[t]-z3[t])^2)^3]+
  340. If[q3 == 0, 0,
  341. (-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]+
  342. (-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]+
  343. (-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]+
  344. (-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]+
  345. (-q3*q6/(4Pi ε0 )/m3 (z6[t]-z3[t]))/Sqrt[((x6[t]-x3[t])^2+(y6[t]-y3[t])^2+(z6[t]-z3[t])^2)^3]+
  346. (-q3*q7/(4Pi ε0 )/m3 (z7[t]-z3[t]))/Sqrt[((x7[t]-x3[t])^2+(y7[t]-y3[t])^2+(z7[t]-z3[t])^2)^3]+
  347. (-q3*q8/(4Pi ε0 )/m3 (z8[t]-z3[t]))/Sqrt[((x8[t]-x3[t])^2+(y8[t]-y3[t])^2+(z8[t]-z3[t])^2)^3]+
  348. (-q3*q9/(4Pi ε0 )/m3 (z9[t]-z3[t]))/Sqrt[((x9[t]-x3[t])^2+(y9[t]-y3[t])^2+(z9[t]-z3[t])^2)^3]]+
  349. Λ*c^2*z3[t]^2/Sqrt[x3[t]^2+y3[t]^2+z3[t]^2],
  350.  
  351. vx4'[t] ==
  352. (G m1 (x1[t]-x4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  353. (G m2 (x2[t]-x4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  354. (G m3 (x3[t]-x4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  355. (G m5 (x5[t]-x4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]+
  356. (G m6 (x6[t]-x4[t]))/Sqrt[((x6[t]-x4[t])^2+(y6[t]-y4[t])^2+(z6[t]-z4[t])^2)^3]+
  357. (G m7 (x7[t]-x4[t]))/Sqrt[((x7[t]-x4[t])^2+(y7[t]-y4[t])^2+(z7[t]-z4[t])^2)^3]+
  358. (G m8 (x8[t]-x4[t]))/Sqrt[((x8[t]-x4[t])^2+(y8[t]-y4[t])^2+(z8[t]-z4[t])^2)^3]+
  359. (G m9 (x9[t]-x4[t]))/Sqrt[((x9[t]-x4[t])^2+(y9[t]-y4[t])^2+(z9[t]-z4[t])^2)^3]+
  360. If[q4 == 0, 0,
  361. (-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]+
  362. (-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]+
  363. (-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]+
  364. (-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]+
  365. (-q4*q6/(4Pi ε0 )/m4 (x6[t]-x4[t]))/Sqrt[((x6[t]-x4[t])^2+(y6[t]-y4[t])^2+(z6[t]-z4[t])^2)^3]+
  366. (-q4*q7/(4Pi ε0 )/m4 (x7[t]-x4[t]))/Sqrt[((x7[t]-x4[t])^2+(y7[t]-y4[t])^2+(z7[t]-z4[t])^2)^3]+
  367. (-q4*q8/(4Pi ε0 )/m4 (x8[t]-x4[t]))/Sqrt[((x8[t]-x4[t])^2+(y8[t]-y4[t])^2+(z8[t]-z4[t])^2)^3]+
  368. (-q4*q9/(4Pi ε0 )/m4 (x9[t]-x4[t]))/Sqrt[((x9[t]-x4[t])^2+(y9[t]-y4[t])^2+(z9[t]-z4[t])^2)^3]]+
  369. Λ*c^2*x4[t]^2/Sqrt[x4[t]^2+y4[t]^2+z4[t]^2],
  370.  
  371. vy4'[t] ==
  372. (G m1 (y1[t]-y4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  373. (G m2 (y2[t]-y4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  374. (G m3 (y3[t]-y4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  375. (G m5 (y5[t]-y4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]+
  376. (G m6 (y6[t]-y4[t]))/Sqrt[((x6[t]-x4[t])^2+(y6[t]-y4[t])^2+(z6[t]-z4[t])^2)^3]+
  377. (G m7 (y7[t]-y4[t]))/Sqrt[((x7[t]-x4[t])^2+(y7[t]-y4[t])^2+(z7[t]-z4[t])^2)^3]+
  378. (G m8 (y8[t]-y4[t]))/Sqrt[((x8[t]-x4[t])^2+(y8[t]-y4[t])^2+(z8[t]-z4[t])^2)^3]+
  379. (G m9 (y9[t]-y4[t]))/Sqrt[((x9[t]-x4[t])^2+(y9[t]-y4[t])^2+(z9[t]-z4[t])^2)^3]+
  380. If[q4 == 0, 0,
  381. (-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]+
  382. (-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]+
  383. (-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]+
  384. (-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]+
  385. (-q4*q6/(4Pi ε0 )/m4 (y6[t]-y4[t]))/Sqrt[((x6[t]-x4[t])^2+(y6[t]-y4[t])^2+(z6[t]-z4[t])^2)^3]+
  386. (-q4*q7/(4Pi ε0 )/m4 (y7[t]-y4[t]))/Sqrt[((x7[t]-x4[t])^2+(y7[t]-y4[t])^2+(z7[t]-z4[t])^2)^3]+
  387. (-q4*q8/(4Pi ε0 )/m4 (y8[t]-y4[t]))/Sqrt[((x8[t]-x4[t])^2+(y8[t]-y4[t])^2+(z8[t]-z4[t])^2)^3]+
  388. (-q4*q9/(4Pi ε0 )/m4 (y9[t]-y4[t]))/Sqrt[((x9[t]-x4[t])^2+(y9[t]-y4[t])^2+(z9[t]-z4[t])^2)^3]]+
  389. Λ*c^2*y4[t]^2/Sqrt[x4[t]^2+y4[t]^2+z4[t]^2],
  390.  
  391. vz4'[t] ==
  392. (G m1 (z1[t]-z4[t]))/Sqrt[((x1[t]-x4[t])^2+(y1[t]-y4[t])^2+(z1[t]-z4[t])^2)^3]+
  393. (G m2 (z2[t]-z4[t]))/Sqrt[((x2[t]-x4[t])^2+(y2[t]-y4[t])^2+(z2[t]-z4[t])^2)^3]+
  394. (G m3 (z3[t]-z4[t]))/Sqrt[((x3[t]-x4[t])^2+(y3[t]-y4[t])^2+(z3[t]-z4[t])^2)^3]+
  395. (G m5 (z5[t]-z4[t]))/Sqrt[((x5[t]-x4[t])^2+(y5[t]-y4[t])^2+(z5[t]-z4[t])^2)^3]+
  396. (G m6 (z6[t]-z4[t]))/Sqrt[((x6[t]-x4[t])^2+(y6[t]-y4[t])^2+(z6[t]-z4[t])^2)^3]+
  397. (G m7 (z7[t]-z4[t]))/Sqrt[((x7[t]-x4[t])^2+(y7[t]-y4[t])^2+(z7[t]-z4[t])^2)^3]+
  398. (G m8 (z8[t]-z4[t]))/Sqrt[((x8[t]-x4[t])^2+(y8[t]-y4[t])^2+(z8[t]-z4[t])^2)^3]+
  399. (G m9 (z9[t]-z4[t]))/Sqrt[((x9[t]-x4[t])^2+(y9[t]-y4[t])^2+(z9[t]-z4[t])^2)^3]+
  400. If[q4 == 0, 0,
  401. (-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]+
  402. (-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]+
  403. (-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]+
  404. (-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]+
  405. (-q4*q6/(4Pi ε0 )/m4 (z6[t]-z4[t]))/Sqrt[((x6[t]-x4[t])^2+(y6[t]-y4[t])^2+(z6[t]-z4[t])^2)^3]+
  406. (-q4*q7/(4Pi ε0 )/m4 (z7[t]-z4[t]))/Sqrt[((x7[t]-x4[t])^2+(y7[t]-y4[t])^2+(z7[t]-z4[t])^2)^3]+
  407. (-q4*q8/(4Pi ε0 )/m4 (z8[t]-z4[t]))/Sqrt[((x8[t]-x4[t])^2+(y8[t]-y4[t])^2+(z8[t]-z4[t])^2)^3]+
  408. (-q4*q9/(4Pi ε0 )/m4 (z9[t]-z4[t]))/Sqrt[((x9[t]-x4[t])^2+(y9[t]-y4[t])^2+(z9[t]-z4[t])^2)^3]]+
  409. Λ*c^2*z4[t]^2/Sqrt[x4[t]^2+y4[t]^2+z4[t]^2],
  410.  
  411. vx5'[t] ==
  412. (G m1 (x1[t]-x5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  413. (G m2 (x2[t]-x5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  414. (G m3 (x3[t]-x5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  415. (G m4 (x4[t]-x5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]+
  416. (G m6 (x6[t]-x5[t]))/Sqrt[((x6[t]-x5[t])^2+(y6[t]-y5[t])^2+(z6[t]-z5[t])^2)^3]+
  417. (G m7 (x7[t]-x5[t]))/Sqrt[((x7[t]-x5[t])^2+(y7[t]-y5[t])^2+(z7[t]-z5[t])^2)^3]+
  418. (G m8 (x8[t]-x5[t]))/Sqrt[((x8[t]-x5[t])^2+(y8[t]-y5[t])^2+(z8[t]-z5[t])^2)^3]+
  419. (G m9 (x9[t]-x5[t]))/Sqrt[((x9[t]-x5[t])^2+(y9[t]-y5[t])^2+(z9[t]-z5[t])^2)^3]+
  420. If[q5 == 0, 0,
  421. (-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]+
  422. (-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]+
  423. (-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]+
  424. (-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]+
  425. (-q5*q6/(4Pi ε0 )/m5 (x6[t]-x5[t]))/Sqrt[((x6[t]-x5[t])^2+(y6[t]-y5[t])^2+(z6[t]-z5[t])^2)^3]+
  426. (-q5*q7/(4Pi ε0 )/m5 (x7[t]-x5[t]))/Sqrt[((x7[t]-x5[t])^2+(y7[t]-y5[t])^2+(z7[t]-z5[t])^2)^3]+
  427. (-q5*q8/(4Pi ε0 )/m5 (x8[t]-x5[t]))/Sqrt[((x8[t]-x5[t])^2+(y8[t]-y5[t])^2+(z8[t]-z5[t])^2)^3]+
  428. (-q5*q9/(4Pi ε0 )/m5 (x9[t]-x5[t]))/Sqrt[((x9[t]-x5[t])^2+(y9[t]-y5[t])^2+(z9[t]-z5[t])^2)^3]]+
  429. Λ*c^2*x5[t]^2/Sqrt[x5[t]^2+y5[t]^2+z5[t]^2],
  430.  
  431. vy5'[t] ==
  432. (G m1 (y1[t]-y5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  433. (G m2 (y2[t]-y5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  434. (G m3 (y3[t]-y5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  435. (G m4 (y4[t]-y5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]+
  436. (G m6 (y6[t]-y5[t]))/Sqrt[((x6[t]-x5[t])^2+(y6[t]-y5[t])^2+(z6[t]-z5[t])^2)^3]+
  437. (G m7 (y7[t]-y5[t]))/Sqrt[((x7[t]-x5[t])^2+(y7[t]-y5[t])^2+(z7[t]-z5[t])^2)^3]+
  438. (G m8 (y8[t]-y5[t]))/Sqrt[((x8[t]-x5[t])^2+(y8[t]-y5[t])^2+(z8[t]-z5[t])^2)^3]+
  439. (G m9 (y9[t]-y5[t]))/Sqrt[((x9[t]-x5[t])^2+(y9[t]-y5[t])^2+(z9[t]-z5[t])^2)^3]+
  440. If[q5 == 0, 0,
  441. (-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]+
  442. (-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]+
  443. (-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]+
  444. (-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]+
  445. (-q5*q6/(4Pi ε0 )/m5 (y6[t]-y5[t]))/Sqrt[((x6[t]-x5[t])^2+(y6[t]-y5[t])^2+(z6[t]-z5[t])^2)^3]+
  446. (-q5*q7/(4Pi ε0 )/m5 (y7[t]-y5[t]))/Sqrt[((x7[t]-x5[t])^2+(y7[t]-y5[t])^2+(z7[t]-z5[t])^2)^3]+
  447. (-q5*q8/(4Pi ε0 )/m5 (y8[t]-y5[t]))/Sqrt[((x8[t]-x5[t])^2+(y8[t]-y5[t])^2+(z8[t]-z5[t])^2)^3]+
  448. (-q5*q9/(4Pi ε0 )/m5 (y9[t]-y5[t]))/Sqrt[((x9[t]-x5[t])^2+(y9[t]-y5[t])^2+(z9[t]-z5[t])^2)^3]]+
  449. Λ*c^2*y5[t]^2/Sqrt[x5[t]^2+y5[t]^2+z5[t]^2],
  450.  
  451. vz5'[t] ==
  452. (G m1 (z1[t]-z5[t]))/Sqrt[((x1[t]-x5[t])^2+(y1[t]-y5[t])^2+(z1[t]-z5[t])^2)^3]+
  453. (G m2 (z2[t]-z5[t]))/Sqrt[((x2[t]-x5[t])^2+(y2[t]-y5[t])^2+(z2[t]-z5[t])^2)^3]+
  454. (G m3 (z3[t]-z5[t]))/Sqrt[((x3[t]-x5[t])^2+(y3[t]-y5[t])^2+(z3[t]-z5[t])^2)^3]+
  455. (G m4 (z4[t]-z5[t]))/Sqrt[((x4[t]-x5[t])^2+(y4[t]-y5[t])^2+(z4[t]-z5[t])^2)^3]+
  456. (G m6 (z6[t]-z5[t]))/Sqrt[((x6[t]-x5[t])^2+(y6[t]-y5[t])^2+(z6[t]-z5[t])^2)^3]+
  457. (G m7 (z7[t]-z5[t]))/Sqrt[((x7[t]-x5[t])^2+(y7[t]-y5[t])^2+(z7[t]-z5[t])^2)^3]+
  458. (G m8 (z8[t]-z5[t]))/Sqrt[((x8[t]-x5[t])^2+(y8[t]-y5[t])^2+(z8[t]-z5[t])^2)^3]+
  459. (G m9 (z9[t]-z5[t]))/Sqrt[((x9[t]-x5[t])^2+(y9[t]-y5[t])^2+(z9[t]-z5[t])^2)^3]+
  460. If[q5 == 0, 0,
  461. (-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]+
  462. (-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]+
  463. (-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]+
  464. (-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]+
  465. (-q5*q6/(4Pi ε0 )/m5 (z6[t]-z5[t]))/Sqrt[((x6[t]-x5[t])^2+(y6[t]-y5[t])^2+(z6[t]-z5[t])^2)^3]+
  466. (-q5*q7/(4Pi ε0 )/m5 (z7[t]-z5[t]))/Sqrt[((x7[t]-x5[t])^2+(y7[t]-y5[t])^2+(z7[t]-z5[t])^2)^3]+
  467. (-q5*q8/(4Pi ε0 )/m5 (z8[t]-z5[t]))/Sqrt[((x8[t]-x5[t])^2+(y8[t]-y5[t])^2+(z8[t]-z5[t])^2)^3]+
  468. (-q5*q9/(4Pi ε0 )/m5 (z9[t]-z5[t]))/Sqrt[((x9[t]-x5[t])^2+(y9[t]-y5[t])^2+(z9[t]-z5[t])^2)^3]]+
  469. Λ*c^2*z5[t]^2/Sqrt[x5[t]^2+y5[t]^2+z5[t]^2],
  470.  
  471. vx6'[t] ==
  472. (G m1 (x1[t]-x6[t]))/Sqrt[((x1[t]-x6[t])^2+(y1[t]-y6[t])^2+(z1[t]-z6[t])^2)^3]+
  473. (G m2 (x2[t]-x6[t]))/Sqrt[((x2[t]-x6[t])^2+(y2[t]-y6[t])^2+(z2[t]-z6[t])^2)^3]+
  474. (G m3 (x3[t]-x6[t]))/Sqrt[((x3[t]-x6[t])^2+(y3[t]-y6[t])^2+(z3[t]-z6[t])^2)^3]+
  475. (G m4 (x4[t]-x6[t]))/Sqrt[((x4[t]-x6[t])^2+(y4[t]-y6[t])^2+(z4[t]-z6[t])^2)^3]+
  476. (G m5 (x5[t]-x6[t]))/Sqrt[((x5[t]-x6[t])^2+(y5[t]-y6[t])^2+(z5[t]-z6[t])^2)^3]+
  477. (G m7 (x7[t]-x6[t]))/Sqrt[((x7[t]-x6[t])^2+(y7[t]-y6[t])^2+(z7[t]-z6[t])^2)^3]+
  478. (G m8 (x8[t]-x6[t]))/Sqrt[((x8[t]-x6[t])^2+(y8[t]-y6[t])^2+(z8[t]-z6[t])^2)^3]+
  479. (G m9 (x9[t]-x6[t]))/Sqrt[((x9[t]-x6[t])^2+(y9[t]-y6[t])^2+(z9[t]-z6[t])^2)^3]+
  480. If[q6 == 0, 0,
  481. (-q6*q1/(4Pi ε0 )/m6 (x1[t]-x6[t]))/Sqrt[((x1[t]-x6[t])^2+(y1[t]-y6[t])^2+(z1[t]-z6[t])^2)^3]+
  482. (-q6*q2/(4Pi ε0 )/m6 (x2[t]-x6[t]))/Sqrt[((x2[t]-x6[t])^2+(y2[t]-y6[t])^2+(z2[t]-z6[t])^2)^3]+
  483. (-q6*q3/(4Pi ε0 )/m6 (x3[t]-x6[t]))/Sqrt[((x3[t]-x6[t])^2+(y3[t]-y6[t])^2+(z3[t]-z6[t])^2)^3]+
  484. (-q6*q4/(4Pi ε0 )/m6 (x4[t]-x6[t]))/Sqrt[((x4[t]-x6[t])^2+(y4[t]-y6[t])^2+(z4[t]-z6[t])^2)^3]+
  485. (-q6*q5/(4Pi ε0 )/m6 (x5[t]-x6[t]))/Sqrt[((x5[t]-x6[t])^2+(y5[t]-y6[t])^2+(z5[t]-z6[t])^2)^3]+
  486. (-q6*q7/(4Pi ε0 )/m6 (x7[t]-x6[t]))/Sqrt[((x7[t]-x6[t])^2+(y7[t]-y6[t])^2+(z7[t]-z6[t])^2)^3]+
  487. (-q6*q8/(4Pi ε0 )/m6 (x8[t]-x6[t]))/Sqrt[((x8[t]-x6[t])^2+(y8[t]-y6[t])^2+(z8[t]-z6[t])^2)^3]+
  488. (-q6*q9/(4Pi ε0 )/m6 (x9[t]-x6[t]))/Sqrt[((x9[t]-x6[t])^2+(y9[t]-y6[t])^2+(z9[t]-z6[t])^2)^3]]+
  489. Λ*c^2*x6[t]^2/Sqrt[x6[t]^2+y6[t]^2+z6[t]^2],
  490.  
  491. vy6'[t] ==
  492. (G m1 (y1[t]-y6[t]))/Sqrt[((x1[t]-x6[t])^2+(y1[t]-y6[t])^2+(z1[t]-z6[t])^2)^3]+
  493. (G m2 (y2[t]-y6[t]))/Sqrt[((x2[t]-x6[t])^2+(y2[t]-y6[t])^2+(z2[t]-z6[t])^2)^3]+
  494. (G m3 (y3[t]-y6[t]))/Sqrt[((x3[t]-x6[t])^2+(y3[t]-y6[t])^2+(z3[t]-z6[t])^2)^3]+
  495. (G m4 (y4[t]-y6[t]))/Sqrt[((x4[t]-x6[t])^2+(y4[t]-y6[t])^2+(z4[t]-z6[t])^2)^3]+
  496. (G m5 (y5[t]-y6[t]))/Sqrt[((x5[t]-x6[t])^2+(y5[t]-y6[t])^2+(z5[t]-z6[t])^2)^3]+
  497. (G m7 (y7[t]-y6[t]))/Sqrt[((x7[t]-x6[t])^2+(y7[t]-y6[t])^2+(z7[t]-z6[t])^2)^3]+
  498. (G m8 (y8[t]-y6[t]))/Sqrt[((x8[t]-x6[t])^2+(y8[t]-y6[t])^2+(z8[t]-z6[t])^2)^3]+
  499. (G m9 (y9[t]-y6[t]))/Sqrt[((x9[t]-x6[t])^2+(y9[t]-y6[t])^2+(z9[t]-z6[t])^2)^3]+
  500. If[q6 == 0, 0,
  501. (-q6*q1/(4Pi ε0 )/m6 (y1[t]-y6[t]))/Sqrt[((x1[t]-x6[t])^2+(y1[t]-y6[t])^2+(z1[t]-z6[t])^2)^3]+
  502. (-q6*q2/(4Pi ε0 )/m6 (y2[t]-y6[t]))/Sqrt[((x2[t]-x6[t])^2+(y2[t]-y6[t])^2+(z2[t]-z6[t])^2)^3]+
  503. (-q6*q3/(4Pi ε0 )/m6 (y3[t]-y6[t]))/Sqrt[((x3[t]-x6[t])^2+(y3[t]-y6[t])^2+(z3[t]-z6[t])^2)^3]+
  504. (-q6*q4/(4Pi ε0 )/m6 (y4[t]-y6[t]))/Sqrt[((x4[t]-x6[t])^2+(y4[t]-y6[t])^2+(z4[t]-z6[t])^2)^3]+
  505. (-q6*q5/(4Pi ε0 )/m6 (y5[t]-y6[t]))/Sqrt[((x5[t]-x6[t])^2+(y5[t]-y6[t])^2+(z5[t]-z6[t])^2)^3]+
  506. (-q6*q7/(4Pi ε0 )/m6 (y7[t]-y6[t]))/Sqrt[((x7[t]-x6[t])^2+(y7[t]-y6[t])^2+(z7[t]-z6[t])^2)^3]+
  507. (-q6*q8/(4Pi ε0 )/m6 (y8[t]-y6[t]))/Sqrt[((x8[t]-x6[t])^2+(y8[t]-y6[t])^2+(z8[t]-z6[t])^2)^3]+
  508. (-q6*q9/(4Pi ε0 )/m6 (y9[t]-y6[t]))/Sqrt[((x9[t]-x6[t])^2+(y9[t]-y6[t])^2+(z9[t]-z6[t])^2)^3]]+
  509. Λ*c^2*y6[t]^2/Sqrt[x6[t]^2+y6[t]^2+z6[t]^2],
  510.  
  511. vz6'[t] ==
  512. (G m1 (z1[t]-z6[t]))/Sqrt[((x1[t]-x6[t])^2+(y1[t]-y6[t])^2+(z1[t]-z6[t])^2)^3]+
  513. (G m2 (z2[t]-z6[t]))/Sqrt[((x2[t]-x6[t])^2+(y2[t]-y6[t])^2+(z2[t]-z6[t])^2)^3]+
  514. (G m3 (z3[t]-z6[t]))/Sqrt[((x3[t]-x6[t])^2+(y3[t]-y6[t])^2+(z3[t]-z6[t])^2)^3]+
  515. (G m4 (z4[t]-z6[t]))/Sqrt[((x4[t]-x6[t])^2+(y4[t]-y6[t])^2+(z4[t]-z6[t])^2)^3]+
  516. (G m5 (z5[t]-z6[t]))/Sqrt[((x5[t]-x6[t])^2+(y5[t]-y6[t])^2+(z5[t]-z6[t])^2)^3]+
  517. (G m7 (z7[t]-z6[t]))/Sqrt[((x7[t]-x6[t])^2+(y7[t]-y6[t])^2+(z7[t]-z6[t])^2)^3]+
  518. (G m8 (z8[t]-z6[t]))/Sqrt[((x8[t]-x6[t])^2+(y8[t]-y6[t])^2+(z8[t]-z6[t])^2)^3]+
  519. (G m9 (z9[t]-z6[t]))/Sqrt[((x9[t]-x6[t])^2+(y9[t]-y6[t])^2+(z9[t]-z6[t])^2)^3]+
  520. If[q6 == 0, 0,
  521. (-q6*q1/(4Pi ε0 )/m6 (z1[t]-z6[t]))/Sqrt[((x1[t]-x6[t])^2+(y1[t]-y6[t])^2+(z1[t]-z6[t])^2)^3]+
  522. (-q6*q2/(4Pi ε0 )/m6 (z2[t]-z6[t]))/Sqrt[((x2[t]-x6[t])^2+(y2[t]-y6[t])^2+(z2[t]-z6[t])^2)^3]+
  523. (-q6*q3/(4Pi ε0 )/m6 (z3[t]-z6[t]))/Sqrt[((x3[t]-x6[t])^2+(y3[t]-y6[t])^2+(z3[t]-z6[t])^2)^3]+
  524. (-q6*q4/(4Pi ε0 )/m6 (z4[t]-z6[t]))/Sqrt[((x4[t]-x6[t])^2+(y4[t]-y6[t])^2+(z4[t]-z6[t])^2)^3]+
  525. (-q6*q5/(4Pi ε0 )/m6 (z5[t]-z6[t]))/Sqrt[((x5[t]-x6[t])^2+(y5[t]-y6[t])^2+(z5[t]-z6[t])^2)^3]+
  526. (-q6*q7/(4Pi ε0 )/m6 (z7[t]-z6[t]))/Sqrt[((x7[t]-x6[t])^2+(y7[t]-y6[t])^2+(z7[t]-z6[t])^2)^3]+
  527. (-q6*q8/(4Pi ε0 )/m6 (z8[t]-z6[t]))/Sqrt[((x8[t]-x6[t])^2+(y8[t]-y6[t])^2+(z8[t]-z6[t])^2)^3]+
  528. (-q6*q9/(4Pi ε0 )/m6 (z9[t]-z6[t]))/Sqrt[((x9[t]-x6[t])^2+(y9[t]-y6[t])^2+(z9[t]-z6[t])^2)^3]]+
  529. Λ*c^2*z6[t]^2/Sqrt[x6[t]^2+y6[t]^2+z6[t]^2],
  530.  
  531. vx7'[t] ==
  532. (G m1 (x1[t]-x7[t]))/Sqrt[((x1[t]-x7[t])^2+(y1[t]-y7[t])^2+(z1[t]-z7[t])^2)^3]+
  533. (G m2 (x2[t]-x7[t]))/Sqrt[((x2[t]-x7[t])^2+(y2[t]-y7[t])^2+(z2[t]-z7[t])^2)^3]+
  534. (G m3 (x3[t]-x7[t]))/Sqrt[((x3[t]-x7[t])^2+(y3[t]-y7[t])^2+(z3[t]-z7[t])^2)^3]+
  535. (G m4 (x4[t]-x7[t]))/Sqrt[((x4[t]-x7[t])^2+(y4[t]-y7[t])^2+(z4[t]-z7[t])^2)^3]+
  536. (G m5 (x5[t]-x7[t]))/Sqrt[((x5[t]-x7[t])^2+(y5[t]-y7[t])^2+(z5[t]-z7[t])^2)^3]+
  537. (G m6 (x6[t]-x7[t]))/Sqrt[((x6[t]-x7[t])^2+(y6[t]-y7[t])^2+(z6[t]-z7[t])^2)^3]+
  538. (G m8 (x8[t]-x7[t]))/Sqrt[((x8[t]-x7[t])^2+(y8[t]-y7[t])^2+(z8[t]-z7[t])^2)^3]+
  539. (G m9 (x9[t]-x7[t]))/Sqrt[((x9[t]-x7[t])^2+(y9[t]-y7[t])^2+(z9[t]-z7[t])^2)^3]+
  540. If[q7 == 0, 0,
  541. (-q7*q1/(4Pi ε0 )/m7 (x1[t]-x7[t]))/Sqrt[((x1[t]-x7[t])^2+(y1[t]-y7[t])^2+(z1[t]-z7[t])^2)^3]+
  542. (-q7*q2/(4Pi ε0 )/m7 (x2[t]-x7[t]))/Sqrt[((x2[t]-x7[t])^2+(y2[t]-y7[t])^2+(z2[t]-z7[t])^2)^3]+
  543. (-q7*q3/(4Pi ε0 )/m7 (x3[t]-x7[t]))/Sqrt[((x3[t]-x7[t])^2+(y3[t]-y7[t])^2+(z3[t]-z7[t])^2)^3]+
  544. (-q7*q4/(4Pi ε0 )/m7 (x4[t]-x7[t]))/Sqrt[((x4[t]-x7[t])^2+(y4[t]-y7[t])^2+(z4[t]-z7[t])^2)^3]+
  545. (-q7*q5/(4Pi ε0 )/m7 (x5[t]-x7[t]))/Sqrt[((x5[t]-x7[t])^2+(y5[t]-y7[t])^2+(z5[t]-z7[t])^2)^3]+
  546. (-q7*q6/(4Pi ε0 )/m7 (x6[t]-x7[t]))/Sqrt[((x6[t]-x7[t])^2+(y6[t]-y7[t])^2+(z6[t]-z7[t])^2)^3]+
  547. (-q7*q8/(4Pi ε0 )/m7 (x8[t]-x7[t]))/Sqrt[((x8[t]-x7[t])^2+(y8[t]-y7[t])^2+(z8[t]-z7[t])^2)^3]+
  548. (-q7*q9/(4Pi ε0 )/m7 (x9[t]-x7[t]))/Sqrt[((x9[t]-x7[t])^2+(y9[t]-y7[t])^2+(z9[t]-z7[t])^2)^3]]+
  549. Λ*c^2*x7[t]^2/Sqrt[x7[t]^2+y7[t]^2+z7[t]^2],
  550.  
  551. vy7'[t] ==
  552. (G m1 (y1[t]-y7[t]))/Sqrt[((x1[t]-x7[t])^2+(y1[t]-y7[t])^2+(z1[t]-z7[t])^2)^3]+
  553. (G m2 (y2[t]-y7[t]))/Sqrt[((x2[t]-x7[t])^2+(y2[t]-y7[t])^2+(z2[t]-z7[t])^2)^3]+
  554. (G m3 (y3[t]-y7[t]))/Sqrt[((x3[t]-x7[t])^2+(y3[t]-y7[t])^2+(z3[t]-z7[t])^2)^3]+
  555. (G m4 (y4[t]-y7[t]))/Sqrt[((x4[t]-x7[t])^2+(y4[t]-y7[t])^2+(z4[t]-z7[t])^2)^3]+
  556. (G m5 (y5[t]-y7[t]))/Sqrt[((x5[t]-x7[t])^2+(y5[t]-y7[t])^2+(z5[t]-z7[t])^2)^3]+
  557. (G m6 (y6[t]-y7[t]))/Sqrt[((x6[t]-x7[t])^2+(y6[t]-y7[t])^2+(z6[t]-z7[t])^2)^3]+
  558. (G m8 (y8[t]-y7[t]))/Sqrt[((x8[t]-x7[t])^2+(y8[t]-y7[t])^2+(z8[t]-z7[t])^2)^3]+
  559. (G m9 (y9[t]-y7[t]))/Sqrt[((x9[t]-x7[t])^2+(y9[t]-y7[t])^2+(z9[t]-z7[t])^2)^3]+
  560. If[q7 == 0, 0,
  561. (-q7*q1/(4Pi ε0 )/m7 (y1[t]-y7[t]))/Sqrt[((x1[t]-x7[t])^2+(y1[t]-y7[t])^2+(z1[t]-z7[t])^2)^3]+
  562. (-q7*q2/(4Pi ε0 )/m7 (y2[t]-y7[t]))/Sqrt[((x2[t]-x7[t])^2+(y2[t]-y7[t])^2+(z2[t]-z7[t])^2)^3]+
  563. (-q7*q3/(4Pi ε0 )/m7 (y3[t]-y7[t]))/Sqrt[((x3[t]-x7[t])^2+(y3[t]-y7[t])^2+(z3[t]-z7[t])^2)^3]+
  564. (-q7*q4/(4Pi ε0 )/m7 (y4[t]-y7[t]))/Sqrt[((x4[t]-x7[t])^2+(y4[t]-y7[t])^2+(z4[t]-z7[t])^2)^3]+
  565. (-q7*q5/(4Pi ε0 )/m7 (y5[t]-y7[t]))/Sqrt[((x5[t]-x7[t])^2+(y5[t]-y7[t])^2+(z5[t]-z7[t])^2)^3]+
  566. (-q7*q6/(4Pi ε0 )/m7 (y6[t]-y7[t]))/Sqrt[((x6[t]-x7[t])^2+(y6[t]-y7[t])^2+(z6[t]-z7[t])^2)^3]+
  567. (-q7*q8/(4Pi ε0 )/m7 (y8[t]-y7[t]))/Sqrt[((x8[t]-x7[t])^2+(y8[t]-y7[t])^2+(z8[t]-z7[t])^2)^3]+
  568. (-q7*q9/(4Pi ε0 )/m7 (y9[t]-y7[t]))/Sqrt[((x9[t]-x7[t])^2+(y9[t]-y7[t])^2+(z9[t]-z7[t])^2)^3]]+
  569. Λ*c^2*y7[t]^2/Sqrt[x7[t]^2+y7[t]^2+z7[t]^2],
  570.  
  571. vz7'[t] ==
  572. (G m1 (z1[t]-z7[t]))/Sqrt[((x1[t]-x7[t])^2+(y1[t]-y7[t])^2+(z1[t]-z7[t])^2)^3]+
  573. (G m2 (z2[t]-z7[t]))/Sqrt[((x2[t]-x7[t])^2+(y2[t]-y7[t])^2+(z2[t]-z7[t])^2)^3]+
  574. (G m3 (z3[t]-z7[t]))/Sqrt[((x3[t]-x7[t])^2+(y3[t]-y7[t])^2+(z3[t]-z7[t])^2)^3]+
  575. (G m4 (z4[t]-z7[t]))/Sqrt[((x4[t]-x7[t])^2+(y4[t]-y7[t])^2+(z4[t]-z7[t])^2)^3]+
  576. (G m5 (z5[t]-z7[t]))/Sqrt[((x5[t]-x7[t])^2+(y5[t]-y7[t])^2+(z5[t]-z7[t])^2)^3]+
  577. (G m6 (z6[t]-z7[t]))/Sqrt[((x6[t]-x7[t])^2+(y6[t]-y7[t])^2+(z6[t]-z7[t])^2)^3]+
  578. (G m8 (z8[t]-z7[t]))/Sqrt[((x8[t]-x7[t])^2+(y8[t]-y7[t])^2+(z8[t]-z7[t])^2)^3]+
  579. (G m9 (z9[t]-z7[t]))/Sqrt[((x9[t]-x7[t])^2+(y9[t]-y7[t])^2+(z9[t]-z7[t])^2)^3]+
  580. If[q7 == 0, 0,
  581. (-q7*q1/(4Pi ε0 )/m7 (z1[t]-z7[t]))/Sqrt[((x1[t]-x7[t])^2+(y1[t]-y7[t])^2+(z1[t]-z7[t])^2)^3]+
  582. (-q7*q2/(4Pi ε0 )/m7 (z2[t]-z7[t]))/Sqrt[((x2[t]-x7[t])^2+(y2[t]-y7[t])^2+(z2[t]-z7[t])^2)^3]+
  583. (-q7*q3/(4Pi ε0 )/m7 (z3[t]-z7[t]))/Sqrt[((x3[t]-x7[t])^2+(y3[t]-y7[t])^2+(z3[t]-z7[t])^2)^3]+
  584. (-q7*q4/(4Pi ε0 )/m7 (z4[t]-z7[t]))/Sqrt[((x4[t]-x7[t])^2+(y4[t]-y7[t])^2+(z4[t]-z7[t])^2)^3]+
  585. (-q7*q5/(4Pi ε0 )/m7 (z5[t]-z7[t]))/Sqrt[((x5[t]-x7[t])^2+(y5[t]-y7[t])^2+(z5[t]-z7[t])^2)^3]+
  586. (-q7*q6/(4Pi ε0 )/m7 (z6[t]-z7[t]))/Sqrt[((x6[t]-x7[t])^2+(y6[t]-y7[t])^2+(z6[t]-z7[t])^2)^3]+
  587. (-q7*q8/(4Pi ε0 )/m7 (z8[t]-z7[t]))/Sqrt[((x8[t]-x7[t])^2+(y8[t]-y7[t])^2+(z8[t]-z7[t])^2)^3]+
  588. (-q7*q9/(4Pi ε0 )/m7 (z9[t]-z7[t]))/Sqrt[((x9[t]-x7[t])^2+(y9[t]-y7[t])^2+(z9[t]-z7[t])^2)^3]]+
  589. Λ*c^2*z7[t]^2/Sqrt[x7[t]^2+y7[t]^2+z7[t]^2],
  590.  
  591. vx8'[t] ==
  592. (G m1 (x1[t]-x8[t]))/Sqrt[((x1[t]-x8[t])^2+(y1[t]-y8[t])^2+(z1[t]-z8[t])^2)^3]+
  593. (G m2 (x2[t]-x8[t]))/Sqrt[((x2[t]-x8[t])^2+(y2[t]-y8[t])^2+(z2[t]-z8[t])^2)^3]+
  594. (G m3 (x3[t]-x8[t]))/Sqrt[((x3[t]-x8[t])^2+(y3[t]-y8[t])^2+(z3[t]-z8[t])^2)^3]+
  595. (G m4 (x4[t]-x8[t]))/Sqrt[((x4[t]-x8[t])^2+(y4[t]-y8[t])^2+(z4[t]-z8[t])^2)^3]+
  596. (G m5 (x5[t]-x8[t]))/Sqrt[((x5[t]-x8[t])^2+(y5[t]-y8[t])^2+(z5[t]-z8[t])^2)^3]+
  597. (G m6 (x6[t]-x8[t]))/Sqrt[((x6[t]-x8[t])^2+(y6[t]-y8[t])^2+(z6[t]-z8[t])^2)^3]+
  598. (G m7 (x7[t]-x8[t]))/Sqrt[((x7[t]-x8[t])^2+(y7[t]-y8[t])^2+(z7[t]-z8[t])^2)^3]+
  599. (G m9 (x9[t]-x8[t]))/Sqrt[((x9[t]-x8[t])^2+(y9[t]-y8[t])^2+(z9[t]-z8[t])^2)^3]+
  600. If[q8 == 0, 0,
  601. (-q8*q1/(4Pi ε0 )/m8 (x1[t]-x8[t]))/Sqrt[((x1[t]-x8[t])^2+(y1[t]-y8[t])^2+(z1[t]-z8[t])^2)^3]+
  602. (-q8*q2/(4Pi ε0 )/m8 (x2[t]-x8[t]))/Sqrt[((x2[t]-x8[t])^2+(y2[t]-y8[t])^2+(z2[t]-z8[t])^2)^3]+
  603. (-q8*q3/(4Pi ε0 )/m8 (x3[t]-x8[t]))/Sqrt[((x3[t]-x8[t])^2+(y3[t]-y8[t])^2+(z3[t]-z8[t])^2)^3]+
  604. (-q8*q4/(4Pi ε0 )/m8 (x4[t]-x8[t]))/Sqrt[((x4[t]-x8[t])^2+(y4[t]-y8[t])^2+(z4[t]-z8[t])^2)^3]+
  605. (-q8*q5/(4Pi ε0 )/m8 (x5[t]-x8[t]))/Sqrt[((x5[t]-x8[t])^2+(y5[t]-y8[t])^2+(z5[t]-z8[t])^2)^3]+
  606. (-q8*q6/(4Pi ε0 )/m8 (x6[t]-x8[t]))/Sqrt[((x6[t]-x8[t])^2+(y6[t]-y8[t])^2+(z6[t]-z8[t])^2)^3]+
  607. (-q8*q7/(4Pi ε0 )/m8 (x7[t]-x8[t]))/Sqrt[((x7[t]-x8[t])^2+(y7[t]-y8[t])^2+(z7[t]-z8[t])^2)^3]+
  608. (-q8*q9/(4Pi ε0 )/m8 (x9[t]-x8[t]))/Sqrt[((x9[t]-x8[t])^2+(y9[t]-y8[t])^2+(z9[t]-z8[t])^2)^3]]+
  609. Λ*c^2*x8[t]^2/Sqrt[x8[t]^2+y8[t]^2+z8[t]^2],
  610.  
  611. vy8'[t] ==
  612. (G m1 (y1[t]-y8[t]))/Sqrt[((x1[t]-x8[t])^2+(y1[t]-y8[t])^2+(z1[t]-z8[t])^2)^3]+
  613. (G m2 (y2[t]-y8[t]))/Sqrt[((x2[t]-x8[t])^2+(y2[t]-y8[t])^2+(z2[t]-z8[t])^2)^3]+
  614. (G m3 (y3[t]-y8[t]))/Sqrt[((x3[t]-x8[t])^2+(y3[t]-y8[t])^2+(z3[t]-z8[t])^2)^3]+
  615. (G m4 (y4[t]-y8[t]))/Sqrt[((x4[t]-x8[t])^2+(y4[t]-y8[t])^2+(z4[t]-z8[t])^2)^3]+
  616. (G m5 (y5[t]-y8[t]))/Sqrt[((x5[t]-x8[t])^2+(y5[t]-y8[t])^2+(z5[t]-z8[t])^2)^3]+
  617. (G m6 (y6[t]-y8[t]))/Sqrt[((x6[t]-x8[t])^2+(y6[t]-y8[t])^2+(z6[t]-z8[t])^2)^3]+
  618. (G m7 (y7[t]-y8[t]))/Sqrt[((x7[t]-x8[t])^2+(y7[t]-y8[t])^2+(z7[t]-z8[t])^2)^3]+
  619. (G m9 (y9[t]-y8[t]))/Sqrt[((x9[t]-x8[t])^2+(y9[t]-y8[t])^2+(z9[t]-z8[t])^2)^3]+
  620. If[q8 == 0, 0,
  621. (-q8*q1/(4Pi ε0 )/m8 (y1[t]-y8[t]))/Sqrt[((x1[t]-x8[t])^2+(y1[t]-y8[t])^2+(z1[t]-z8[t])^2)^3]+
  622. (-q8*q2/(4Pi ε0 )/m8 (y2[t]-y8[t]))/Sqrt[((x2[t]-x8[t])^2+(y2[t]-y8[t])^2+(z2[t]-z8[t])^2)^3]+
  623. (-q8*q3/(4Pi ε0 )/m8 (y3[t]-y8[t]))/Sqrt[((x3[t]-x8[t])^2+(y3[t]-y8[t])^2+(z3[t]-z8[t])^2)^3]+
  624. (-q8*q4/(4Pi ε0 )/m8 (y4[t]-y8[t]))/Sqrt[((x4[t]-x8[t])^2+(y4[t]-y8[t])^2+(z4[t]-z8[t])^2)^3]+
  625. (-q8*q5/(4Pi ε0 )/m8 (y5[t]-y8[t]))/Sqrt[((x5[t]-x8[t])^2+(y5[t]-y8[t])^2+(z5[t]-z8[t])^2)^3]+
  626. (-q8*q6/(4Pi ε0 )/m8 (y6[t]-y8[t]))/Sqrt[((x6[t]-x8[t])^2+(y6[t]-y8[t])^2+(z6[t]-z8[t])^2)^3]+
  627. (-q8*q7/(4Pi ε0 )/m8 (y7[t]-y8[t]))/Sqrt[((x7[t]-x8[t])^2+(y7[t]-y8[t])^2+(z7[t]-z8[t])^2)^3]+
  628. (-q8*q9/(4Pi ε0 )/m8 (y9[t]-y8[t]))/Sqrt[((x9[t]-x8[t])^2+(y9[t]-y8[t])^2+(z9[t]-z8[t])^2)^3]]+
  629. Λ*c^2*y8[t]^2/Sqrt[x8[t]^2+y8[t]^2+z8[t]^2],
  630.  
  631. vz8'[t] ==
  632. (G m1 (z1[t]-z8[t]))/Sqrt[((x1[t]-x8[t])^2+(y1[t]-y8[t])^2+(z1[t]-z8[t])^2)^3]+
  633. (G m2 (z2[t]-z8[t]))/Sqrt[((x2[t]-x8[t])^2+(y2[t]-y8[t])^2+(z2[t]-z8[t])^2)^3]+
  634. (G m3 (z3[t]-z8[t]))/Sqrt[((x3[t]-x8[t])^2+(y3[t]-y8[t])^2+(z3[t]-z8[t])^2)^3]+
  635. (G m4 (z4[t]-z8[t]))/Sqrt[((x4[t]-x8[t])^2+(y4[t]-y8[t])^2+(z4[t]-z8[t])^2)^3]+
  636. (G m5 (z5[t]-z8[t]))/Sqrt[((x5[t]-x8[t])^2+(y5[t]-y8[t])^2+(z5[t]-z8[t])^2)^3]+
  637. (G m6 (z6[t]-z8[t]))/Sqrt[((x6[t]-x8[t])^2+(y6[t]-y8[t])^2+(z6[t]-z8[t])^2)^3]+
  638. (G m7 (z7[t]-z8[t]))/Sqrt[((x7[t]-x8[t])^2+(y7[t]-y8[t])^2+(z7[t]-z8[t])^2)^3]+
  639. (G m9 (z9[t]-z8[t]))/Sqrt[((x9[t]-x8[t])^2+(y9[t]-y8[t])^2+(z9[t]-z8[t])^2)^3]+
  640. If[q8 == 0, 0,
  641. (-q8*q1/(4Pi ε0 )/m8 (z1[t]-z8[t]))/Sqrt[((x1[t]-x8[t])^2+(y1[t]-y8[t])^2+(z1[t]-z8[t])^2)^3]+
  642. (-q8*q2/(4Pi ε0 )/m8 (z2[t]-z8[t]))/Sqrt[((x2[t]-x8[t])^2+(y2[t]-y8[t])^2+(z2[t]-z8[t])^2)^3]+
  643. (-q8*q3/(4Pi ε0 )/m8 (z3[t]-z8[t]))/Sqrt[((x3[t]-x8[t])^2+(y3[t]-y8[t])^2+(z3[t]-z8[t])^2)^3]+
  644. (-q8*q4/(4Pi ε0 )/m8 (z4[t]-z8[t]))/Sqrt[((x4[t]-x8[t])^2+(y4[t]-y8[t])^2+(z4[t]-z8[t])^2)^3]+
  645. (-q8*q5/(4Pi ε0 )/m8 (z5[t]-z8[t]))/Sqrt[((x5[t]-x8[t])^2+(y5[t]-y8[t])^2+(z5[t]-z8[t])^2)^3]+
  646. (-q8*q6/(4Pi ε0 )/m8 (z6[t]-z8[t]))/Sqrt[((x6[t]-x8[t])^2+(y6[t]-y8[t])^2+(z6[t]-z8[t])^2)^3]+
  647. (-q8*q7/(4Pi ε0 )/m8 (z7[t]-z8[t]))/Sqrt[((x7[t]-x8[t])^2+(y7[t]-y8[t])^2+(z7[t]-z8[t])^2)^3]+
  648. (-q8*q9/(4Pi ε0 )/m8 (z9[t]-z8[t]))/Sqrt[((x9[t]-x8[t])^2+(y9[t]-y8[t])^2+(z9[t]-z8[t])^2)^3]]+
  649. Λ*c^2*z8[t]^2/Sqrt[x8[t]^2+y8[t]^2+z8[t]^2],
  650.  
  651. vx9'[t] ==
  652. (G m1 (x1[t]-x9[t]))/Sqrt[((x1[t]-x9[t])^2+(y1[t]-y9[t])^2+(z1[t]-z9[t])^2)^3]+
  653. (G m2 (x2[t]-x9[t]))/Sqrt[((x2[t]-x9[t])^2+(y2[t]-y9[t])^2+(z2[t]-z9[t])^2)^3]+
  654. (G m3 (x3[t]-x9[t]))/Sqrt[((x3[t]-x9[t])^2+(y3[t]-y9[t])^2+(z3[t]-z9[t])^2)^3]+
  655. (G m4 (x4[t]-x9[t]))/Sqrt[((x4[t]-x9[t])^2+(y4[t]-y9[t])^2+(z4[t]-z9[t])^2)^3]+
  656. (G m5 (x5[t]-x9[t]))/Sqrt[((x5[t]-x9[t])^2+(y5[t]-y9[t])^2+(z5[t]-z9[t])^2)^3]+
  657. (G m6 (x6[t]-x9[t]))/Sqrt[((x6[t]-x9[t])^2+(y6[t]-y9[t])^2+(z6[t]-z9[t])^2)^3]+
  658. (G m7 (x7[t]-x9[t]))/Sqrt[((x7[t]-x9[t])^2+(y7[t]-y9[t])^2+(z7[t]-z9[t])^2)^3]+
  659. (G m8 (x8[t]-x9[t]))/Sqrt[((x8[t]-x9[t])^2+(y8[t]-y9[t])^2+(z8[t]-z9[t])^2)^3]+
  660. If[q9 == 0, 0,
  661. (-q9*q1/(4Pi ε0 )/m9 (x1[t]-x9[t]))/Sqrt[((x1[t]-x9[t])^2+(y1[t]-y9[t])^2+(z1[t]-z9[t])^2)^3]+
  662. (-q9*q2/(4Pi ε0 )/m9 (x2[t]-x9[t]))/Sqrt[((x2[t]-x9[t])^2+(y2[t]-y9[t])^2+(z2[t]-z9[t])^2)^3]+
  663. (-q9*q3/(4Pi ε0 )/m9 (x3[t]-x9[t]))/Sqrt[((x3[t]-x9[t])^2+(y3[t]-y9[t])^2+(z3[t]-z9[t])^2)^3]+
  664. (-q9*q4/(4Pi ε0 )/m9 (x4[t]-x9[t]))/Sqrt[((x4[t]-x9[t])^2+(y4[t]-y9[t])^2+(z4[t]-z9[t])^2)^3]+
  665. (-q9*q5/(4Pi ε0 )/m9 (x5[t]-x9[t]))/Sqrt[((x5[t]-x9[t])^2+(y5[t]-y9[t])^2+(z5[t]-z9[t])^2)^3]+
  666. (-q9*q6/(4Pi ε0 )/m9 (x6[t]-x9[t]))/Sqrt[((x6[t]-x9[t])^2+(y6[t]-y9[t])^2+(z6[t]-z9[t])^2)^3]+
  667. (-q9*q7/(4Pi ε0 )/m9 (x7[t]-x9[t]))/Sqrt[((x7[t]-x9[t])^2+(y7[t]-y9[t])^2+(z7[t]-z9[t])^2)^3]+
  668. (-q9*q8/(4Pi ε0 )/m9 (x8[t]-x9[t]))/Sqrt[((x8[t]-x9[t])^2+(y8[t]-y9[t])^2+(z8[t]-z9[t])^2)^3]]+
  669. Λ*c^2*x9[t]^2/Sqrt[x9[t]^2+y9[t]^2+z9[t]^2],
  670.  
  671. vy9'[t] ==
  672. (G m1 (y1[t]-y9[t]))/Sqrt[((x1[t]-x9[t])^2+(y1[t]-y9[t])^2+(z1[t]-z9[t])^2)^3]+
  673. (G m2 (y2[t]-y9[t]))/Sqrt[((x2[t]-x9[t])^2+(y2[t]-y9[t])^2+(z2[t]-z9[t])^2)^3]+
  674. (G m3 (y3[t]-y9[t]))/Sqrt[((x3[t]-x9[t])^2+(y3[t]-y9[t])^2+(z3[t]-z9[t])^2)^3]+
  675. (G m4 (y4[t]-y9[t]))/Sqrt[((x4[t]-x9[t])^2+(y4[t]-y9[t])^2+(z4[t]-z9[t])^2)^3]+
  676. (G m5 (y5[t]-y9[t]))/Sqrt[((x5[t]-x9[t])^2+(y5[t]-y9[t])^2+(z5[t]-z9[t])^2)^3]+
  677. (G m6 (y6[t]-y9[t]))/Sqrt[((x6[t]-x9[t])^2+(y6[t]-y9[t])^2+(z6[t]-z9[t])^2)^3]+
  678. (G m7 (y7[t]-y9[t]))/Sqrt[((x7[t]-x9[t])^2+(y7[t]-y9[t])^2+(z7[t]-z9[t])^2)^3]+
  679. (G m8 (y8[t]-y9[t]))/Sqrt[((x8[t]-x9[t])^2+(y8[t]-y9[t])^2+(z8[t]-z9[t])^2)^3]+
  680. If[q9 == 0, 0,
  681. (-q9*q1/(4Pi ε0 )/m9 (y1[t]-y9[t]))/Sqrt[((x1[t]-x9[t])^2+(y1[t]-y9[t])^2+(z1[t]-z9[t])^2)^3]+
  682. (-q9*q2/(4Pi ε0 )/m9 (y2[t]-y9[t]))/Sqrt[((x2[t]-x9[t])^2+(y2[t]-y9[t])^2+(z2[t]-z9[t])^2)^3]+
  683. (-q9*q3/(4Pi ε0 )/m9 (y3[t]-y9[t]))/Sqrt[((x3[t]-x9[t])^2+(y3[t]-y9[t])^2+(z3[t]-z9[t])^2)^3]+
  684. (-q9*q4/(4Pi ε0 )/m9 (y4[t]-y9[t]))/Sqrt[((x4[t]-x9[t])^2+(y4[t]-y9[t])^2+(z4[t]-z9[t])^2)^3]+
  685. (-q9*q5/(4Pi ε0 )/m9 (y5[t]-y9[t]))/Sqrt[((x5[t]-x9[t])^2+(y5[t]-y9[t])^2+(z5[t]-z9[t])^2)^3]+
  686. (-q9*q6/(4Pi ε0 )/m9 (y6[t]-y9[t]))/Sqrt[((x6[t]-x9[t])^2+(y6[t]-y9[t])^2+(z6[t]-z9[t])^2)^3]+
  687. (-q9*q7/(4Pi ε0 )/m9 (y7[t]-y9[t]))/Sqrt[((x7[t]-x9[t])^2+(y7[t]-y9[t])^2+(z7[t]-z9[t])^2)^3]+
  688. (-q9*q8/(4Pi ε0 )/m9 (y8[t]-y9[t]))/Sqrt[((x8[t]-x9[t])^2+(y8[t]-y9[t])^2+(z8[t]-z9[t])^2)^3]]+
  689. Λ*c^2*y9[t]^2/Sqrt[x9[t]^2+y9[t]^2+z9[t]^2],
  690.  
  691. vz9'[t] ==
  692. (G m1 (z1[t]-z9[t]))/Sqrt[((x1[t]-x9[t])^2+(y1[t]-y9[t])^2+(z1[t]-z9[t])^2)^3]+
  693. (G m2 (z2[t]-z9[t]))/Sqrt[((x2[t]-x9[t])^2+(y2[t]-y9[t])^2+(z2[t]-z9[t])^2)^3]+
  694. (G m3 (z3[t]-z9[t]))/Sqrt[((x3[t]-x9[t])^2+(y3[t]-y9[t])^2+(z3[t]-z9[t])^2)^3]+
  695. (G m4 (z4[t]-z9[t]))/Sqrt[((x4[t]-x9[t])^2+(y4[t]-y9[t])^2+(z4[t]-z9[t])^2)^3]+
  696. (G m5 (z5[t]-z9[t]))/Sqrt[((x5[t]-x9[t])^2+(y5[t]-y9[t])^2+(z5[t]-z9[t])^2)^3]+
  697. (G m6 (z6[t]-z9[t]))/Sqrt[((x6[t]-x9[t])^2+(y6[t]-y9[t])^2+(z6[t]-z9[t])^2)^3]+
  698. (G m7 (z7[t]-z9[t]))/Sqrt[((x7[t]-x9[t])^2+(y7[t]-y9[t])^2+(z7[t]-z9[t])^2)^3]+
  699. (G m8 (z8[t]-z9[t]))/Sqrt[((x8[t]-x9[t])^2+(y8[t]-y9[t])^2+(z8[t]-z9[t])^2)^3]+
  700. If[q9 == 0, 0,
  701. (-q9*q1/(4Pi ε0 )/m9 (z1[t]-z9[t]))/Sqrt[((x1[t]-x9[t])^2+(y1[t]-y9[t])^2+(z1[t]-z9[t])^2)^3]+
  702. (-q9*q2/(4Pi ε0 )/m9 (z2[t]-z9[t]))/Sqrt[((x2[t]-x9[t])^2+(y2[t]-y9[t])^2+(z2[t]-z9[t])^2)^3]+
  703. (-q9*q3/(4Pi ε0 )/m9 (z3[t]-z9[t]))/Sqrt[((x3[t]-x9[t])^2+(y3[t]-y9[t])^2+(z3[t]-z9[t])^2)^3]+
  704. (-q9*q4/(4Pi ε0 )/m9 (z4[t]-z9[t]))/Sqrt[((x4[t]-x9[t])^2+(y4[t]-y9[t])^2+(z4[t]-z9[t])^2)^3]+
  705. (-q9*q5/(4Pi ε0 )/m9 (z5[t]-z9[t]))/Sqrt[((x5[t]-x9[t])^2+(y5[t]-y9[t])^2+(z5[t]-z9[t])^2)^3]+
  706. (-q9*q6/(4Pi ε0 )/m9 (z6[t]-z9[t]))/Sqrt[((x6[t]-x9[t])^2+(y6[t]-y9[t])^2+(z6[t]-z9[t])^2)^3]+
  707. (-q9*q7/(4Pi ε0 )/m9 (z7[t]-z9[t]))/Sqrt[((x7[t]-x9[t])^2+(y7[t]-y9[t])^2+(z7[t]-z9[t])^2)^3]+
  708. (-q9*q8/(4Pi ε0 )/m9 (z8[t]-z9[t]))/Sqrt[((x8[t]-x9[t])^2+(y8[t]-y9[t])^2+(z8[t]-z9[t])^2)^3]]+
  709. Λ*c^2*z9[t]^2/Sqrt[x9[t]^2+y9[t]^2+z9[t]^2],
  710.  
  711. x1[0] == x1x, y1[0] == y1y, z1[0] == z1z,
  712. x2[0] == x2x, y2[0] == y2y, z2[0] == z2z,
  713. x3[0] == x3x, y3[0] == y3y, z3[0] == z3z,
  714. x4[0] == x4x, y4[0] == y4y, z4[0] == z4z,
  715. x5[0] == x5x, y5[0] == y5y, z5[0] == z5z,
  716. x6[0] == x6x, y6[0] == y6y, z6[0] == z6z,
  717. x7[0] == x7x, y7[0] == y7y, z7[0] == z7z,
  718. x8[0] == x8x, y8[0] == y8y, z8[0] == z8z,
  719. x9[0] == x9x, y9[0] == y9y, z9[0] == z9z,
  720.  
  721. vx1[0] == v1x, vy1[0] == v1y, vz1[0] == v1z,
  722. vx2[0] == v2x, vy2[0] == v2y, vz2[0] == v2z,
  723. vx3[0] == v3x, vy3[0] == v3y, vz3[0] == v3z,
  724. vx4[0] == v4x, vy4[0] == v4y, vz4[0] == v4z,
  725. vx5[0] == v5x, vy5[0] == v5y, vz5[0] == v5z,
  726. vx6[0] == v6x, vy6[0] == v6y, vz6[0] == v6z,
  727. vx7[0] == v7x, vy7[0] == v7y, vz7[0] == v7z,
  728. vx8[0] == v8x, vy8[0] == v8y, vz8[0] == v8z,
  729. vx9[0] == v9x, vy9[0] == v9y, vz9[0] == v9z},
  730.  
  731. {x1, x2, x3, x4, x5, x6, x7, x8, x9, y1, y2, y3, y4, y5, y6, y7, y8, y9, z1, z2, z3, z4, z5, z6, z7, z8, z9,
  732. vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8, vx9, vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8, vy9, vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8, vz9},
  733.  
  734. {t, 0, Tmax},
  735.  
  736. WorkingPrecision-> wp,
  737. MaxSteps-> Infinity,
  738. Method-> mta,
  739. InterpolationOrder-> All,
  740. StepMonitor :> (laststep=plunge; plunge=t;
  741. stepsize=plunge-laststep;), Method->{"EventLocator",
  742. "Event" :> (If[stepsize<1*^-4, 0, 1])}];
  743.  
  744. (* Position, Geschwindigkeit *)
  745.  
  746. 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]}, {x6[t], y6[t], z6[t]}, {x7[t], y7[t], z7[t]}, {x8[t], y8[t], z8[t]}, {x9[t], y9[t], z9[t]}}/.nds[[1]];
  747. 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]}, {vx6[t], vy6[t], vz6[t]}, {vx7[t], vy7[t], vz7[t]}, {vx8[t], vy8[t], vz8[t]}, {vx9[t], vy9[t], vz9[t]}}/.nds[[1]];
  748. 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]]]+m6 Evaluate[f2p[t][[6]]]+m7 Evaluate[f2p[t][[7]]]+m8 Evaluate[f2p[t][[8]]]+m9 Evaluate[f2p[t][[9]]])/(m1+m2+m3+m4+m5+m6+m7+m8+m9);
  749.  
  750. (* Formatierung *)
  751.  
  752. s[text_]=Style[text, FontSize->11];
  753. sw[text_]=Style[text, White, FontSize->11];
  754. colorfunc[n_]=Function[{x, y, z, t},
  755. Hue[0, n, 0.5,
  756. If[Tmax<0, Max[Min[(+T+(-t+trail))/trail, 1], 0],
  757. Max[Min[(-T+(t+trail))/trail, 1], 0]]]];
  758.  
  759. (* Animation *)
  760.  
  761. Do[Print[Rasterize[
  762. Grid[{{
  763. Show[
  764.  
  765. If[T == 0, {},
  766.  
  767. ParametricPlot3D[Evaluate[f2p[t]],
  768. {t, Max[0, T-trail], T},
  769.  
  770. PlotStyle->{
  771. {Thickness[thk], Red},
  772. {Thickness[thk], Blue},
  773. {Thickness[thk], Green},
  774. {Thickness[thk], Magenta},
  775. {Thickness[thk], Cyan},
  776. {Thickness[thk], Orange},
  777. {Thickness[thk], Purple},
  778. {Thickness[thk], Pink},
  779. {Thickness[thk], Brown}},
  780.  
  781. PlotRange->plotrange, AspectRatio->1, MaxRecursion->15, Axes->True, ImageSize->imagesize]],
  782.  
  783. Graphics3D[
  784. If[startpos==1, {
  785. {PointSize[2point/3], Lighter[Red], Point[{x1x, y1y, z1z}]},
  786. {PointSize[2point/3], Lighter[Blue], Point[{x2x, y2y, z2z}]},
  787. {PointSize[2point/3], Lighter[Green], Point[{x3x, y3y, z3z}]},
  788. {PointSize[2point/3], Lighter[Magenta], Point[{x4x, y4y, z4z}]},
  789. {PointSize[2point/3], Lighter[Cyan], Point[{x5x, y5y, z5z}]},
  790. {PointSize[2point/3], Lighter[Orange], Point[{x6x, y6y, z6z}]},
  791. {PointSize[2point/3], Lighter[Purple], Point[{x7x, y7y, z7z}]},
  792. {PointSize[2point/3], Lighter[Pink], Point[{x8x, y8y, z8z}]},
  793. {PointSize[2point/3], Lighter[Brown], Point[{x9x, y9y, z9z}]}
  794. }, {}],
  795.  
  796. PlotRange->plotrange, AspectRatio->1, Axes->True, ImageSize->imagesize],
  797.  
  798. Graphics3D[{PointSize[point], Red, Point[Evaluate[f2p[T]][[1]]]}],
  799. Graphics3D[{PointSize[point], Blue, Point[Evaluate[f2p[T]][[2]]]}],
  800. Graphics3D[{PointSize[point], Green, Point[Evaluate[f2p[T]][[3]]]}],
  801. Graphics3D[{PointSize[point], Magenta, Point[Evaluate[f2p[T]][[4]]]}],
  802. Graphics3D[{PointSize[point], Cyan, Point[Evaluate[f2p[T]][[5]]]}],
  803. Graphics3D[{PointSize[point], Orange, Point[Evaluate[f2p[T]][[6]]]}],
  804. Graphics3D[{PointSize[point], Purple, Point[Evaluate[f2p[T]][[7]]]}],
  805. Graphics3D[{PointSize[point], Pink, Point[Evaluate[f2p[T]][[8]]]}],
  806. Graphics3D[{PointSize[point], Brown, Point[Evaluate[f2p[T]][[9]]]}],
  807.  
  808. ViewPoint->viewpoint]},
  809.  
  810. { },
  811. {s["t"->N[T]], sw[1/2]},
  812. { },
  813. {s["p1{x,y,z}"-> Evaluate[f2p[T][[1]]]], sw[1/2]},
  814. {s["v1{x,y,z}"-> Evaluate[f2v[T][[1]]]], sw[1/2]},
  815. {s["v1{total}"->{Evaluate[Chop@Norm[f2v[T][[1]]]]}], sw[1/2]},
  816. { },
  817. {s["p2{x,y,z}"-> Evaluate[f2p[T][[2]]]], sw[1/2]},
  818. {s["v2{x,y,z}"-> Evaluate[f2v[T][[2]]]], sw[1/2]},
  819. {s["v2{total}"->{Evaluate[Chop@Norm[f2v[T][[2]]]]}], sw[1/2]},
  820. { },
  821. {s["p3{x,y,z}"-> Evaluate[f2p[T][[3]]]], sw[1/2]},
  822. {s["v3{x,y,z}"-> Evaluate[f2v[T][[3]]]], sw[1/2]},
  823. {s["v3{total}"->{Evaluate[Chop@Norm[f2v[T][[3]]]]}], sw[1/2]},
  824. { },
  825. {s["p4{x,y,z}"-> Evaluate[f2p[T][[4]]]], sw[1/2]},
  826. {s["v4{x,y,z}"-> Evaluate[f2v[T][[4]]]], sw[1/2]},
  827. {s["v4{total}"->{Evaluate[Chop@Norm[f2v[T][[4]]]]}], sw[1/2]},
  828. { },
  829. {s["p5{x,y,z}"-> Evaluate[f2p[T][[5]]]], sw[1/2]},
  830. {s["v5{x,y,z}"-> Evaluate[f2v[T][[5]]]], sw[1/2]},
  831. {s["v5{total}"->{Evaluate[Chop@Norm[f2v[T][[5]]]]}], sw[1/2]},
  832. { },
  833. {s["p6{x,y,z}"-> Evaluate[f2p[T][[6]]]], sw[1/2]},
  834. {s["v6{x,y,z}"-> Evaluate[f2v[T][[6]]]], sw[1/2]},
  835. {s["v6{total}"->{Evaluate[Chop@Norm[f2v[T][[6]]]]}], sw[1/2]},
  836. { },
  837. {s["p7{x,y,z}"-> Evaluate[f2p[T][[7]]]], sw[1/2]},
  838. {s["v7{x,y,z}"-> Evaluate[f2v[T][[7]]]], sw[1/2]},
  839. {s["v7{total}"->{Evaluate[Chop@Norm[f2v[T][[7]]]]}], sw[1/2]},
  840. { },
  841. {s["p8{x,y,z}"-> Evaluate[f2p[T][[8]]]], sw[1/2]},
  842. {s["v8{x,y,z}"-> Evaluate[f2v[T][[8]]]], sw[1/2]},
  843. {s["v8{total}"->{Evaluate[Chop@Norm[f2v[T][[8]]]]}], sw[1/2]},
  844. { },
  845. {s["p9{x,y,z}"-> Evaluate[f2p[T][[9]]]], sw[1/2]},
  846. {s["v9{x,y,z}"-> Evaluate[f2v[T][[9]]]], sw[1/2]},
  847. {s["v9{total}"->{Evaluate[Chop@Norm[f2v[T][[9]]]]}], sw[1/2]},
  848. { },
  849. {s["ps{x,y,z}"-> swp[T]], sw[1/2]},
  850. {s["vs{x,y,z}"-> swp'[T]], sw[1/2]},
  851. {s["vs{total}"->{Chop@Norm[swp'[T]]}], sw[1/2]}
  852. }, Alignment->Left]]],
  853.  
  854. (* Zeitregler *)
  855.  
  856. {T, 0, tMax, tMax/5}]
  857.  
  858. (* Export als HTML Dokument *)
  859. (* Export["dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
  860. (* Export direkt als Bildsequenz *)
  861. (* ParallelDo[Export["dateiname" <> ToString[T] <> ".png", Rasterize[...] ], {T, 0, 10, 5}] *)
  862.  
  863.  
  864.  
  865.  
  866.  
  867.  
  868.  
  869.  
  870.  
RAW Paste Data