(* Plot Optionen *)
Tmax=3yr; trail=yr/24; r=1.2 Au; point=0.01; font=11;

(* Sonne *)
m1=mSol;
q1=0; 
x1x=+4.767411813724087*^+5 km;
y1y=-4.881392129931528*^+4 km;
z1z=-2.184128405553382*^+4 km;
v1x=+5.632851183659290*^-3 km/sek;
v1y=+1.020293715934504*^-2 km/sek;
v1z=-1.481762989445804*^-4 km/sek;

(* Merkur *)
m2= 3.302*^23 kg;
q2=0;
x2x=+5.357892741028007*^+7 km;
y2y=-1.147221840127503*^+6 km;
z2z=-4.983469907397524*^+6 km;
v2x=-8.396224055586826*^+0 km/sek;
v2y=+5.089273470880620*^+1 km/sek;
v2z=+4.928244620166743*^+0 km/sek;

(* Venus *)
m3=4.8685*^24 kg;
q3=0;
x3x=-3.828375486384919*^+7 km;
y3y=+1.001776253505679*^+8 km;
z3z=+3.588675833902176*^+6 km;
v3x=-3.277640635284602*^+1 km/sek;
v3y=-1.281675204600856*^+1 km/sek;
v3z=+1.715864558132450*^+0 km/sek;

(* Erde *)
m4= 5.97219*^24kg; 
q4=0;
x4x=-1.436305243395423*^+8 km;
y4y=-4.046527983363845*^+7 km;
z4z=-2.022577828499544*^+4 km;
v4x=+7.558464309586846*^+0 km/sek;
v4y=-2.877403274780793*^+1 km/sek;
v4z=-2.203233704649544*^-4 km/sek;

(* Schweif *)
Tmin=yr/12;

(* Eckdaten *)
G=667384/10^16; (*Gravitationskonstante*) 
ε0=8.854187817*^-12;  (* Feldkonstante *) 
Au=149597870690; (*Astronomische Einheit*) 
mSol=1988*10^27; (*Sonnenmasse*) 
kg=1; m=1; sek=1;km=1000m;  yr=36525*24*36sek;

(* 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)]))+


((-q1*q2/(4Pi e0 )/m1 (x2[t]-x1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q1*q3/(4Pi e0 )/m1(x3[t]-x1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((-q1*q4/(4Pi e0 )/m1 (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)]))+


((-q1*q2/(4Pi e0 )/m1(y2[t]-y1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q1*q3/(4Pi e0 )/m1(y3[t]-y1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((-q1*q4/(4Pi e0 )/m1 (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)]))+


((-q1*q2/(4Pi e0 )/m1(z2[t]-z1[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q1*q3/(4Pi e0 )/m1(z3[t]-z1[t]))/(Sqrt[((((((x3[t]-x1[t]))^2)+(((y3[t]-y1[t]))^2)+(((z3[t]-z1[t]))^2)))^3)]))+((-q1*q4/(4Pi e0 )/m1 (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)]))+


((-q2*q1/(4Pi e0 )/m2  (x1[t]-x2[t]))/(Sqrt[((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)+(((z1[t]-z2[t]))^2)))^3)]))+((-q2*q3/(4Pi e0 )/m2  (x3[t]-x2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((-q2*q4/(4Pi e0 )/m2  (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)]))+


((-q2*q1/(4Pi e0 )/m2  (y1[t]-y2[t]))/(Sqrt[((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)+(((z1[t]-z2[t]))^2)))^3)]))+((-q2*q3/(4Pi e0 )/m2(y3[t]-y2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((-q2*q4/(4Pi e0 )/m2 (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)]))+


((-q2*q1/(4Pi e0 )/m2  (z1[t]-z2[t]))/(Sqrt[((((((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2)+(((z2[t]-z1[t]))^2)))^3)]))+((-q2*q3/(4Pi e0 )/m2 (z3[t]-z2[t]))/(Sqrt[((((((x3[t]-x2[t]))^2)+(((y3[t]-y2[t]))^2)+(((z3[t]-z2[t]))^2)))^3)]))+((-q2*q4/(4Pi e0 )/m2 (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)]))+


((-q3*q1/(4Pi e0 )/m3  (x1[t]-x3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((-q3*q2/(4Pi e0 )/m3 (x2[t]-x3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((-q3*q4/(4Pi e0 )/m3 (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)]))+


((-q3*q1/(4Pi e0 )/m3   (y1[t]-y3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((-q3*q2/(4Pi e0 )/m3   (y2[t]-y3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((-q3*q4/(4Pi e0 )/m3 (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)]))+


((-q3*q1/(4Pi e0 )/m3  (z1[t]-z3[t]))/(Sqrt[((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)+(((z1[t]-z3[t]))^2)))^3)]))+((-q3*q2/(4Pi e0 )/m3  (z2[t]-z3[t]))/(Sqrt[((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)+(((z2[t]-z3[t]))^2)))^3)]))+((-q3*q4/(4Pi e0 )/m3 (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)]))+


((-q4*q1/(4Pi e0 )/m4  (x1[t]-x4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((-q4*q2/(4Pi e0 )/m4   (x2[t]-x4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((-q4*q3/(4Pi e0 )/m4  (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)]))+


((-q4*q1/(4Pi e0 )/m4   (y1[t]-y4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((-q4*q2/(4Pi e0 )/m4   (y2[t]-y4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((-q4*q3/(4Pi e0 )/m4 (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)]))+


((-q4*q1/(4Pi e0 )/m4  (z1[t]-z4[t]))/(Sqrt[((((((x1[t]-x4[t]))^2)+(((y1[t]-y4[t]))^2)+(((z1[t]-z4[t]))^2)))^3)]))+((-q4*q2/(4Pi e0 )/m4  (z2[t]-z4[t]))/(Sqrt[((((((x2[t]-x4[t]))^2)+(((y2[t]-y4[t]))^2)+(((z2[t]-z4[t]))^2)))^3)]))+((-q4*q3/(4Pi e0 )/m4 (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]==v1x,vy1[0]==v1y,vz1[0]==v1z,
vx2[0]==v2x,vy2[0]==v2y,vz2[0]==v2z,
vx3[0]==v3x,vy3[0]==v3y,vz3[0]==v3z,
vx4[0]==v4x,vy4[0]==v4y,vz4[0]==v4z},

{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,Tmax}, MaxSteps->Infinity];

(* Position, Geschwindigkeit *)
f2p={{x1[t],y1[t],z1[t]},{x2[t],y2[t],z2[t]},{x3[t],y3[t],z3[t]},{x4[t],y4[t],z4[t]}}/.nds[[1]];
f2v={{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]];

(* Formatierung *)
s[text_]:=Style[text,FontSize->font]

(* Animation *)
Do[Print[Rasterize[
Grid[{{
Show[

ParametricPlot3D[Evaluate[f2p],
{t,T-trail,T},
PlotStyle->{{Red,Thick},{Blue,Thick},{Green,Thick},{Orange,Thick}},
(*Plot Range*)
PlotRange->{{-r,r},{-r,r},{-r, r}},AspectRatio->1,MaxRecursion->ControlActive[3,15],Axes->True, ImageSize->440],

Graphics3D[{
{PointSize[point/2],Lighter[Red],Point[{x1x,y1y,z1z}]},
{PointSize[point/2],Lighter[Blue],Point[{x2x,y2y,z2z}]},
{PointSize[point/2],Lighter[Green],Point[{x3x,y3y,z3z}]},
{PointSize[point/2],Lighter[Orange],Point[{x4x,y4y,z4z}]}}],

Graphics3D[{PointSize[point], Red,Point[With[{t=T},Evaluate[f2p]][[1]]]}],
Graphics3D[{PointSize[point], Blue,Point[With[{t=T},Evaluate[f2p]][[2]]]}],
Graphics3D[{PointSize[point], Green,Point[With[{t=T},Evaluate[f2p]][[3]]]}],
Graphics3D[{PointSize[point], Orange,Point[With[{t=T},Evaluate[f2p]][[4]]]}],

ViewPoint->{4,4,4}]},

{ },
{s["p1{x,y,z}"->With[{t=T},Evaluate[f2p]][[1]]/Au]},
{s["v1{x,y,z}"->With[{t=T},Evaluate[f2v]][[1]]]},
{ },
{s["p2{x,y,z}"->With[{t=T},Evaluate[f2p]][[2]]/Au]},
{s["v2{x,y,z}"->With[{t=T},Evaluate[f2v]][[2]]]},
{ },
{s["p3{x,y,z}"->With[{t=T},Evaluate[f2p]][[3]]/Au]},
{s["v3{x,y,z}"->With[{t=T},Evaluate[f2v]][[3]]]},
{ },
{s["p4{x,y,z}"->With[{t=T},Evaluate[f2p]][[4]]/Au]},
{s["v4{x,y,z}"->With[{t=T},Evaluate[f2v]][[4]]]}
}, Alignment->Left]]],

(* Zeitregler *)
{T,trail,Tmax,trail/10}]