Simulation of the motion of planets under the effect of the gravitational force

[OmegaForce]

Hello.

I tried to make a program which simulated the movement of planets under the effect of the gravitational force. Here is my program:

n=5;
mass=1;
speed=0;
scale=3000;
distance=1;
time=200;
p=0.001;
G=1;

x=distance*(2*rand(1,n)-1);
y=distance*(2*rand(1,n)-1);
m=mass*rand(1,n);
vx=speed*(2*rand(1,n)-1);
vy=speed*(2*rand(1,n)-1);
fx=zeros(1,n);
fy=zeros(1,n);

for t=1:time
    
for a=1:n
    
    for b=1:n
        
        if b~=a
           
           fx(a)=0;
           fy(a)=0; 
            
           fx(a)=(G*m(b)/((y(b)-y(a))^2+(x(b)-x(a))^2))*((x(b)-x(a))/sqrt((x(b)-x(a))^2+(y(b)-y(a))^2));
           fy(a)=(G*m(b)/((y(b)-y(a))^2+(x(b)-x(a))^2))*((y(b)-y(a))/sqrt((x(b)-x(a))^2+(y(b)-y(a))^2));
           
           vx(a)=vx(a)+fx(a);
           vy(a)=vy(a)+fy(a);
        
           x(a)=x(a)+vx(a);
           y(a)=y(a)+vy(a);
           
        end
    
    end
    
    plot(x(a),y(a),'o')
    xlim([-scale,scale])
    ylim([-scale,scale])
    hold on
    
end

pause(p)
hold off

end

The beginning is just some parameters which we can adjust. The parameter 'n' is the number of bodies in the simulation. 'p' is the time of pause/break in seconds between two graphs. 'G' is the value of the gravitational constant. 'fx' is the force along the axis x, 'fy' is the force along the axis y. 'vx' is the speed along the axis x, 'vy' is the speed along the axis y. 'x' is the position along the axis x, 'y' is the position along the axis y. The initial position, the initial speed and the initial mass of each body is taken/generated randomly.

Then I intend to put values near the reality: G = 6,67 . 10^(-11) which is the true value of the gravitational constant, distance = 10^11 which is a typical order of magnitude in the Solar System (in meters), mass = 10^25 which is an order of magnitude for the mass of the planets (in kg), etc.

But the problem is that my program doesn't work at all: instead of having bodies which attract each other, they repulse each other... My program simulates anti-gravity!...

I really don't understand why. I have read and re-read my program for hours and I haven't found where are the errors... Note that I am only a beginner in programming, so maybe there's just something that escapes me...

Thanks in advance for your answers.

 

[SomptingGuy]

You may need an m(a) in the force calc. And a minus sign.


Steve