# Gravity & Charge, 2 Body Simulator

Feb 16th, 2019 (edited)
55
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. (* |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| *)
2. (* ||| Mathematica Syntax || yukterez.net || 2 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 = 24 sek;
21. tMax = Min[Tmax, plunge];
22. trail = 12 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 = 1000000000 kg;
43. q1 = 0;
44.
45. x1x = 1/2 m;
46. y1y = 0 m;
47. z1z = 0 m;
48.
49. v1x = 0 m/sek;
50. v1y = 0 m/sek;
51. v1z = Sqrt[G m2] Sqrt[1/2];
52.
53. (* Körper 2 *)
54.
55. m2 = m1/2;
56. q2 = 0 Amp sek;
57.
58. x2x = -1/2 m;
59. y2y = 0 m;
60. z2z = 0 m;
61.
62. v2x = 0 m/sek;
63. v2y = 0 m/sek;
64. v2z = -Sqrt[G m1] Sqrt[1/2];
65.
66. (* Differentialgleichung *)
67.
68. nds=NDSolve[{
69.
70. x1'[t] == vx1[t], y1'[t] == vy1[t], z1'[t] == vz1[t],
71. x2'[t] == vx2[t], y2'[t] == vy2[t], z2'[t] == vz2[t],
72.
73. vx1'[t] ==
74. (G m2 (x2[t]-x1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
75. If[q1 == 0, 0,
76. (-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]]+
77. Λ*c^2*x1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
78.
79. vy1'[t] ==
80. (G m2 (y2[t]-y1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
81. If[q1 == 0, 0,
82. (-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]]+
83. Λ*c^2*y1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
84.
85. vz1'[t] ==
86. (G m2 (z2[t]-z1[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
87. If[q1 == 0, 0,
88. (-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]]+
89. Λ*c^2*z1[t]^2/Sqrt[x1[t]^2+y1[t]^2+z1[t]^2],
90.
91. vx2'[t] ==
92. (G m1 (x1[t]-x2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
93. If[q2 == 0, 0,
94. (-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]]+
95. Λ*c^2*x2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
96.
97. vy2'[t] ==
98. (G m1 (y1[t]-y2[t]))/Sqrt[((x1[t]-x2[t])^2+(y1[t]-y2[t])^2+(z1[t]-z2[t])^2)^3]+
99. If[q2 == 0, 0,
100. (-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]]+
101. Λ*c^2*y2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
102.
103. vz2'[t] ==
104. (G m1 (z1[t]-z2[t]))/Sqrt[((x2[t]-x1[t])^2+(y2[t]-y1[t])^2+(z2[t]-z1[t])^2)^3]+
105. If[q2 == 0, 0,
106. (-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]]+
107. Λ*c^2*z2[t]^2/Sqrt[x2[t]^2+y2[t]^2+z2[t]^2],
108.
109. x1[0] == x1x, y1[0] == y1y, z1[0] == z1z,
110. x2[0] == x2x, y2[0] == y2y, z2[0] == z2z,
111.
112. vx1[0] == v1x, vy1[0] == v1y, vz1[0] == v1z,
113. vx2[0] == v2x, vy2[0] == v2y, vz2[0] == v2z},
114.
115. {x1, x2, y1, y2, z1, z2,
116. vx1, vx2, vy1, vy2, vz1, vz2},
117.
118. {t, 0, Tmax},
119.
120. WorkingPrecision-> wp,
121. MaxSteps-> Infinity,
122. Method-> mta,
123. InterpolationOrder-> All,
124. StepMonitor :> (laststep=plunge; plunge=t;
125. stepsize=plunge-laststep;), Method->{"EventLocator",
126. "Event" :> (If[stepsize<1*^-4, 0, 1])}];
127.
128. (* Position, Geschwindigkeit *)
129.
130. f2p[t_]={{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}}/.nds[[1]];
131. f2v[t_]={{vx1[t], vy1[t], vz1[t]}, {vx2[t], vy2[t], vz2[t]}}/.nds[[1]];
132. swp[t_]=(m1 Evaluate[f2p[t][[1]]]+m2 Evaluate[f2p[t][[2]]])/(m1+m2);
133.
134. (* Formatierung *)
135.
136. s[text_]=Style[text, FontSize->11];
137. sw[text_]=Style[text, White, FontSize->11];
138. colorfunc[n_]=Function[{x, y, z, t},
139. Hue[0, n, 0.5,
140. If[Tmax<0, Max[Min[(+T+(-t+trail))/trail, 1], 0],
141. Max[Min[(-T+(t+trail))/trail, 1], 0]]]];
142.
143. (* Animation *)
144.
145. Do[Print[Rasterize[
146. Grid[{{
147. Show[
148.
149. If[T == 0, {},
150.
151. ParametricPlot3D[Evaluate[f2p[t]],
152. {t, Max[0, T-trail], T},
153.
154. PlotStyle->{
155. {Thickness[thk], Red},
156. {Thickness[thk], Blue}},
157.
158. PlotRange->plotrange, AspectRatio->1, MaxRecursion->15, Axes->True, ImageSize->imagesize]],
159.
160. Graphics3D[
161. If[startpos==1, {
162. {PointSize[2point/3], Lighter[Red], Point[{x1x, y1y, z1z}]},
163. {PointSize[2point/3], Lighter[Blue],Point[{x2x, y2y, z2z}]}
164. }, {}],
165.
166. PlotRange->plotrange, AspectRatio->1, Axes->True, ImageSize->imagesize],
167.
168. Graphics3D[{PointSize[point], Red, Point[Evaluate[f2p[T]][[1]]]}],
169. Graphics3D[{PointSize[point], Blue, Point[Evaluate[f2p[T]][[2]]]}],
170.
171. ViewPoint->viewpoint]},
172.
173. { },
174. {s["t"->N[T]], sw[1/2]},
175. { },
176. {s["p1{x,y,z}"-> Evaluate[f2p[T][[1]]]], sw[1/2]},
177. {s["v1{x,y,z}"-> Evaluate[f2v[T][[1]]]], sw[1/2]},
178. {s["v1{total}"->{Evaluate[Chop@Norm[f2v[T][[1]]]]}], sw[1/2]},
179. { },
180. {s["p2{x,y,z}"-> Evaluate[f2p[T][[2]]]], sw[1/2]},
181. {s["v2{x,y,z}"-> Evaluate[f2v[T][[2]]]], sw[1/2]},
182. {s["v2{total}"->{Evaluate[Chop@Norm[f2v[T][[2]]]]}], sw[1/2]},
183. { },
184. {s["ps{x,y,z}"-> swp[T]], sw[1/2]},
185. {s["vs{x,y,z}"-> swp'[T]], sw[1/2]},
186. {s["vs{total}"->{Chop@Norm[swp'[T]]}], sw[1/2]}
187. }, Alignment->Left]]],
188.
189. (* Zeitregler *)
190.
191. {T, 0, tMax, tMax/5}]
192.
193. (* Export als HTML Dokument *)
194. (* Export["dateiname.html", EvaluationNotebook[], "GraphicsOutput" -> "PNG"] *)
195. (* Export direkt als Bildsequenz *)
196. (* ParallelDo[Export["dateiname" <> ToString[T] <> ".png", Rasterize[...] ], {T, 0, 10, 5}] *)
197.
198.
199.
200.
RAW Paste Data