Path: news.mathworks.com!not-for-mail
From: "Iris " <my.velvet.faerie.tale@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Creating a movie for n-body simulations
Date: Tue, 27 Nov 2012 20:38:12 +0000 (UTC)
Organization: FCUP
Lines: 168
Message-ID: <k938bk$pa9$1@newscl01ah.mathworks.com>
Reply-To: "Iris " <my.velvet.faerie.tale@gmail.com>
NNTP-Posting-Host: www-00-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1354048692 25929 172.30.248.45 (27 Nov 2012 20:38:12 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 27 Nov 2012 20:38:12 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 3755865
Xref: news.mathworks.com comp.soft-sys.matlab:783647

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