Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Creating a movie for n-body simulations

Subject: Creating a movie for n-body simulations

From: Iris

Date: 27 Nov, 2012 20:38:12

Message: 1 of 4

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
    

Subject: Creating a movie for n-body simulations

From: Nasser M. Abbasi

Date: 28 Nov, 2012 01:24:09

Message: 2 of 4

On 11/27/2012 2:38 PM, Iris wrote:
> Hello guys! I've done a n-body simulation program in matlab,

You should really fix all the problems in the code first.
A simple use of code anlyzer shows these


8: The value assigned to variable 'R' might be unused.

46: The variable 'Rab' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

47: The variable 'Vab' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

53: Assignment to 'dt' might be unnecessary.

57: The variable 'cm' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

69: The variable 'cmv' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

92: Use of brackets [] is unnecessary. Use parentheses to group, if needed.

103: The variable 'b' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

108: The variable 'pbt' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

109: The variable 'accelerations' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

129: The variable 'Rab' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

130: The variable 'Vab' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

131: The variable 'Aab' appears to change size on every loop iteration (within a script). Consider preallocating for speed.

137: The variable 'p' appears to change size on every loop iteration (within a script). Consider preallocating for speed

Subject: Creating a movie for n-body simulations

From: Iris

Date: 28 Nov, 2012 11:27:06

Message: 3 of 4

Hello Nasser!
The code is working fine, dont worry about that.
Any idea of how can I do a movie like that?
Thank you

Subject: Creating a movie for n-body simulations

From: Nasser M. Abbasi

Date: 28 Nov, 2012 11:47:53

Message: 4 of 4

On 11/28/2012 5:27 AM, Iris wrote:
> Hello Nasser!
> The code is working fine, dont worry about that.
> Any idea of how can I do a movie like that?
> Thank you
>


I am not an expert on this. BUt I'll show what I did once
to make animated gif file.

--------------------------------
%initialize section, first time only
    
    f = getframe(gcf);
    [im,map] = rgb2ind(f.cdata,256,'nodither');
    im(1,1,1,20) = 0;
    frame_number = 0;

    for i=1:100 %simulation loop, whatever your loop is
        %... after making all you update to figure

        f = getframe(gcf);
        frame_number = frame_number+1;
        im(:,:,1,frame_number)=rgb2ind(f.cdata,map,'nodither');
    end

    %done simulation
    imwrite(im,map,'my_animation.gif','DelayTime',0,'LoopCount',inf)

--------------------------------------

The above for animated gif file. For other types of 'movies'
I am not sure now since I only done it for animated gif.

I can't guarantee this will work, I just copied pieces out
of old code I have.

--Nasser

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us