Error in loop in Matlab, calculate position

Question

GOAL:

In this project, we calculate the gravitational forces exerted by each body in a solar system on the other bodies in the system. Based upon those forces and initial position/velocities of the bodies, one can predict their motion using Newton’s second law. We assume that the following information is given:

The masses of all bodies involved in a solar system.
The positions and velocities of all planets at a given time.

Then we can calculate the gravitational forces using Newton’s law of universal gravitation. Each body is acted upon by gravitational forces from all the other bodies, and therefore the total force on the body is the sum of the forces induced by all the other bodies. Once the force on a body is known, the acceleration of each body can be determined using Newton’s second law of motion. Based upon this acceleration at any time t, we can calculate the new velocity and position of the bodies at time t + ∆t. We start with known positions and velocities at time t = 0, then calculate those at ∆t, then at 2∆t and so on.

ASSIGNMENT:

Write a MATLAB function that takes as input a struct array (where each element represents one body in your solar system), a time step ∆t and a final time T. Each struct in the struct array should have the following ﬁelds:

name: A string that holds the name of the planet.
x: The x-coordinate of the initial position of the planet.
y: The y-coordinate of the initial position of the planet.
z: The z-coordinate of the initial position of the planet.
vx: The x-component of the initial velocity of the planet.
vy: The y-component of the initial velocity of the planet.
vz: The z-component of the initial velocity of the planet.

Then your function should use the model of a gravitational system described above to compute the positions and velocities of all planets at all times. Your function should return these values in an appropriate way and also produce a plot of the positions of all planets in a 3D plot. The plot should at least contain a title and a legend. The legend should use the ﬁeld name in order to name each planet. Further, each planet should be shown in a diﬀerent color where the color is determined randomly. A sample plot is given in Figure 1.

MY ERROR

There is an error in my last loop with time in it, has n=1:T/t; regarding matrix dimensions. I'm trying to calculate new forces using the updated positions from each iteration, however I am getting and error, and my plot shows a bunch of straight lines, and not a solar system orbit.**

function [x, y, z, vx, vy, vz] = solarsystemsim(F,t,T)

m = size(F,2);
G = 6.67384e-11;
j=1;
x = 1:10;
% r = zeros(m-1,1);
% gF = zeros(m-1,1);
for(i=1:m)
    for(j=1:m)
        r(i,j) = distform2(F(i).x,F(j).x,F(i).y,F(j).y,F(i).z,F(j).z);
        gF(i,j) = (G*(F(i).mass)*(F(j).mass))/((r(i,j))^2);
        gFx(i,j) = -((gF(i,j))*(F(i).x-F(j).x))/(r(i,j));
        gFy(i,j) = -((gF(i,j))*(F(i).y-F(j).y))/(r(i,j));
        gFz(i,j) = -((gF(i,j))*(F(i).z-F(j).z))/(r(i,j));

        if(i==j)
            gF(i,j) = 0; 
            gFx(i,j) = 0;
            gFy(i,j) = 0;
            gFz(i,j) = 0;
        end
    end
end

for(i=1:m)
    gFxT(i) = sum(gFx(i,:)); 
    gFyT(i) = sum(gFy(i,:)); 
    gFzT(i) = sum(gFz(i,:)); 
end


tn = 0;
n = 1;

x = zeros(T/t,m);
y = zeros(T/t,m);
z = zeros(T/t,m);
vx = zeros(T/t,m);
vy = zeros(T/t,m);
vz = zeros(T/t,m);

for(i=1:m)
    vx(1,i) = F(i).vx;
    vy(1,i) = F(i).vy;
    vz(1,i) = F(i).vz;
    x(1,i) = F(i).x;
    y(1,i) = F(i).y;
    z(1,i) = F(i).z;
end

for(n=1:T/t)
    for(i=1:m)
        for(j=1:m)
            r(i+1,j+1) = distform2(x(i,j),x(i,j),y(i,j),y(i,j),z(i,j),z(i,j));
            gF(i+1,j+1) = (G*(F(i).mass)*(F(j).mass))/((r(i,j))^2);
            gFx(i+1,j+1) = -((gF(i,j))*(x(i,j)-x(i,j+1)))/(r(i,j));
            gFy(i+1,j+1) = -((gF(i,j))*(y(i,j)-y(i,j+1)))/(r(i,j));
            gFz(i+1,j+1) = -((gF(i,j))*(z(i,j)-z(i,j+1)))/(r(i,j));

            if(i==j)
                gF(i,j) = 0; 
                gFx(i,j) = 0;
                gFy(i,j) = 0;
                gFz(i,j) = 0;
            end
        end
    end

gFxT(i) = sum(gFx(i,:)); 
gFyT(i) = sum(gFy(i,:)); 
gFzT(i) = sum(gFz(i,:)); 

    for(i=1:T/t)
        for(j=1:m)
            vx(i,j) = vx(i,j) + t*(gFxT(j))/(F(j).mass);
            vy(i,j) = vy(i,j) + t*(gFyT(j))/(F(j).mass);
            vz(i,j) = vz(i,j) + t*(gFzT(j))/(F(j).mass);
            x(i,j) = x(i,j) + t*(vx(j));
            y(i,j) = y(i,j) + t*(vy(j));
            z(i,j) = z(i,j) + t*(vz(j));
        end
    end 

plot3(x,y,z);
end

Dave Durbin Dave Durbin · Accepted Answer · 2015-05-06T07:44:51

F is not a matrix and so you can't index it using parentheses. It's a cell array and must be indexed using braces {}

EDITED to correct struct to cell as per lmillefiori's comment

Error in loop in Matlab, calculate position

1 Answers