精通MATLAB的请帮忙

jonwang 发表于 2004-12-30 21:42:41

行星和太阳的质量，行星到太阳的平均距离，行星的平均速度
采用牛顿万有引力定律和加速度定律模拟太阳系的运动轨迹。
运行过程中观察4个行星的运动。修改该程序：
 % 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,'o',Col);
axis([-RANGE RANGE -RANGE RANGE]);
hold off;

drawnow;
end
 end
A添加一个彗星（需要增加Nsphere 数组等）.
假设这个和彗星以一个初速度从Jupiter附近开始运动。
B在地球附近添加月球
C能够演示出彗星撞地球或其他行星；找到月球对地球应力最大的时刻
请问怎么修改?

页: [1]

数模论坛's Archiver

精通MATLAB的请帮忙