Matlab动画模拟太阳系行星运动.doc

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 ona=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 phongmaterial shiny t=0:0.01*pi:50*pi;t;num=length(a); for n=1:numx(:,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);%画出所有轨迹线endset(p(1),facecolor,y);set(p(2),facecolor,m);set(p(3),facecolor,b);set(p(4),facecolor,m);set(p(5),facecolor,r);set(p(6),facecolor,c);set(p(7),facecolor,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); %所有行星的即时位置设置 endpause(pausetime); %暂停一会drawnowend下面是更加复杂的动画模拟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 ona=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 phongmaterial shiny t=0:0.01*pi:50*pi;t;num=length(a); for n=1:numx(:,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=9p(n)=surf(r(n)*X+x(1,n),r(n)*Y+y(1,n),r(n)*Z+z(1,n);shading interpplot3(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);%月球轨道设置 endend 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); endset(p(1),facecolor,y);set(p(2),facecolor,m);set(p(3),facecolor,b);set(p(4),facecolor,m);set(p(5),facecolor,r);set(p(6),facecolor,c);set(p(7),facecolor,b);set(p(8),facecolor,b);set(p(9),facecolor,r);for n=10:13 set(p(n),facecolor,b);endset(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