近来整理硬盘发现了一个远古旧物(如下所示)，注意看最后修改时间。掐指一算这个程序已经有16岁了……
那时，MATLAB的版本是6.0，现在都是2019了……
那时，MATLAB安装完占磁盘空间只有600M+，现在2019版完整安装占磁盘28G+……
那时，MATLAB最多的运用还是最传统的数值计算，现在，江湖上的传言是：除生孩子外，MATLAB啥都能做……
那时，我还只个大二学生，现在，我教过的大二学生都已经博士毕业在高校教大二学生了……
这是当年我在学习《理论力学》这门课程时提交的研究型大作业，当时全年级就我一人提交了这个作业而获得了免试的资格。回想起来，这应该是我人生中最后的高光时刻了吧o(╥﹏╥)o
讲真，这个程序现在仍然能跑起来，虽然无比简陋，但是作为一个初学入门者的作品，包含了MATLAB的矩阵运算、数值计算、GUI设计、绘图、动画等功能，应该算是还可以了吧。
——————华丽分割线——————
这个程序解决的是这么一个问题：有一如下图所示的卫星，过其质心的对称轴任意指向某一点。现要求通过控制卫星上的喷嘴的喷速和方向，使卫星在最短的时间内，令其对称轴指向某一指定点。试分析卫星在此过程中的运动情况。(是的，你没看错，下面这个东东它就是个卫星(￣３￣)a ￣ω￣=)
对于如上图所示的卫星模型，其对称轴和指定点可以确定一个平面。如果卫星在这个平面内绕质心作定点转动而去指向指定点，显然可以使其转过的角度最小，从而时间最短。但由于卫星喷嘴的方向不确定性，它们与确定平面间可能有夹角，故卫星的定点转动又可分为两种情况：一种情况是通过控制喷嘴的喷速，使卫星在与确定平面垂直的方向上受力平衡，而在确定平面内直接旋转过去，这时各个喷嘴的喷速大小是不一样的；另一种情况是先控制喷嘴的喷出方向，使卫星绕自身的对称轴作定轴转动，从而将与确定平面夹角最小的那个喷嘴转到确定平面内，然后再在确定平面内定点旋转过去，这时由于喷嘴对确定面具有对称性，故沿确定面对称的喷嘴的喷速大小一样。
根据以上的分析作如下的假设：
在建立的空间三维坐标系下，初始时刻卫星的对称轴指向空间点N(0，0，-1)，指定点P的坐标为(a，b，c)；
喷嘴方向与确定面的夹角分别为α1，α2，α3；
对称轴方向与质心和指定点的连线方向夹角为β；
各个喷嘴的喷速分别为Q1，Q2，Q3，对应的力为F1，F2，F3；
假设喷嘴的最大相对喷速为Vm,相应的最大喷量为Qm=KVm,则对应的最大力为Fmax=QmVm=KVm2;
卫星上下喷嘴之间的距离为L，喷嘴杆的长度为r；
卫星绕对称轴作定轴转动的转动惯量为J1，绕质心作定点转动的转动惯量为J2。
接下来是两个模型：(好吧，我知道应该没人愿意看下面这堆数学公式的ㄟ( ▔, ▔ )ㄏ ㄟ( ▔, ▔ )ㄏ ㄟ( ▔, ▔ )ㄏ )
——————华丽分割线——————
总之，按照这样的设计，构建了一个MATLAB程序，包含3个函数和1个主程序：
模型一：animation1函数
function animation1(action) cla resetxa=0;ya=1;za=0;xb=sqrt(3)/2;yb=-1/2;zb=0;xc=-sqrt(3)/2;yc=-1/2;zc=0; %%A,B,C坐标xn=0;yn=0;zn=-1; %%初始卫星所对准的点Na=str2num(get(findobj(gcf,'tag','h_e1'),'string'));b=str2num(get(findobj(gcf,'tag','h_e2'),'string'));c=str2num(get(findobj(gcf,'tag','h_e3'),'string'));xp=a;yp=b;zp=c; %%要求卫星所对准的点Palpha1=acos((xa*xp+ya*yp)/(sqrt(xa^2+ya^2)*sqrt(xp^2+yp^2))); alpha2=acos((xb*xp+yb*yp)/(sqrt(xb^2+yb^2)*sqrt(xp^2+yp^2))); alpha3=acos((xc*xp+yc*yp)/(sqrt(xc^2+yc^2)*sqrt(xp^2+yp^2))); A=[alpha1; alpha2; alpha3];alpha=min(A,[],1); %%计算alpha beta=acos((xn*xp+yn*yp+zn*zp)/sqrt(xn^2+yn^2+zn^2)/sqrt(xp^2+yp^2+zp^2)); %%计算beta v=50; %%喷气对卫星本身的相对速度,单位是m/s Q=0.004*v; %%喷气喷量最大值,单位是Kg/s J1=30000; %%卫星对于过其质心一轴的转动惯量,单位是Kg*m^2 J2=20000; %%卫星对于过其质心另一轴的转动惯量，单位是Kg*m^2 r=1; %%6个喷嘴到卫星转动轴的距离,单位是m l=10; %%上下两组喷嘴平面之间的距离，单位是m a1=6*Q*v*r/J2;Time1=sqrt(alpha/a1); %%计算第一次加速减速总时间 a2=2*Q*v*l/J1;Time2=sqrt(beta/a2); %%计算第二次加速减速总时间 %------------------------------------------------------------------------------------------ H0=plot3([0,0],[0,0],[-4,4],'m'); %%画卫星实体hold on[h1x,h1y,h1z]=ellipsoid(0,0,0,0.5,0.5,0.5,40); H1=[h1x,h1y,h1z];surfl(h1x,h1y,h1z);shading interp;hold on[h2x,h2y,h2z]=ellipsoid(0,0,0,0.75,0.75,0.25,40); H2=[h2x,h2y,h2z];surfl(h2x,h2y,h2z);shading flat;hold onH4=plot3([0,0],[0,1],[0,0],'r--',... [0,sqrt(3)/2],[0,-1/2],[0,0],'r--',[0,-sqrt(3)/2],[0,-1/2],[0,0],'r--'); %%画OA,OB,OC hold on H5=plot3([0,0],[0,1],[0.75,0.75],'b',... [0,sqrt(3)/2],[0,-1/2],[0.75,0.75],'b',[0,-sqrt(3)/2],[0,-1/2],[0.75,0.75],'b'); %%画OA1,OB1,OC1hold onH6=plot3([0,0],[0,1],[-0.75,-0.75],'b',... [0,sqrt(3)/2],[0,-1/2],[-0.75,-0.75],'b',[0,-sqrt(3)/2],[0,-1/2],[-0.75,-0.75],'b'); %%画OA2,OB2,OC2 H=get(gca,'children'); %%获得当前图形对象的句柄hold onH3=plot3([-1.5,1.5],[0,0],[0,0],'k ',... [0,0],[-1.5,1.5],[0,0],'k',... [0,0],[0,0],[-2.5,2.5],'k'); text(1.7,0,0,'x','fontsize',14); text(0,1.7,0,'y','fontsize',14); text(0,0,2.7,'z','fontsize',14); %%画x,y,z轴hold onplot3(a,b,c,'r+'); text(a+0.1,b+0.1,c+0.1,'P','fontsize',14); %%标定P点N=30;t=Time1/N;for k=1:N alpha(k)=a1*(k^2*t^2-(k-1)^2*t^2)/2; rotate(H,[0,0,-1],180/pi*alpha(k),[0,0,0]) %%[0,0,-1]是第一次旋转的转轴 pause(0.1)end %%旋转动画 for k=1:N alpha(k)=a1*((N-k+1)^2*t^2-(N-k)^2*t^2)/2; rotate(H,[0,0,-1],180/pi*alpha(k),[0,0,0]) pause(0.1)end tt=Time2/N;for k=1:N beta(k)=a2*(k^2*tt^2-(k-1)^2*tt^2)/2; rotate(H,[yn*zp-zn*yp,zn*xp-xn*zp,xn*yp-yn*xp],180/pi*beta(k),[0,0,0]) %%[yn*zp-zn*yp,xn*zp-zn*xp,xn*yp-yn*xp]是第二次旋转的转轴 pause(0.1)endfor k=1:N beta(k)=a2*((N-k+1)^2*tt^2-(N-k)^2*tt^2)/2; rotate(H,[yn*zp-zn*yp,zn*xp-xn*zp,xn*yp-yn*xp],180/pi*beta(k),[0,0,0]) pause(0.1)end模型二：animation2函数
function animation1(action)cla resetxa=0;ya=1;za=0;xb=sqrt(3)/2;yb=-1/2;zb=0;xc=-sqrt(3)/2;yc=-1/2;zc=0; %%A,B,C坐标xn=0;yn=0;zn=-1; %%初始卫星所对准的点Na=str2num(get(findobj(gcf,'tag','h_e1'),'string'));b=str2num(get(findobj(gcf,'tag','h_e2'),'string'));c=str2num(get(findobj(gcf,'tag','h_e3'),'string'));xp=a;yp=b;zp=c; %%要求卫星所对准的点Palpha1=acos((xa*xp+ya*yp)/(sqrt(xa^2+ya^2)*sqrt(xp^2+yp^2))); alpha2=acos((xb*xp+yb*yp)/(sqrt(xb^2+yb^2)*sqrt(xp^2+yp^2))); alpha3=acos((xc*xp+yc*yp)/(sqrt(xc^2+yc^2)*sqrt(xp^2+yp^2))); A=[alpha1; alpha2; alpha3];alpha=min(A,[],1); %%计算alpha beta=acos((xn*xp+yn*yp+zn*zp)/sqrt(xn^2+yn^2+zn^2)/sqrt(xp^2+yp^2+zp^2)); %%计算beta v=50; %%喷气对卫星本身的相对速度,单位是m/s Q=0.004*v; %%喷气喷量最大值,单位是Kg/s J1=30000; %%卫星对于过其质心一轴的转动惯量,单位是Kg*m^2 J2=20000; %%卫星对于过其质心另一轴的转动惯量，单位是Kg*m^2 r=1; %%6个喷嘴到卫星转动轴的距离,单位是m l=10; %%上下两组喷嘴平面之间的距离，单位是m %-------------------------------------------------------------------------------------aa=(cos(alpha)+cos(pi/3-alpha)+(sin(pi/3-alpha)-sin(alpha))*cos(pi/3+alpha)/sin(pi/3+alpha))*Q*v*l/J1;Time=sqrt(beta/aa); %%计算加速减速总时间 %------------------------------------------------------------------------------------------H0=plot3([0,0],[0,0],[-4,4],'m'); %%画卫星实体hold on[h1x,h1y,h1z]=ellipsoid(0,0,0,0.5,0.5,0.5,40); H1=[h1x,h1y,h1z];surfl(h1x,h1y,h1z);shading interp;hold on[h2x,h2y,h2z]=ellipsoid(0,0,0,0.75,0.75,0.25,40); H2=[h2x,h2y,h2z];surfl(h2x,h2y,h2z);shading flat;hold onH4=plot3([0,0],[0,1],[0,0],'r--',... [0,sqrt(3)/2],[0,-1/2],[0,0],'r--',[0,-sqrt(3)/2],[0,-1/2],[0,0],'r--'); %%画OA,OB,OC hold on H5=plot3([0,0],[0,1],[0.75,0.75],'b',... [0,sqrt(3)/2],[0,-1/2],[0.75,0.75],'b',[0,-sqrt(3)/2],[0,-1/2],[0.75,0.75],'b'); %%画OA1,OB1,OC1hold onH6=plot3([0,0],[0,1],[-0.75,-0.75],'b',... [0,sqrt(3)/2],[0,-1/2],[-0.75,-0.75],'b',[0,-sqrt(3)/2],[0,-1/2],[-0.75,-0.75],'b'); %%画OA2,OB2,OC2 H=get(gca,'children'); %%获得当前图形对象的句柄hold onH3=plot3([-1.5,1.5],[0,0],[0,0],'k ',... [0,0],[-1.5,1.5],[0,0],'k',... [0,0],[0,0],[-2.5,2.5],'k'); text(1.7,0,0,'x','fontsize',14); text(0,1.7,0,'y','fontsize',14); text(0,0,2.7,'z','fontsize',14); %%画x,y,z轴hold onplot3(a,b,c,'r+'); text(a+0.1,b+0.1,c+0.1,'P','fontsize',14); %%标定P点 %%标定P点 N=30;t=Time/N;for k=1:N beta(k)=aa*(k^2*t^2-(k-1)^2*t^2)/2; rotate(H,[yn*zp-zn*yp,zn*xp-xn*zp,xn*yp-yn*xp],180/pi*beta(k),[0,0,0]) %%[yn*zp-zn*yp,xn*zp-zn*xp,xn*yp-yn*xp]是旋转的转轴 pause(0.1)endfor k=1:N beta(k)=aa*((N-k+1)^2*t^2-(N-k)^2*t^2)/2; rotate(H,[yn*zp-zn*yp,zn*xp-xn*zp,xn*yp-yn*xp],180/pi*beta(k),[0,0,0]) pause(0.1)end绘图：curve函数
function curve(action)xa=0;ya=1;za=0;xb=sqrt(3)/2;yb=-1/2;zb=0;xc=-sqrt(3)/2;yc=-1/2;zc=0; %%A,B,C坐标xn=0;yn=0;zn=-1; %%初始卫星所对准的点Na=str2num(get(findobj(gcf,'tag','h_e1'),'string'));b=str2num(get(findobj(gcf,'tag','h_e2'),'string'));c=str2num(get(findobj(gcf,'tag','h_e3'),'string'));xp=a;yp=b;zp=c; %%要求卫星所对准的点Palpha1=acos((xa*xp+ya*yp)/(sqrt(xa^2+ya^2)*sqrt(xp^2+yp^2))); alpha2=acos((xb*xp+yb*yp)/(sqrt(xb^2+yb^2)*sqrt(xp^2+yp^2))); alpha3=acos((xc*xp+yc*yp)/(sqrt(xc^2+yc^2)*sqrt(xp^2+yp^2))); A=[alpha1; alpha2; alpha3];alpha=min(A,[],1); %%计算alpha beta=acos((xn*xp+yn*yp+zn*zp)/sqrt(xn^2+yn^2+zn^2)/sqrt(xp^2+yp^2+zp^2)); %%计算beta v=50; %%喷气对卫星本身的相对速度,单位是m/s Q=0.004*v; %%喷气喷量最大值,单位是Kg/s J1=30000; %%卫星对于过其质心一轴的转动惯量,单位是Kg*m^2 J2=20000; %%卫星对于过其质心另一轴的转动惯量，单位是Kg*m^2 r=1; %%6个喷嘴到卫星转动轴的距离,单位是m l=10; %%上下两组喷嘴平面之间的距离，单位是m a1=6*Q*v*r/J2;Time1=sqrt(alpha/a1); %%计算第一次加速减速总时间 a2=2*Q*v*l/J1;Time2=sqrt(beta/a2); %%计算第二次加速减速总时间 aa=(cos(alpha)+cos(pi/3-alpha)+(sin(pi/3-alpha)-sin(alpha))*cos(pi/3+alpha)/sin(pi/3+alpha))*Q*v*l/J1;Time=sqrt(beta/aa); %%计算加速减速总时间 t11=linspace(0,Time1,100);t12=linspace(Time1,2*Time1,100);y11=a1.*t11;y12=2*a1*Time1-a1.*t12;t21=linspace(0,Time2,100);t22=linspace(Time2,2*Time2,100);y21=a2.*t21;y22=2*a2*Time2-a2.*t22;beta11=a2.*t21.^2/2;beta12=a2*Time2^2/2+a2*Time2.*(t22-Time2)-a2.*(t22-Time2).^2/2;tt1=linspace(0,Time,100);tt2=linspace(Time,2*Time,100);yy1=aa.*tt1;yy2=2*aa*Time-aa.*tt2;beta21=aa.*tt1.^2/2;beta22=aa*Time^2/2+aa*Time.*(tt2-Time)-aa.*(tt2-Time).^2/2;figure set(gcf,'numbertitle','off','name','角速度与转角曲线图');set(gcf,'defaultuicontrolfontsize',11,'color',[0.4 0.8 0.6]); uicontrol(gcf,'style','text',... 'unit','normalized',... 'position',[0.24 0.94 0.1 0.03],... 'string','模型一',... 'fontsize',10); uicontrol(gcf,'style','text',... 'unit','normalized',... 'position',[0.70 0.64 0.1 0.03],... 'string','模型二',... 'fontsize',10); subplot(3,2,1) plot(t11,y11,'r',t12,y12,'r') ylabel('自转角速度') xlabel('时间') hold on subplot(3,2,3) plot(2*Time1,0,'k',2*Time1+t21,y21,'k',2*Time1+t22,y22,'k') ylabel('角速度') xlabel('时间') hold on subplot(3,2,5) plot(2*Time1,0,'r',2*Time1+t21,beta11,'r',2*Time1+t22,beta12,'r') ylabel('转角') xlabel('时间') hold on subplot(3,2,4) plot(tt1,yy1,'k',tt2,yy2,'k') ylabel('角速度') xlabel('时间') hold on subplot(3,2,6) plot(tt1,beta21,'r',tt2,beta22,'r') ylabel('转角') xlabel('时间')主程序：包括GUI设计和程序运行
%----------------------------------------------------------------------------------figureset(gcf,'numbertitle','off','name','卫星调整方向仿真');set(gcf,'defaultuicontrolfontsize',11,'color',[0.4 0.8 0.96]);h_t1=uicontrol(gcf,'style','text',... 'unit','normalized',... 'position',[0.8 0.85 0.15 0.06],... 'string','P点x坐标值',... 'fontsize',12,... 'backgroundcolor',[0.4 0.8 0.96]); h_t2=uicontrol(gcf,'style','text',... 'unit','normalized',... 'position',[0.8 0.68 0.15 0.06],... 'string','P点y坐标值',... 'fontsize',12,... 'backgroundcolor',[0.4 0.8 0.96]); h_t3=uicontrol(gcf,'style','text',... 'unit','normalized',... 'position',[0.8 0.51 0.15 0.06],... 'string','P点z坐标值',... 'fontsize',12,... 'backgroundcolor',[0.4 0.8 0.96]);h_e1=uicontrol(gcf,'style','edit',... 'unit','normalized',... 'position',[0.80 0.80,0.1 0.05],... 'string','1.5',... 'tag','h_e1',... 'BackgroundColor',[1 1 0]);h_e2=uicontrol(gcf,'style','edit',... 'unit','normalized',... 'position',[0.80 0.63,0.1 0.05],... 'string','2.0',... 'tag','h_e2',... 'BackgroundColor',[1 1 0]); h_e3=uicontrol(gcf,'style','edit',... 'unit','normalized',... 'position',[0.80 0.47,0.1 0.05],... 'string','2.5',... 'tag','h_e3',... 'BackgroundColor',[1 1 0]); %---------------------------------------------------------------btnwidth=0.15;btnheight=0.09;btnnumber=4; yposition=(btnnumber-1)*btnheight; xposition=0; btnPos=[xposition yposition btnwidth btnheight]; hrotate=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position',btnPos, ... 'fontsize',12,... 'String','3D旋转', ... 'Callback','rotate3d'); btnnumber=3; yposition=(btnnumber-1)*btnheight; xposition=0; btnPos=[xposition yposition btnwidth btnheight]; hzoomin=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position',btnPos, ... 'fontsize',12,... 'String','放大', ... 'Callback','zoom(1.1)'); btnnumber=2; yposition=(btnnumber-1)*btnheight; xposition=0; btnPos=[xposition yposition btnwidth btnheight]; hzoomout=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position',btnPos, ... 'fontsize',12,... 'String','缩小', ... 'Callback','zoom(0.9)'); btnnumber=1; yposition=(btnnumber-1)*btnheight; xposition=0; btnPos=[xposition yposition btnwidth btnheight]; hclose=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position',btnPos, ... 'fontsize',12,... 'String','关闭', ... 'Callback','close(gcf)'); % ---------------------------------------------------------------------- btnnumber=4; yposition=(btnnumber-1)*btnheight; xposition=1-btnwidth; btnPos=[xposition yposition btnwidth btnheight]; hintegration=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position', btnPos, ... 'fontsize',12,... 'String','曲线图', ... 'Callback','curve'); btnnumber=3; yposition=(btnnumber-1)*btnheight; xposition=1-btnwidth; btnPos=[xposition yposition btnwidth btnheight]; hmodel1=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position', btnPos, ... 'fontsize',12,... 'String','运行模型一', ... 'Callback','animation1'); btnnumber=2; yposition=(btnnumber-1)*btnheight; xposition=1-btnwidth; btnPos=[xposition yposition btnwidth btnheight]; hmodel2=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position', btnPos, ... 'fontsize',12,... 'String','运行模型二', ... 'Callback','animation2'); btnnumber=1; yposition=(btnnumber-1)*btnheight; xposition=1-btnwidth; btnPos=[xposition yposition btnwidth btnheight]; hedit=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position', btnPos, ... 'fontsize',12,... 'String','源程序', ... 'Callback','edit lllxgui');——————华丽分割线——————
把主程序代码存为一个MABLAB脚本文件(如maingui.m)，与其它3个函数的m文件放在一个目录，在MATLAB命令行窗口直接运行maingui即可。
GUI界面是这样的：
很帅对不对！！很有美感对不对！！∠( °ω°)／
接着点击“运行模型一”或“运行模型二”可以分别运行两个模型，点击“曲线图”可以查看观察变量曲线，其它按钮如字面意思。
运行模型一：
运行模型二：
查看变量变化曲线：
其它功能就不赘述了。
——————华丽分割线——————
最后的最后，我忍不住想要吐槽一下，现在我经常批学生提交的文档和制作的PPT不规范且太丑，当我看到我当年的作业后，我选择了……闭嘴……ヾ(*ΦωΦ)ツヾ(*ΦωΦ)ツヾ(*ΦωΦ)ツ
扫一扫
8863
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
