太阳系模拟求助

yemulg

下面是一个已经写好了的源程序,里面有行星和太阳的质量，行星到太阳的平均距离，行星的平均速度。采用牛顿万有引力定律和加速度定律模拟太阳系的运动轨迹。运行过程中观察4个行星的运动。
  现在想添加一个彗星（需要增加Nsphere 数组等）.假设这个和彗星以一个初速度从Jupiter
附近开始运动，能够演示出彗星撞地球或其他行星。
    想了很久不知道如何实现,感觉要通过各个天体之间的相互作用最后再相撞,对彗星的质量,初始位置和处速度都必须有一定的要求才行,但不知道如何确定这些量,还望指教.
% total number of big objects in solar system
Nsphere  = 10;

% color of planets
Col = [1:Nsphere];

% Gravitational constant in kg^-2*m^3*s^-2
Grav     = 6.67259*(10^(-11));

% mass of planets in kg
Mass     = zeros(Nsphere,1);
Mass(1)  = 0.33*(10^24);       % Mercury
Mass(2)  = 4.87*(10^24);       % Venus   
Mass(3)  = 5.97*(10^24);       % Earth   
Mass(4)  = 0.642*(10^24);      % Mars   
Mass(5)  = 1899*(10^24);       % Jupiter
Mass(6)  = 568*(10^24);        % Saturn  
Mass(7)  = 86.8*(10^24);       % Uranus  
Mass(8)  = 102*(10^24);        % Neptune
Mass(9)  = 0.0125*(10^24);     % Pluto   
Mass(10) = 1047.3486*Mass(5);  % The sun

% initial position in m
X        = zeros(Nsphere,1);
Y        = zeros(Nsphere,1);
X(1)     =   57.9*(10^9);      % Mercury
X(2)     =  108.2*(10^9);         % Venus  
X(3)     =  149.6*(10^9);        % Earth  
X(4)     =  227.9*(10^9);         % Mars   
X(5)     =  778.6*(10^9);        % Jupiter
X(6)     = 1433.5*(10^9);        % Saturn
X(7)     = 2872.5*(10^9);        % Uranus
X(8)     = 4495.1*(10^9);        % Neptune
X(9)     = 5870.0*(10^9);        % Pluto  
X(10)    = 0;                        % The sun

scatter(X,Y,'or');

% initial velocity in m*s^-1 (using mean orbital velocity)
VX        = zeros(Nsphere,1);
VY        = zeros(Nsphere,1);

VY(1)   = 47.89*(10^3);  % Mercury
VY(2)   = 35.03*(10^3);  % Venus   
VY(3)   = 29.79*(10^3);  % Earth   
VY(4)   = 24.13*(10^3);  % Mars   
VY(5)   = 13.06*(10^3);  % Jupiter
VY(6)   =  9.64*(10^3);  % Saturn  
VY(7)   =  6.80*(10^3);  % Uranus  
VY(8)   =  5.40*(10^3);  % Neptune
VY(9)   =  4.70*(10^3);  % Pluto   
VY(10)  =  0;                  % The sun

% acceleration of spheres
AX        = zeros(Nsphere,1);
AY        = zeros(Nsphere,1);

% Use Newton's law of gravitation to move the planets in space
dt = 300;
day = 24*60*60/dt;

for tstep=1:100000000

   for sphere=1:Nsphere
     rest  = setdiff( (1:Nsphere),sphere);
     dists = (X(rest)-X(sphere)).^2+(Y(rest)-Y(sphere)).^2;
     dists = dists.^(-3/2);

     AX(sphere) = Grav*Mass(rest)'*((X(rest)-X(sphere)).*dists);
     AY(sphere) = Grav*Mass(rest)'*((Y(rest)-Y(sphere)).*dists);
   end
   
   X  =  X + dt*VX;
   Y  =  Y + dt*VY;
   
   VX = VX + dt*AX;
   VY = VY + dt*AY;

   if(~mod(tstep,day))
     title(sprintf('day=%d', tstep/day));

     % only plot sun and first four planets
     RANGE = max(max(abs(X(1:4))),max(abs(Y(1:4))));
     hold on;
     scatter(X,Y,Col);
     axis([-RANGE RANGE -RANGE RANGE]);
     hold off;
     drawnow;
     end

end

xjzuo

“三体”问题就已经不能解析求解了, “多体”问题即使计算机模拟也不好做。但如果你利用"万有引力定律",并考虑"力的叠加原理"进行模拟,也许可行.　＂对彗星的质量,初始位置和处速度＂这可能要查一下专业的天体资料，比如哈雷慧星的记录资料．
下面是一个 模拟地、日、月运动 的程序，转自萝卜。
作个参考：

% 模拟太阳系运动
t=linspace(0,2*pi,100);
fill(cos(t),sin(t),'r');
hold on;
plot(4*cos(t),sin(t)*4,'k');
set(gca,'position',[0 0.11 0.775 0.815])
a=0.1;b=0;
xe=4*cos(a)+cos(t)*0.6;
ye=4*sin(a)+sin(t)*0.6;
He=fill(xe,ye,'b');
xm=4*cos(a)+cos(b);
ym=4*sin(a)+sin(b);
set(gcf,'doublebuffer','on');
Hm=plot(xm,ym,'c.','markersize',24);
aa=gca;
axis([-6,6,-6,6]);
axis square;
k=1;da=0.1;db=0.5;
xlabel('Please press "space" key and stop this program!',...
  'fontsize',12,'color','r');
title('simulate solar system')
axes('position',[0.75,0.11,0.25,0.8]);
fill(0.2+cos(t)*0.18,0.75+sin(t)*0.08,'r');
ylim([0,1]);xlim([0,0.9]);
text(0.5,0.75,'Sun');hold on;
fill(0.2+cos(t)*0.11,0.5+sin(t)*0.05,'b');
text(0.5,0.5,'Earth');
plot(0.2,0.3,'c.','markersize',24);
text(0.5,0.3,'Moon');
axis off
axes(aa);
while k;
   s=get(gcf,'currentkey');
   if strcmp(s,'space');
       clc;k=0;
   end
   a=a+da;
   b=b+db;
   xe=4*cos(a)+cos(t)*0.6;
   ye=4*sin(a)+sin(t)*0.6;
   xm=4*cos(a)+cos(b);
   ym=4*sin(a)+sin(b);
   set(He,'xdata',xe,'ydata',ye);
   set(Hm,'xdata',xm,'ydata',ym);
   pause(0.1);
   if a<80;
       plot(xm,ym);
   end
end
figure(gcf);

[ 本帖最后由 xjzuo 于 2006-11-12 13:30 编辑 ]

yemulg

恩,谢谢指导了