Iris
Posts:
8
Registered:
11/27/12
|
|
Creating a movie for n-body simulations
Posted:
Nov 27, 2012 3:38 PM
|
|
Hello guys! I've done a n-body simulation program in matlab, now is defined for our solar system (till Saturn) and i've added a solar twin in certain position to see what happens. Now i need a movie from it, and my idea is to do a scatter for the present position and a line till the beginning of the simulation (following the path of the particles). This is all the code, i wish i could do the upload of the .m files (hope that runs ok)... Thank you for your help :)
if true
%code
N = 8;
yr = 365.25*24*3600;
T = input('Define the interval of time:');
fprintf('T: %d\n', T)
G = 6.67300e-20;%km^3*kg^-1*s^-2
R = [0 0 0; %--> sun
6.7e7 0 0; %--> mercury
1.08e8 0 0; % --> venus
1.5e8 0 0; % --> earth
2.28e8 0 0; % --> mars
7.78e8 0 0; % --> jupiter
1.43e9 0 0; % --> saturn
-1e9 1.5e9 1e8]; % --> 2nd sun
R = [0 0 0; %--> sun
6.7e7 0 0; %--> mercury
1.08e8 0 0; % --> venus
1.5e8 0 0; % --> earth
2.28e8 0 0; % --> mars
% 7.761173e8 0 0; %--> callisto
% 7.7773291e8 0 0; % --> europa
% 7.7775782e8 0 0; % --> io
% 7.776929600 0 0; % --> ganymede
7.78e8 0 0; % --> jupiter
1.43e9 0 0; % --> saturn
-1e9 1.5e9 1e8]; % --> 2nd sun
V = [0 0 0; %--> sun
0 47.9 0; %--> mercury
0 35 0; % --> venus
0 30 0; % --> earth
0 24 0; % --> mars
% 0 (8.2 + vj) 0; %--> callisto
% 0 (13.74 + vj) 0; % --> europa
% 0 (17 + vj) 0; % --> io
% 0 (10.9 +vj) 0; % --> ganymede
0 13 0; % --> jupiter
0 9.7 0; % --> saturn
10 -20 0]; % --> 2nd sun
m = [2e30 3.3e23 4.9e24 6e24 6.4e23 1.9e27 5.7e26 2e30]; % 10.7e22 4.7e22 8.9e22 14.8e22 1.9e27];
M = sum(m);
% computing the center of mass
%dt = (1/2)min(ri/vi;vi/ai)
dt = yr/36500;
for i = 1:N
Rab(i) = (R(i,1).^2 + R(i,2).^2 + R(i,3).^2);
Vab(i) = (V(i,1).^2 + V(i,2).^2 + V(i,3).^2);
end
Min = ((1/2)*min(Rab./Vab))*0.0000000001;
if dt > Min
dt = Min;
else
dt = dt;
end
for j = 1:3
for i = 1:N
cm(i,j) = m(i)*R(i,j);
end
end
CM = (1/M)*sum(cm);
for i = 1:N
for j = 1:3
R(i,j) = R(i,j) - CM(j); %--> new possitons for R with (0,0,0) in CM
end
end
% velocity
for j = 1:3
for i = 1:N
cmv(i,j) = m(i)*V(i,j);
end
end
CMv = (1/M)*sum(cmv);
for i = 1:N
for j = 1:3
V(i,j) = V(i,j) - CMv(j); %--> new velocities for R with (0,0,0) in CM
end
end
Rr = reshape(R.',1,3*N).';
Rd = reshape(Rr,1,3*N);
Ra(1, :) = Rd;
Vr = reshape(V.',1,3*N).';
Vd = reshape(Vr,1,3*N);
Va(1, :) = Vd;
o = 1;
t = dt;
while t < T
o = o + 1;
count1 = 0;
for i = 1:N
count1 = count1 + 1;
ac = zeros(N,3);
n = [1:N];
n(n==i) = [];
count2 = 0;
pt = zeros(N-1,1);
for v = n(1:(N-1))
count2 = count2 + 1;
dist = ((R(i,1) - R(v,1))^2 + (R(i,2) - R(v,2))^2 + (R(i,3) - R(v,3))^2)^(1/2);
c = dist^3;
pa = m(i)*m(v)/dist;
pt(count2,1) = pa;
for j = 1:3
b(i,j) = (R(i,j) - R(v,j))/c;
ac(v,j) = -G*m(v).*b(i,j);
end
end
pb = sum(pt);
pbt(count1,1) = pb;
accelerations(i,:) = sum(ac);
end
pc = sum(pbt);
Ar = reshape(accelerations.',1,3*N).';
Ad = reshape(Ar,1,3*N);
for i = 1:N
for j = 1:3
V(i,j) = V(i,j) + accelerations(i,j)*dt;
end
end
Vr = reshape(V.',1,3*N).';
Vd = reshape(Vr,1,3*N);
for i = 1:N
for j = 1:3
R(i,j) = R(i,j) + V(i,j)*dt + (1/2)*accelerations(i,j)*(dt)^2;
end
end
Rr = reshape(R.',1,3*N).';
Rd = reshape(Rr,1,3*N);
for i = 1:N
Rab(:, i) = (Rd(:,(3*i-2))^2 + Rd(:,(3*i-1))^2 + Rd(:,(3*i))^2)^(1/2);
Vab(:, i) = (Vd(:,(3*i-2))^2 + Vd(:,(3*i-1))^2 + Vd(:,(3*i))^2)^(1/2);
Aab(:, i) = (Ad(:,(3*i-2))^2 + Ad(:,(3*i-1))^2 + Ad(:,(3*i))^2)^(1/2);
end
Min1 = (min(Rab./Vab))*0.001;
Min2 = (min(Vab./Aab))*0.001;
MIN = min([Min1 Min2]);
dt = MIN;
p(o, :) = pc;
Ra(o, :) = Rd;
Va(o, :) = Vd;
t = t + dt;
end
i_t = size(Ra);
i_t = i_t(1);
hold on
color = ['b' 'c' 'g' 'y' 'r' 'm' 'b' 'c' 'g' 'y'];
%area = [pi*(695500000*39.37001*72)^2 pi*(2440000*39.37001*72)^2 pi*(6502000*39.37001*72)^2 pi*(6378000*39.37001*72)^2 pi*(3396000*39.37001*72)^2 pi*(71492000*39.37001*72)^2 pi*(60268000*39.37001*72)^2 pi*(695500000*39.37001*72)^2];
for i = 1:N
X = Ra(:,(3*i-2));
Y = Ra(:,(3*i-1));
Z = Ra(:,3*i);
scatter3(Ra(1,(3*i-2)), Ra(1,(3*i-1)), Ra(1,(3*i)), color(i))
plot3(X, Y, Z, color(i))
scatter3(Ra(i_t,(3*i-2)), Ra(i_t,(3*i-1)), Ra(i_t,(3*i)), color(i))
%M(i) = getframe;
axis square
end
view(3)
grid on
xlabel('x')
ylabel('y')
zlabel('z')
xlim([-1.5e9 1.5e9]);
ylim([-1.5e9 1.5e9]);
zlim([-1e8 1e8]);
hold off
|
|
