Creating a movie for n-body simulations

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

Answer 1

[Babak]

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.

[Babak]

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)

Answer 2

[Iris]

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 *