Having trouble writing a multiple gravitational physics simulation code, as every planet effects one another.

Here is my code so far:

function [p, v] = solarsystem(p, v, mass, stop_time, hide_animation) 
%p and v are matrices that describe the position and velocity for a set of
%planets
if nargin < 5
    hide_animation = false;
end
G=6.67*10^(-11); %Gravitational Constant (Nm^2/kg^2)
ts=250; %timestep (s)
jump = 0;

o=size(p,1);  %extraction number of planets
n=1;
F=zeros(o*(o-1),2);
r=zeros(o*(o-1),2);
a=zeros(o*(o-1),2);

if hide_animation==0
    figure()
    s=plot(p(1,1),p(1,2),'yo');
    hold on
    e=plot(p(2,1),p(2,2),'bo');
    axis equal
    xlim([-1.6*10^11,1.6*10^11])
    ylim([-1.6*10^11,1.6*10^11])
end 

while (n*ts)<stop_time
   c=1;
    for i = 1;o;        %These nested for loops are to try and account for all combinations of each planet
        for j = 1;o;
            if i~=j
                r(c,:)=p(j,:)-p(i,:); %calulate the distance between each of the planets
                F(c,:)=((G*mass(i)*mass(j))/(norm(r(c))^3)).*r(c); %Calulate the forces acting on each of the planets
                c=c+1;
            end
        end
        for q = 1:(o-1):(o*(o-1))
            a(i,:)=(sum(F(q:(q+o-1),:)))./mass(i); %Calculate acceleration of each planet
            v(i,:)=v(i,:)+a(i,:)*ts;    %calulate new velocity using acceleration found in previous line
            p(i,:) = p(i,:)+v(i,:)*ts;  %calulate new position
        end
    end

    
        if hide_animation == 0
            if jump == 40
                set(s, 'Xdata', p(1,1), 'Ydata', p(1,2));       %Update each of the plot points for the animation
                set(e, 'Xdata', p(2,1), 'Ydata', p(2,2));
                drawnow limitrate;
                jump=0;
            else
                jump=jump+1;
            end
        end
   
 n=n+1;
end

Right now I'm getting errors saying Index in position 1 exceeds array bounds. Index must not exceed 2.

Error in solarsystem (line 83)

a(i,:)=(sum(F(q:(q+o-1),:)))./mass(i);

Any help would be appreciated :).

(I also haven't entirely finished this yet so I can still only plot two planets interacting at the moment).

[Walter Roberson]

When I read your question I suddenly realized that for acurrate simulations, you should not be calculating mutual gravitational influence based upon the "current" location of the objects: that you should be calculating based upon the location where it was a time in the past corresponding to the time light would take to travel between the objects. Gravity has a speed, which is the speed of light.

Answer 1

[VBBV]

a(q,:)=(sum(F(q:(q+o-1),:)))./mass(q); %Calculate acceleration of each planet
          v(q,:)=v(q,:)+a(q,:)*ts;    %calulate new velocity using acceleration found in previous line
          p(q,:) = p(q,:)+v(q,:)*ts;  %calulate

Use q as the index for these matrices.

[VBBV]

for i = 1:o       %These nested for loops are to try and account for all combinations of each planet
      for j = 1:o

[Walter Roberson]

Notice that : needs to be used in the for loops.