%what can be changed to fit different computers:
%clockmax,dt,step
clear;
figure(1);
clf;
G=1;
N=500;%the whole body consist of 500 bodies
k=2500;%a parameter
clockmax=5000;% that is too small, not acurate but fast..recommand:>2000
dt=1/2000 ;% cannot be larger, so this program cannot run faster
time=clockmax*dt;
r0=1;
r1=N^(1/3)*r0;% so N bodies have the same volume as the whole body
L=3*r1;
m0=1000;
W=sqrt(G*m0*N/r1^3);% W must be large enough so that it can be seen, so m0 have to be large
%r(1)=0;a(1)=0;b(1)=0;x(1)=0;y(1)=0;z(1)=0;u(1)=0;v(1)=0;w(1)=0;m(1)=m0;
%it's first designed as a massive center, but finnally I decided not use this
clock=1;record=1;
step=5;%based on computer's memory,too small will make computer very slow.
while clock<=N % to make them in random place, but have to be symatric so that the mass center is the origin.
x1=(rand-0.5)*r1*2;
y1=(rand-0.5)*r1*2;
z1=(rand-0.5)*r1*2;
if x1^2+y1^2+z1^2<=r1^2;%in the sphere
x(clock)=x1;
y(clock)=y1;
z(clock)=z1;
x(clock+1)=-x1;
y(clock+1)=-y1;
z(clock+1)=-z1;
%r(clock)=sqrt(x1^2+y1^2+z1^2); useless
%r(clock+1)=sqrt(x1^2+y1^2+z1^2);
clock=clock+2;
end%else nothing happend.
end
%save them, new part.
xsave=x;
ysave=y;
zsave=z;
for i=2:N%I first use spherical coordinate, but that means center have much more bodies, much higher dencity
%r(i)=r1*rand;
%a(i)=2*pi*rand;
%b(i)=2*pi*rand;
%x(i)=r(i)*cos(b(i))*cos(a(i));
%y(i)=r(i)*cos(b(i))*sin(a(i));
%z(i)=r(i)*sin(b(i));
u(i)=-W*y(i);%keep this part to determine velocity and mass
v(i)=W*x(i);
w(i)=0;
m(i)=m0;
end
%2*2 subplot, x-y,y-z,3D and a graph to show the radius.
set(gcf,'double','on');
subplot(2,2,1),h1=plot(x,y,'ro');
axis([-L,L,-L,L]);
title('x-y')
xlabel('x'),ylabel('y')
set(gcf,'double','on');
subplot(2,2,2),h2=plot(y,z,'ro');
axis([-L,L,-L,L])
title('z-y')
xlabel('z'),ylabel('y')
set(gcf,'double','on');
subplot(2,2,4),h=plot3(x,y,z,'ro');
axis([-L,L,-L,L,-L,L]);
title('3D plot')
%I didn't set a hook,but use "hold on" and plot many dot on the graph tohave a line
subplot(2,2,3),hold on
axis([0,time,.5*r1,3*r1])
xlabel('time'),ylabel('maxium on each axis')
plot (0,max(abs(x)),'b.')
plot (0,max(abs(x)),'y.')
plot (0,max(abs(x)),'g.')
legend('X','Y','Z')
drawnow
for clock=1:clockmax
t(clock)=clock*dt;
for i=1:N
for j=i+1:N % a way to make computer compute less
r=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2+(z(i)-z(j))^2);
if(r>2*r0)% not touched
du=G*dt*m0*(x(j)-x(i))/r^3;
dv=G*dt*m0*(y(j)-y(i))/r^3;
dw=G*dt*m0*(z(j)-z(i))/r^3;
u(i)=u(i)+du;
v(i)=v(i)+dv;
w(i)=w(i)+dw;
u(j)=u(j)-du;
v(j)=v(j)-dv;
w(j)=w(j)-dw;
end
if(r<2*r0)% touched
du=-dt*k*(2*r0-r)*(x(j)-x(i))/r;
dv=-dt*k*(2*r0-r)*(y(j)-y(i))/r;
dw=-dt*k*(2*r0-r)*(z(j)-z(i))/r;
u(i)=u(i)+du;
v(i)=v(i)+dv;
w(i)=w(i)+dw;
u(j)=u(j)-du;
v(j)=v(j)-dv;
w(j)=w(j)-dw;
end
end
end
for i=1:N
x(i)=x(i)+dt*u(i);
y(i)=y(i)+dt*v(i);
z(i)=z(i)+dt*w(i);
end
figure(1)
set(h,'xdata',x,'ydata',y,'zdata',z);
set(h1,'xdata',x,'ydata',y);
set(h2,'xdata',z,'ydata',y);
drawnow
%find radius on each axis
xm(clock)=max(abs(x));
ym(clock)=max(abs(y));
zm(clock)=max(abs(z));
subplot(2,2,3)%here I use avrage so that it won't change so dramatically.
plot (t(clock),mean(xm),'b.')
plot (t(clock),mean(ym),'y.')
plot (t(clock),mean(zm),'g.')
drawnow
%new part, to save and run it faster later
%careful! use a lot of memory!!
if (step*record<clock)
xsave=[xsave;x];
ysave=[ysave;y];
zsave=[zsave;z];
record=record+1
end
end
pause
for i=1:record
figure(2)
subplot(2,2,1)
plot3(xsave(i,:),ysave(i,:),zsave(i,:),'ro')
axis([-L,L,-L,L,-L,L]);
subplot(2,2,2)
plot(xsave(i,:),zsave(i,:),'ro')
axis([-L,L,-L,L]);
title('x-z')
xlabel('x'),ylabel('z')
subplot(2,2,3)
plot(ysave(i,:),zsave(i,:),'ro')
axis([-L,L,-L,L]);
title('z-y')
xlabel('z'),ylabel('y')
subplot(2,2,4)
plot(xsave(i,:),ysave(i,:),'ro')
axis([-L,L,-L,L]);
title('x-y')
xlabel('x'),ylabel('y')
pause (0.1)
end