Three body Simulation

t1=0.003;parnum=4;G=0.1;R=80;ϵ=1;
m={1000,1000,1000,10};mass0=Table[m,{i,1,parnum}];mass=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];For[i=1,i<parnum+1,i++,mass[[i]]=Drop[mass0[[i]],{i}]];(*x0={10,-10};y0={0,0};z0={0,0};*)
x0=RandomReal[{-30,30},parnum];y0=RandomReal[{-30,30},parnum];z0=RandomReal[{-30,30},parnum];
x0[[parnum]]=0;y0[[parnum]]=0;y0[[parnum]]=0;
xtrace=List[{x0[[1]],y0[[1]],z0[[1]]}];
(*vx0={0,0};vy0={1,-1};vz0={0,0};*)vx0=1RandomReal[{-1,1},parnum];vy0=1RandomReal[{-1,1},parnum];vz0=1RandomReal[{-1,1},parnum];
vx0[[parnum]]=0;vy0[[parnum]]=0;vz0[[parnum]]=0;

vx0[[1]]=-((m*vx0)[[2;;parnum]]//Total)/m[[1]];
vy0[[1]]=-((m*vy0)[[2;;parnum]]//Total)/m[[1]];
vz0[[1]]=-((m*vz0)[[2;;parnum]]//Total)/m[[1]];
While[True,
x0=x0+vx0*t1;y0=y0+vy0*t1;
z0=z0+vz0*t1;

(*xtrace=Append[xtrace,{x0[[1]],y0[[1]],z0[[1]]}];*)

dx1=Table[Table[x0[[j]],{i,1,parnum}],{j,1,parnum}];dx2=Table[x0,{i,1,parnum}];dx3=dx2-dx1;
dx=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];
dy1=Table[Table[y0[[j]],{i,1,parnum}],{j,1,parnum}];dy2=Table[y0,{i,1,parnum}];dy3=dy2-dy1;
dy=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];
dz1=Table[Table[z0[[j]],{i,1,parnum}],{j,1,parnum}];dz2=Table[z0,{i,1,parnum}];dz3=dz2-dz1;
dz=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];


For[i=1,i<parnum+1,i++,dx[[i]]=Drop[dx3[[i]],{i}]];
For[i=1,i<parnum+1,i++,dy[[i]]=Drop[dy3[[i]],{i}]];
For[i=1,i<parnum+1,i++,dz[[i]]=Drop[dz3[[i]],{i}]];


disx=dx*dx;
disy=dy*dy;
disz=dz*dz;
d=(disx+disy+disz)^(3/2)+ϵ;

forcex=G*Total[mass*dx/d,{2}];forcey=G*Total[mass*dy/d,{2}];forcez=G*Total[mass*dz/d,{2}];
ax=forcex;ay=forcey;az=forcez;

vx0=vx0+ax*t1;vy0=vy0+ay*t1;vz0=vz0+az*t1;

Pause[0];
figure=Table[{Hue[0.5],Ball[{x0[[i]],y0[[i]],z0[[i]]},(m[[i]])^(1/6)]},{i,1,parnum}];fig=Graphics3D[figure(*,PlotRange->{{-50,50},{-50,50},{-50,50}}*),PlotRange->{{Total[m*x0]/Total[m]-R,Total[m*x0]/Total[m]+R},{Total[m*y0]/Total[m]-R,Total[m*y0]/Total[m]+R},{Total[m*z0]/Total[m]-R,Total[m*z0]/Total[m]+R}},Boxed->True,(*PlotRange->{{x0[[parnum]]-R,x0[[parnum]]+R},{y0[[parnum]]-R,y0[[parnum]]+R},{z0[[parnum]]-R,z0[[parnum]]+R}},*)ImageSize->Full];]

Dynamic[fig]

GIF Maker

t1=0.1;parnum=4;G=0.1;R=80;ϵ=1;
m={1000,1000,1000,10};mass0=Table[m,{i,1,parnum}];mass=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];For[i=1,i<parnum+1,i++,mass[[i]]=Drop[mass0[[i]],{i}]];(*x0={10,-10};y0={0,0};z0={0,0};*)
x0=RandomReal[{-30,30},parnum];y0=RandomReal[{-30,30},parnum];z0=RandomReal[{-30,30},parnum];
x0[[parnum]]=0;y0[[parnum]]=0;y0[[parnum]]=0;
xtrace=List[{x0[[1]],y0[[1]],z0[[1]]}];
(*vx0={0,0};vy0={1,-1};vz0={0,0};*)vx0=1RandomReal[{-1,1},parnum];vy0=1RandomReal[{-1,1},parnum];vz0=1RandomReal[{-1,1},parnum];
vx0[[parnum]]=0;vy0[[parnum]]=0;vz0[[parnum]]=0;

vx0[[1]]=-((m*vx0)[[2;;parnum]]//Total)/m[[1]];
vy0[[1]]=-((m*vy0)[[2;;parnum]]//Total)/m[[1]];
vz0[[1]]=-((m*vz0)[[2;;parnum]]//Total)/m[[1]];
threebody=Table[0,{i,1,1000}];
picture=0;
While[True,
x0=x0+vx0*t1;y0=y0+vy0*t1;
z0=z0+vz0*t1;

(*xtrace=Append[xtrace,{x0[[1]],y0[[1]],z0[[1]]}];*)

dx1=Table[Table[x0[[j]],{i,1,parnum}],{j,1,parnum}];dx2=Table[x0,{i,1,parnum}];dx3=dx2-dx1;
dx=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];
dy1=Table[Table[y0[[j]],{i,1,parnum}],{j,1,parnum}];dy2=Table[y0,{i,1,parnum}];dy3=dy2-dy1;
dy=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];
dz1=Table[Table[z0[[j]],{i,1,parnum}],{j,1,parnum}];dz2=Table[z0,{i,1,parnum}];dz3=dz2-dz1;
dz=Table[Table[0,{i,1,parnum-1}],{j,1,parnum}];


For[i=1,i<parnum+1,i++,dx[[i]]=Drop[dx3[[i]],{i}]];
For[i=1,i<parnum+1,i++,dy[[i]]=Drop[dy3[[i]],{i}]];
For[i=1,i<parnum+1,i++,dz[[i]]=Drop[dz3[[i]],{i}]];


disx=dx*dx;
disy=dy*dy;
disz=dz*dz;
d=(disx+disy+disz)^(3/2)+ϵ;

forcex=G*Total[mass*dx/d,{2}];forcey=G*Total[mass*dy/d,{2}];forcez=G*Total[mass*dz/d,{2}];
ax=forcex;ay=forcey;az=forcez;

vx0=vx0+ax*t1;vy0=vy0+ay*t1;vz0=vz0+az*t1;

Pause[0.01];
figure=Table[{Hue[0.5],Ball[{x0[[i]],y0[[i]],z0[[i]]},(m[[i]])^(1/6)]},{i,1,parnum}];fig=Graphics3D[figure(*,PlotRange->{{-50,50},{-50,50},{-50,50}}*),PlotRange->{{Total[m*x0]/Total[m]-R,Total[m*x0]/Total[m]+R},{Total[m*y0]/Total[m]-R,Total[m*y0]/Total[m]+R},{Total[m*z0]/Total[m]-R,Total[m*z0]/Total[m]+R}},Boxed->True,(*PlotRange->{{x0[[parnum]]-R,x0[[parnum]]+R},{y0[[parnum]]-R,y0[[parnum]]+R},{z0[[parnum]]-R,z0[[parnum]]+R}},*)ImageSize->Full];
picture=picture+1;
threebody[[picture]]=fig;
If[picture>999,Break[]];]

Dynamic[fig]

threebody

Export["threebody1.gif",threebody,"DisplayDurations"->0.02,"AnimationRepetitions"->Infinity]