Creating a movie for n-body simulations

1 次查看（过去 30 天）
Iris
Iris 2012-11-27
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
end

请先登录，再回答此问题。

回答（2 个）

Babak
Babak 2012-11-27
use "getframe" and "addframe" to "avifile" object. Look them up in help to learn how to make a movie. if you have the center of the spheres and their radii, you can plot them in 2D or 3D and then getframe and addframe to your movie file. You can also draw the path of the motion of the planets and superimpose it to the figure.
  1 个评论
Babak
Babak 2012-11-28
编辑：Babak 2012-11-28
After getting all the time step's solutions, , you create a for loop for plotting. In the for loop you plot the scatter and then
hold on
and then plot the line (path), then you will do a
getframe
then you should
hold off
This way, in the next iteration of your for loop, you will plot the new scatter and plot the new path then getframe,...
Note that, if you are going to do this, the path (line) you want to plot, in every iteration of the for loop should include all the previous path as well as the path of the current time step. To do this, you need to create a matrix that saves the previous path when you are deriving the solution. Like this:
for j=1:N
% compute line
path = [path , line];
end
This way you can save the lines and when plotting you can plot only
path(1:j)

请先登录，再进行评论。

Iris
Iris 2012-11-28
Hello Babak, thank you for your answer. The problem is that i have 2 plots (the line and the scatter) and at each time i need to erase the previous scatter but remain with the line... Thank u *

类别

在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by