Zitat:
Also innerhalb des Codes einer Methode wüsste ich nicht, wie ich da die Einheiten mitschleifen soll, schließlich soll man ja nicht jede einzelne Zeile kommentieren.
Gannz einfach:
Code:
(* Yukterez Code, Gravity Simulator *)
(* Einstellungen *)
M1 = mSol; M2 = 2 mSol; M3 = 3 mSol; M4 = 4 mSol; (* Massen *)
P1x = Au; P2x = Au; P3x = -Au; P4x = -Au; (* Anfangspositionen *)
P1y = Au; P2y = -Au; P3y = Au; P4y = -Au;
P1z = 0; P2z = 0; P3z = 0; P4z = 0;
v1x = 0; v2x = 0; v3x = 0; v4x = 0; (* Anfangsgeschwindigkeiten *)
v1y = 0; v2y = 0; v3y = 0; v4y = 0;
v1z = 0; v2z = 0; v3z = 0; v4z = 0;
(* Eckdaten *)
G = 667384/10^16 m³/kg/sek²; (* Gravitationskonstante *)
Au = 149597870690 m; (* Astronomische Einheit *)
mSol = 1988*10^27 kg; (* Sonnenmasse *)
(* Dimensionen *)
kg=1; m=1; sek=1;
(* Container *)
Funktion[{{x1x_, y1y_, z1z_}, {x2x_, y2y_, z2z_}, {x3x_, y3y_,
z3z_}, {x4x_, y4y_, z4z_}}, {{vx1x_, vy1y_, vz1z_}, {vx2x_, vy2y_,
vz2z_}, {vx3x_, vy3y_, vz3z_}, {vx4x_, vy4y_, vz4z_}}, {m1_, m2_,
m3_, m4_}, T_, plotType : ("x" | "v"), plotOptions___] :=
Module[{nds, Tmax, funcToPlot},
(* DGL *)
nds = NDSolve[{
x1'[t] == vx1[t], y1'[t] == vy1[t], z1'[t] == vz1[t],
x2'[t] == vx2[t], y2'[t] == vy2[t], z2'[t] == vz2[t],
x3'[t] == vx3[t], y3'[t] == vy3[t], z3'[t] == vz3[t],
x4'[t] == vx4[t], y4'[t] == vy4[t], z4'[t] == vz4[t],
vx1'[t] == (G m2 (x2[t] - x1[t]))/
Sqrt[((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 + (z2[t] -
z1[t])^2)^3] + (G m3 (x3[t] - x1[t]))/
Sqrt[((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2 + (z3[t] -
z1[t])^2)^3] + (G m4 (x4[t] - x1[t]))/
Sqrt[((x4[t] - x1[t])^2 + (y4[t] - y1[t])^2 + (z4[t] -
z1[t])^2)^3],
vy1'[t] == (G m2 (y2[t] - y1[t]))/
Sqrt[((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 + (z2[t] -
z1[t])^2)^3] + (G m3 (y3[t] - y1[t]))/
Sqrt[((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2 + (z3[t] -
z1[t])^2)^3] + (G m4 (y4[t] - y1[t]))/
Sqrt[((x4[t] - x1[t])^2 + (y4[t] - y1[t])^2 + (z4[t] -
z1[t])^2)^3],
vz1'[t] == (G m2 (z2[t] - z1[t]))/
Sqrt[((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 + (z2[t] -
z1[t])^2)^3] + (G m3 (z3[t] - z1[t]))/
Sqrt[((x3[t] - x1[t])^2 + (y3[t] - y1[t])^2 + (z3[t] -
z1[t])^2)^3] + (G m4 (z4[t] - z1[t]))/
Sqrt[((x4[t] - x1[t])^2 + (y4[t] - y1[t])^2 + (z4[t] -
z1[t])^2)^3],
vx2'[t] == (G m1 (x1[t] - x2[t]))/
Sqrt[((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] -
z2[t])^2)^3] + (G m3 (x3[t] - x2[t]))/
Sqrt[((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2 + (z3[t] -
z2[t])^2)^3] + (G m4 (x4[t] - x2[t]))/
Sqrt[((x4[t] - x2[t])^2 + (y4[t] - y2[t])^2 + (z4[t] -
z2[t])^2)^3],
vy2'[t] == (G m1 (y1[t] - y2[t]))/
Sqrt[((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2 + (z1[t] -
z2[t])^2)^3] + (G m3 (y3[t] - y2[t]))/
Sqrt[((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2 + (z3[t] -
z2[t])^2)^3] + (G m4 (y4[t] - y2[t]))/
Sqrt[((x4[t] - x2[t])^2 + (y4[t] - y2[t])^2 + (z4[t] -
z2[t])^2)^3],
vz2'[t] == (G m1 (z1[t] - z2[t]))/
Sqrt[((x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 + (z2[t] -
z1[t])^2)^3] + (G m3 (z3[t] - z2[t]))/
Sqrt[((x3[t] - x2[t])^2 + (y3[t] - y2[t])^2 + (z3[t] -
z2[t])^2)^3] + (G m4 (z4[t] - z2[t]))/
Sqrt[((x4[t] - x2[t])^2 + (y4[t] - y2[t])^2 + (z4[t] -
z2[t])^2)^3],
vx3'[t] == (G m1 (x1[t] - x3[t]))/
Sqrt[((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2 + (z1[t] -
z3[t])^2)^3] + (G m2 (x2[t] - x3[t]))/
Sqrt[((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2 + (z2[t] -
z3[t])^2)^3] + (G m4 (x4[t] - x3[t]))/
Sqrt[((x4[t] - x3[t])^2 + (y4[t] - y3[t])^2 + (z4[t] -
z3[t])^2)^3],
vy3'[t] == (G m1 (y1[t] - y3[t]))/
Sqrt[((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2 + (z1[t] -
z3[t])^2)^3] + (G m2 (y2[t] - y3[t]))/
Sqrt[((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2 + (z2[t] -
z3[t])^2)^3] + (G m4 (y4[t] - y3[t]))/
Sqrt[((x4[t] - x3[t])^2 + (y4[t] - y3[t])^2 + (z4[t] -
z3[t])^2)^3],
vz3'[t] == (G m1 (z1[t] - z3[t]))/
Sqrt[((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2 + (z1[t] -
z3[t])^2)^3] + (G m2 (z2[t] - z3[t]))/
Sqrt[((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2 + (z2[t] -
z3[t])^2)^3] + (G m4 (z4[t] - z3[t]))/
Sqrt[((x4[t] - x3[t])^2 + (y4[t] - y3[t])^2 + (z4[t] -
z3[t])^2)^3],
vx4'[t] == (G m1 (x1[t] - x4[t]))/
Sqrt[((x1[t] - x4[t])^2 + (y1[t] - y4[t])^2 + (z1[t] -
z4[t])^2)^3] + (G m2 (x2[t] - x4[t]))/
Sqrt[((x2[t] - x4[t])^2 + (y2[t] - y4[t])^2 + (z2[t] -
z4[t])^2)^3] + (G m3 (x3[t] - x4[t]))/
Sqrt[((x3[t] - x4[t])^2 + (y3[t] - y4[t])^2 + (z3[t] -
z4[t])^2)^3],
vy4'[t] == (G m1 (y1[t] - y4[t]))/
Sqrt[((x1[t] - x4[t])^2 + (y1[t] - y4[t])^2 + (z1[t] -
z4[t])^2)^3] + (G m2 (y2[t] - y4[t]))/
Sqrt[((x2[t] - x4[t])^2 + (y2[t] - y4[t])^2 + (z2[t] -
z4[t])^2)^3] + (G m3 (y3[t] - y4[t]))/
Sqrt[((x3[t] - x4[t])^2 + (y3[t] - y4[t])^2 + (z3[t] -
z4[t])^2)^3],
vz4'[t] == (G m1 (z1[t] - z4[t]))/
Sqrt[((x1[t] - x4[t])^2 + (y1[t] - y4[t])^2 + (z1[t] -
z4[t])^2)^3] + (G m2 (z2[t] - z4[t]))/
Sqrt[((x2[t] - x4[t])^2 + (y2[t] - y4[t])^2 + (z2[t] -
z4[t])^2)^3] + (G m3 (z3[t] - z4[t]))/
Sqrt[((x3[t] - x4[t])^2 + (y3[t] - y4[t])^2 + (z3[t] -
z4[t])^2)^3],
x1[0] == x1x, y1[0] == y1y, z1[0] == z1z,
x2[0] == x2x, y2[0] == y2y, z2[0] == z2z,
x3[0] == x3x, y3[0] == y3y, z3[0] == z3z,
x4[0] == x4x, y4[0] == y4y, z4[0] == z4z,
vx1[0] == vx1x, vy1[0] == vy1y, vz1[0] == vz1z,
vx2[0] == vx2x, vy2[0] == vy2y, vz2[0] == vz2z,
vx3[0] == vx3x, vy3[0] == vy3y, vz3[0] == vz3z,
vx4[0] == vx4x, vy4[0] == vy4y, vz4[0] == vz4z},
{x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4,
vx1, vx2, vx3, vx4, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4},
{t, 0, T}];
(* Solve *)
If[Head[nds] =!= NDSolve, Tmax = nds[[1, 1, 2, 1, 1, 2]];
funcToPlot =
If[plotType ===
"x", {{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}, {x3[t],
y3[t], z3[t]}, {x4[t], y4[t], z4[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]}}] /. nds[[1]];
(* Plot *)
ParametricPlot3D[Evaluate[funcToPlot], {t, 0, Tmax},
PlotStyle -> {{Red, Thick}, {Blue, Thick}, {Green,
Thick}, {Orange, Thick}},
(* Plot Range *)
PlotRange -> {{-2 Au, 2 Au}, {-2 Au, 2 Au}, {-2 Au, 2 Au}},
AspectRatio -> 1, MaxRecursion -> ControlActive[3, 100],
Axes -> True, plotOptions], Text["Yukterez Mod."]]] // Quiet
(* Manipulate *)
plot = Manipulate[Show[Funktion[
{{P1x, P1y, P1z}, {P2x, P2y, P2z}, {P3x, P3y, P3z}, {P4x, P4y, P4z}},
{{v1x, v1y, v1z}, {v2x, v2y, v2z}, {v3x, v3y, v3z}, {v4x, v4y, v4z}},
{M1, M2, M3, M4},
T, xv,
ImageSize -> {440, 440}],
(* Startpositionen *)
Graphics3D[{
{Red, Point[{P1x, P1y, P1z}]},
{Blue, Point[{P2x, P2y, P2z}]},
{Green, Point[{P3x, P3y, P3z}]},
{Orange, Point[{P4x, P4y, P4z}]}}],
(* Perspektive *)
ViewPoint -> {0, 0, 2}],
{{xv, "x", "Position"}, {"x" -> "Position"}},
(* Zeitregler *)
{{T, 1*^4, "Time"}, 1*^4, 1*10^7},
ControlPlacement -> {Bottom}]
Export["plot.avi", plot, "FrameRate" -> 75]
Zumindest hoffend dass man hier trotz relativ unleserlichem .txt Format die Einheiten alle nachvollziehen kann,
