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

1 view (last 30 days)

Show older comments

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).

##### 1 Comment

Walter Roberson
on 8 May 2022

### Answers (1)

VBBV
on 8 May 2022

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.

##### 2 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!