太阳系模拟求助
下面是一个已经写好了的源程序,里面有行星和太阳的质量，行星到太阳的平均距离，行星的平均速度。采用牛顿万有引力定律和加速度定律模拟太阳系的运动轨迹。运行过程中观察4个行星的运动。现在想添加一个彗星（需要增加Nsphere 数组等）.假设这个和彗星以一个初速度从Jupiter
附近开始运动，能够演示出彗星撞地球或其他行星。
想了很久不知道如何实现,感觉要通过各个天体之间的相互作用最后再相撞,对彗星的质量,初始位置和处速度都必须有一定的要求才行,但不知道如何确定这些量,还望指教.
% total number of big objects in solar system
Nsphere= 10;
% color of planets
Col = ;
% 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
回复
“三体”问题就已经不能解析求解了, “多体”问题即使计算机模拟也不好做。但如果你利用"万有引力定律",并考虑"力的叠加原理"进行模拟,也许可行. ＂对彗星的质量,初始位置和处速度＂这可能要查一下专业的天体资料，比如哈雷慧星的记录资料．下面是一个 模拟地、日、月运动 的程序，转自萝卜。
作个参考：
% 模拟太阳系运动
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',)
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',);
fill(0.2+cos(t)*0.18,0.75+sin(t)*0.08,'r');
ylim();xlim();
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 编辑 ] 恩,谢谢指导了
页:
[1]