1、wz0=20; %偏航轴初始旋转角速度tf=2*pi/（Ix/It-1）*wx0）;ts=0,tf;w0=wx0;wy0;wz0; %初始条件t,w=ode45（dwdt,ts,w0）;m,n=size（w）;plot3（w（:,2）,w（:,3）,w（:,1）;axis（-30,30,-30,30,0,60）;hold ongg（w）;theta=atan（It/Ix）*（sqrt（wy02+wz02）/wx0）;%beta=atan（w（:,3）./w（:,2）;beta=0.35;gama=atan（sqrt（wy02+wz02）/wx0）;k1=tan（theta-gama）/tan（
2、gama）;k2=cos（theta-gama）/cos（gama）;H（:,1）=k2*w（:,1）;,2）=k2*k1*w（:,2）;,3）=k2*k1*w（:,3）;gg1（H）;xlabel（z）,ylabel（y）,zlabel（x） grid on,hold onx=0:80;y=0*x;z=0*x;x1=0:y1=0*x1;%z=（sqrt（wy02+wz02）/wx0）*x;zz=（It/Ix）*（sqrt（wy02+wz02）/wx0）*x1;%plot3（-z,y,x）;%plot3（y1,y1,x1,-）;%plot3（y,z,x）;%plot3（y,-z,x）;plot
3、3（y,z,x,m-,Linewidthplot3（y1,zz,x1,hold on %动量矩H的方向text（0,10,70）,text（30,0,70）,text（-10,10,60）,hold off其中包含如下函数体：（1）function omegad=dwdt（t,w）global Ix It wx0 OmegaOmega=（Ix-It）*wx0/It;wx0=w（1）;wy=w（2）;wz=w（3）;omegad=0;-Omega*wz;Omega*wy;（2）function gg（a）gama=atan（sqrt（a（1,2）2+a（1,3）2）/a（1,1）;for i=0
4、:0.1:1 b=i*a; plot3（b（:,2）,b（:,3）,b（:endm,n=size（a）;for j=1:30 k=fix（1+（j-1）*0.5*m/30）; alpha=atan（a（k,2）/a（k,3）; x=0:a（1,1）;y=sin（alpha）*tan（gama）*x; z=cos（alpha）*tan（gama）*x; plot3（y,z,x）; plot3（-y,-z,x）;（3）function gg1（a）A=cos（-theta）,0,-sin（-theta）; 0,1,0; sin（-theta）,0,cos（-theta）; %转换矩阵 b=i*a*A;,1）,g XX（:,1）=x;XX（:,2）=y,3）=z YY（:YY（:,2）=-y,3）=-z XX=XX*A; YY=YY*A; plot3（XX（:,2）,XX（:,3）,XX（:c plot3（YY（:,2）,YY（:,3）,YY（: 逆进动2、 正进动Ix=60;It=80; 正进动
