Matlab动画模拟太阳系行星运动教学教材

Matlab动画模拟太阳系行星运动
figure('name','星系演示');%设置标题名字
pausetime=.02;%设置暂停时间
set(gca,'xlim',[-50 50],'ylim',[-50 30],'zlim',[-50 50]);
set(gcf,'doublebuffer','on') %消除抖动
xlabel('x轴'),ylabel('y轴'),zlabel('z轴');
axis equal;
grid on;
view([3 5 2]);
hold on
a=[8.5 12.5 20 30 50 60 80 100 90];b=[8 12 18 26 45 55 70 90 30]; omga=[4 1.25 1 0.5 0.1 0.05 0.25 0.125 1];r=[0.35 0.8 0.8 0.5 3 2.5 1.5 1.5 0.35];%长轴，短轴，角速度，球体半径
c=sqrt(a.^2-b.^2);h=pi/18;h1=pi/10;f=pi/9;g=pi/8;
aby=[h h 0;h1 h 0;h h 0;h h 0;h h 0;h h 0;h g 0;h h h;g 0 g];%每个轨道平面倾斜角度，偏移设置
%colo={'y','m','b','m','r','c','b','b'};
[X,Y,Z]=sphere(40);
surf(5*X,5*Y,5*Z);colormap(autumn) %设置太阳
light ('position',[1 0 2],'style','infinite')
lighting phong
material shiny
t=0:0.01*pi:50*pi;
t';
num=length(a);
for n=1:num
x(:,n)=a(n)*cos(omga(n)*t)+c(n);
y(:,n)=b(n)*sin(omga(n)*t);
z(:,n)=0*t; %计算未经轨道平面角度倾斜的轨道位置
xuanz(:,:)=[1 0 0;0 cos(aby(n,1)) -sin(aby(n,1));0 sin(aby(n,1)) cos(aby(n,1))]*[cos(aby(n,2)) 0 sin(aby(n,2));0 1 0;-sin(aby(n,2)) 0 cos(aby(n,2))]*[cos(aby(n,3)) -sin(aby(n,3)) 0;sin(aby(n,3)) cos(aby(n,3)) 0;0 0 1]; %每个轨道平面倾斜计算
xyz(:,:)=[x(:,n) y(:,n) z(:,n)]*xuanz(:,:);
x(:,n)=xyz(:,1);
y(:,n)=xyz(:,2);
z(:,n)=xyz(:,3); %计算轨道平面倾斜后的轨道位置
p(n)=surf(r(n)*X+x(1,n),r(n)*Y+y(1,n),r(n)*Z+z(1,n));shading
interp %画出每个行星
plot3(x(:,n),y(:,n),z(:,n),'-k');%画出所有轨迹线
end
set(p(1),'facecolor','y');
set(p(2),'facecolor','m');set(p(3),'facecolor','b');set(p(4),'fac ecolor','m');
set(p(5),'facecolor','r');set(p(6),'facecolor','c');set(p(7),'fac ecolor','b');set(p(8),'facecolor','b');set(p(9),'facecolor','r'); %设置所有行星的颜色
for m=1:5000 %旋转计算
for n=1:length(a)
set(p(n),'xdata',r(n)*X+x(m,n),'ydata',r(n)*Y+y(m,n),'zdata',r(n) *Z+z(m,n)); %所有行星的即时位置设置
end
pause(pausetime); %暂停一会
drawnow
end
下面是更加复杂的动画模拟
figure('name','星系演示');%设置标题名字
pausetime=.01;%设置暂停时间
set(gca,'xlim',[-50 50],'ylim',[-50 30],'zlim',[-50 50]);
set(gcf,'doublebuffer','on') %消除抖动
xlabel('x 轴'),ylabel('y 轴'),zlabel('z 轴');
axis equal ;
grid on ;
view([3 5 2]);
hold on
a=[8.5 12.5 20 30 50 60 80 100 90 4 4.5 4.9 5 1.5];%长轴
b=[8 12 18 26 45 55 70 90 30 4 4.5 4.9 5 1.5];%短轴前八个为对应行星，第九个为彗星，后面为卫星
omga=[4 1.25 1 0.5 0.1 0.05 0.25 0.125 0.4 4 3.9 3.5 3 6];%角速度
r=[0.35 0.8 0.8 0.5 3 2.5 1.5 1.5 0.5 0.35 0.36 0.5 0.4 0.35];%球体半径
c=sqrt(a.^2-b.^2);h=pi/18;h1=pi/10;f=pi/9;g=pi/8;g1=pi/6;
aby=[h h 0;h1 h 0;h h 0;h h 0;h h 0;h h 0;h g 0;h h h;g 0 g;0 0 0;g1 h 0;0 f 0;0 0 0;0 g1 0];%每个轨道平面偏移设置
runu=35:0.5:40;theta=(0:0.05*pi:2*pi)';runa=2.8:0.4:5.6;
xx=cos(theta)*runu+20;yy=0.9*sin(theta)*runu;zz=-0.17*xx-0.17*yy; plot3(xx,yy,zz,':k'); %小行带设置
hx=cos(theta)*runa;hy=sin(theta)*runa;hz=-0.1*hx-0.2*hy;
%colo={'y','m','b','m','r','c','b','b'};
[X,Y,Z]=sphere(40);
surf(5*X,5*Y,5*Z);colormap(autumn) %设置
light ('position',[1 0 2],'style','infinite')
lighting phong
material shiny
t=0:0.01*pi:50*pi;
t';
num=length(a);
for n=1:num
x(:,n)=a(n)*cos(omga(n)*t)+c(n);
y(:,n)=b(n)*sin(omga(n)*t);
z(:,n)=0*t;
xuanz(:,:)=[1 0 0;0 cos(aby(n,1)) -sin(aby(n,1));0 sin(aby(n,1)) cos(aby(n,1))]*[cos(aby(n,2)) 0 sin(aby(n,2));0 1 0;-sin(aby(n,2)) 0 cos(aby(n,2))]*[cos(aby(n,3)) -sin(aby(n,3)) 0;sin(aby(n,3)) cos(aby(n,3)) 0;0 0 1];
xyz(:,:)=[x(:,n) y(:,n) z(:,n)]*xuanz(:,:);
x(:,n)=xyz(:,1);
y(:,n)=xyz(:,2);
z(:,n)=xyz(:,3);
if n<=9
p(n)=surf(r(n)*X+x(1,n),r(n)*Y+y(1,n),r(n)*Z+z(1,n));shading
interp
plot3(x(:,n),y(:,n),z(:,n),'-k');%画出所有轨迹线
else if n<=13
p(n)=surf(r(n)*X+x(1,n)+x(1,5),r(n)*Y+y(1,n)+y(1,5),r(n)*Z+z(1,n)
+z(1,5));shading interp
pmuw(n-9)=plot3(x(:,n)+x(1,5),y(:,n)+y(1,5),z(:,n)+z(1,5),'-k');%木卫1,2,3,4轨道初位置
else
p(n)=surf(r(n)*X+x(1,n)+x(1,3),r(n)*Y+y(1,n)+y(1,3),r(n)*Z+z(1,n) +z(1,3));shading interp
pmuw(n-9)=plot3(x(:,n)+x(1,3),y(:,n)+y(1,3),z(:,n)+z(1,3),'-k');%月球轨道设置
end
end
end
for n=1:length(runa)
ph(n)=plot3(hx(:,n)+x(1,6),hy(:,n)+y(1,6),hz(:,n)+z(1,6),'-c');
end
set(p(1),'facecolor','y');
set(p(2),'facecolor','m');set(p(3),'facecolor','b');set(p(4),'fac ecolor','m');
set(p(5),'facecolor','r');set(p(6),'facecolor','c');set(p(7),'fac ecolor','b');set(p(8),'facecolor','b');set(p(9),'facecolor','r'); for n=10:13
set(p(n),'facecolor','b');
end
set(p(14),'facecolor','k');
for m=1:5000
for n=1:num
if n<=9
set(p(n),'xdata',r(n)*X+x(m,n),'ydata',r(n)*Y+y(m,n),'zdata',r(n) *Z+z(m,n));%所有的即时位置
else if n<=13
set(p(n),'xdata',r(n)*X+x(m,n)+x(m,5),'ydata',r(n)*Y+y(m,n)+y(m,5 ),'zdata',r(n)*Z+z(m,n)+z(m,5));
set(pmuw(n-9),'xdata',x(:,n)+x(m,5),'ydata',y(:,n)+y(m,5),'zdata' ,z(:,n)+z(m,5));
else
set(p(n),'xdata',r(n)*X+x(m,n)+x(m,3),'ydata',r(n)*Y+y(m,n)+y(m,3 ),'zdata',r(n)*Z+z(m,n)+z(m,3));
set(pmuw(n-9),'xdata',x(:,n)+x(m,3),'ydata',y(:,n)+y(m,3),'zdata' ,z(:,n)+z(m,3));
end
end
end
for n=1:length(runa)
set(ph(n),'xdata',hx(:,n)+x(m,6),'ydata',hy(:,n)+y(m,6),'zdata',h z(:,n)+z(m,6));%光环即时位置
end
pause(pausetime); %暂停一会
drawnow
end。