% N-body gravity simulation using random initial conditions.
% Positional units are AU, time units are Yr.
function random_stars
K = 20; tmax = 10; dt = 0.0005; L = 2;
x0 = r(K)*L/2; y0 = r(K)*L/2; z0 = r(K)*L/2;
vx0 = 0.3*r(K); vy0 = 0.3*r(K); vz0 = 0.3*r(K);
M = rand(1,K) + 1;
[t,x,y,z,vx,vy,vz] = gravity( dt,tmax,x0,y0,z0,vx0,vy0,vz0,M );
plot3(x,y,z); axis equal; axis vis3d; grid on; box on;
xlabel('x /AU'); ylabel('y /AU'); zlabel('z /AU'); 
xlim([-L,L]); ylim([-L,L]); zlim([-L,L]);
set( gcf, 'units','normalized','position',[0.05,0.05,0.85,0.85] );
print(gcf,'random stars.png','-dpng','-r300');

%%

%Generates N random numbers between -1 and 1
function n = r(K)
n = (-1+2*rand(1,K));

%End of code
