[roshan]

万有引力公式：F=Gm1m2/d12^2
现在有若干个星球在太空中由于引力相互作用，但是我的代码显示出来的若干个星球只能沿一个方向作直线运动，而不是相互环绕运行。如何才能让他们正确运行，请教各位大侠指教，等待指点迷津。万分感谢。
———————————————————题目大意：————————————————————————

the motion of bodies in 2-d space is to be investigated. the bodies
together constitute a small universe, and as such they both create and
are influenced by a gravitional field

Newton/Kepler demostrated that the total force Fi in the general case
of 50 masses is equal to
Fi=Gm1m2/d12^2+Gm1m3/d13^2+.......+Gm1m50/d150^2

create initial positions, velocities and displacement for 50 masse , display the motion graphically

when gets together,the force is big, the remidies include: place a
limit on max. force between masses. decrese the time step. set a
collision distance, when 2 masses are closer than this to each other,
they are assumed to collide, and become 1 body, with mass and momentum
equal to the sum of the pre-collision values.

———————————————————以下代码：————————————————————————

clear all;
close all;

G=6.67*10^-11;
% The universal gravitational constant

m = [1.989e30,3.5844e23,4.89868e24,5.974e24,6.5714e23,1.89854e27,5.68725e26,8.72204e25,1.02753e26];
% An array of masses

n = length(m); Number_of_planets = n
% The number of masses

x_p = [0,58340100000,1.07705e11,1.4959e11,2.27377e11,7.77868e11,1.42709e12,2.87512e12,4.49668e12];
% An array of x positions
y_p= [0,0,0,0,0,0,0,0,0];
% An array of y positions

x_v = [0,0,0,0,0,0,0,0,0];
% An array of x velocities
y_v = [0,47856.46,34961.72,29780,23883.56,12924.52,9618.94,6789.84,5419.96];
% An array of y velocities

T= 365*24*60*60;%(2*pi*G*m(1))/(x_v(1)^3);
dt=360*60*60*10;

colordef black;
ph = plot(x_p,y_p,'.','MarkerSize',30);
xlabel('distance(km)');
ylabel('distance(km)');
title('Planetary motion');

for a=0:dt:1000*T
%loop for all masses with respect to time

for k=1:n
%loop for individual masses calculating neccessary quantities

dx = x_p - x_p(k);
% difference in x positions
dy = y_p - y_p(k);
% difference in y positions
mag = (dx.^2 + dy.^2).^0.5;
% magnitude of the distance between the 2 masses

f = ((G*m(k)*m)./(mag.^2));
% total force

fx_old =(f.*dx./mag);
fx_old(k)= 0;
fx = sum (fx_old);
%Summing the total force in x direction for each planet

fy_old =(f.*dy./mag);
fy_old(k)= 0;
fy = sum (fy_old);
%Summing the total force in y direction for each planet


a_x = fx/m(k);
% calculating acceleration change in x direction
a_y = fy/m(k);
% calculating acceleration change in y direction

x_v(k) = x_v(k) + dt*a_x;
% calculating velocity change in x direction
y_v(k) = y_v(k) + dt*a_y;
% calculating velocity change in y direction

x_p_new(k)= x_p(k)+dt*x_v(k);
% calculating new x positions
y_p_new(k)= y_p(k)+dt*y_v(k);
% calculating new y positions
end

x_p=x_p_new;
y_p=y_p_new;
% overwriting old positions with new positions

pause(0.1)
plot(x_p,y_p,'.','MarkerSize',30)
axis([-(10^15) 10^15 -(10^15) 10^15])
xlabel('distance(km)');
ylabel('distance(km)');
title('Planetary motion');
drawnow
% Plotting graph
end

———————————————————华丽的分割线————————————————————————