Matlab动画模拟太阳系行星运动.docx
《Matlab动画模拟太阳系行星运动.docx》由会员分享，可在线阅读，更多相关《Matlab动画模拟太阳系行星运动.docx（13页珍藏版）》请在冰点文库上搜索。
Matlab动画模拟太阳系行星运动
Matlab动画模拟太阳系行星运动
figure（'name','星系演示'）;%设置标题名字
pausetime=.02;%设置暂停时间
set（gca,'xlim',[-5050],'ylim',[-5030],'zlim',[-5050]）;
set（gcf,'doublebuffer','on'）%消除抖动
xlabel（'x轴'）,ylabel（'y轴'）,zlabel（'z轴'）;
axisequal;
gridon;
view（[352]）;
holdon
a=[203050608010090];b=[81218264555709030];
omga=[411];r=[3];%长轴，短轴，角速度，球体半径
c=sqrt（a.^2-b.^2）;h=pi/18;h1=pi/10;f=pi/9;g=pi/8;
aby=[hh0;h1h0;hh0;hh0;hh0;hh0;hg0;hhh;g0g];%每个轨道平面倾斜角度，偏移设置
%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',[102],'style','infinite'）
lightingphong
materialshiny
t=0:
*pi:
50*pi;
t';
num=length（a）;
forn=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（:
:
）=[100;0cos（aby（n,1））-sin（aby（n,1））;0sin（aby（n,1））cos（aby（n,1））]*[cos（aby（n,2））0sin（aby（n,2））;010;-sin（aby（n,2））0cos（aby（n,2））]*[cos（aby（n,3））-sin（aby（n,3））0;sin（aby（n,3））cos（aby（n,3））0;001];%每个轨道平面倾斜计算
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））;shadinginterp%画出每个行星
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）,'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'）;%设置所有行星的颜色
form=1:
5000%旋转计算
forn=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',[-5050],'ylim',[-5030],'zlim',[-5050]）;
set（gcf,'doublebuffer','on'）%消除抖动
xlabel（'x轴'）,ylabel（'y轴'）,zlabel（'z轴'）;
axisequal;
gridon;
view（[352]）;
holdon
a=[20305060801009045];%长轴
b=[8121826455570903045];%短轴前八个为对应行星，第九个为彗星，后面为卫星
omga=[41436];%角速度
r=[3];%球体半径
c=sqrt（a.^2-b.^2）;h=pi/18;h1=pi/10;f=pi/9;g=pi/8;g1=pi/6;
aby=[hh0;h1h0;hh0;hh0;hh0;hh0;hg0;hhh;g0g;000;g1h0;0f0;000;0g10];%每个轨道平面偏移设置
runu=35:
:
40;theta=（0:
*pi:
2*pi）';runa=:
:
;
xx=cos（theta）*runu+20;yy=*sin（theta）*runu;zz=**yy;
plot3（xx,yy,zz,':
k'）;%小行带设置
hx=cos（theta）*runa;hy=sin（theta）*runa;hz=**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',[102],'style','infinite'）
lightingphong
materialshiny
t=0:
*pi:
50*pi;
t';
num=length（a）;
forn=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（:
:
）=[100;0cos（aby（n,1））-sin（aby（n,1））;0sin（aby（n,1））cos（aby（n,1））]*[cos（aby（n,2））0sin（aby（n,2））;010;-sin（aby（n,2））0cos（aby（n,2））]*[cos（aby（n,3））-sin（aby（n,3））0;sin（aby（n,3））cos（aby（n,3））0;001];
xyz（:
:
）=[x（:
n）y（:
n）z（:
n）]*xuanz（:
:
）;
x（:
n）=xyz（:
1）;
y（:
n）=xyz（:
2）;
z（:
n）=xyz（:
3）;
ifn<=9
p（n）=surf（r（n）*X+x（1,n）,r（n）*Y+y（1,n）,r（n）*Z+z（1,n））;shadinginterp
plot3（x（:
n）,y（:
n）,z（:
n）,'-k'）;%画出所有轨迹线
elseifn<=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））;shadinginterp
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））;shadinginterp
pmuw（n-9）=plot3（x（:
n）+x（1,3）,y（:
n）+y（1,3）,z（:
n）+z（1,3）,'-k'）;%月球轨道设置
end
end
end
forn=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）,'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'）;
forn=10:
13
set（p（n）,'facecolor','b'）;
end
set（p（14）,'facecolor','k'）;
form=1:
5000
forn=1:
num
ifn<=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））;%所有的即时位置
elseifn<=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
forn=1:
length（runa）
set（ph（n）,'xdata',hx（:
n）+x（m,6）,'ydata',hy（:
n）+y（m,6）,'zdata',hz（:
n）+z（m,6））;%光环即时位置
end
pause（pausetime）;%暂停一会
drawnow
end
升级会员 