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

1 view (last 30 days)
Kit Leeming
Kit Leeming on 8 May 2022
Commented: Walter Roberson on 8 May 2022
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
Walter Roberson on 8 May 2022
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.

Sign in to comment.

Answers (1)

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

Community Treasure Hunt

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

Start Hunting!

Translated by