/sci/ - Science & Math

OP here again. I really want someone else to post some code!

Here is a cute little n-body gravity simulation that I know you will love:

clear all, close all; clc
figure('units','normalized','position',[.1 .1 .7 .7])

to = 0;
tf = 10;
dt = .01;
T = to:dt:tf;
n = 100;
a = 1e3*.2*dt;
r = 20;

Xold1 = [r*2*(rand(n,3) - .5)] + 50;
Xold2 = [r*2*(rand(n,3) - .5)] - 50;

Xold = [Xold1;Xold2;0 0 0];
Xc = Xold*[cos(a*pi/180) -sin(a*pi/180) 0;
sin(a*pi/180) cos(a*pi/180) 0;
0 0 1];

M = [1e2*(rand(2*n,1).^5); 1e7];

N = length(M);
dt2 = dt^2;
Xn = zeros(N,3);

for ti = 1:length(T)
F = zeros(N,3);
for n = 1:N
I = [1:n-1,n+1:N];
for i = 1:N-1
V = Xc(I(i),:) - Xc(n,:);
norm_V = norm(V);
F(n,:) = F(n,:) + (M(I(i))/(norm_V^3))*V;
end
end
Xn = 2*Xc - Xold + dt2*F;
Xold = Xc;
Xc = Xn;
plot3(Xn(:,1),Xn(:,2),Xn(:,3),'k*'), grid on
axis([-150 150 -150 150 -150 150]), getframe;
end