# 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