一个用matlab编的可以动的飞机动画制作小程序

MVH

一个用matlab编的可以动的飞机动画制作小程序，注意六个状态值是需要外部输入的，比较适合做飞行器仿真试验用。

该程序可以运行，但是尚存在一定的问题，希望大家讨论解决


总共有8个m文件

PLANE_.M

1. % plane_ same as List_B4.
   
2. % Plots a wire frame airplane.
   
3. % Copyright S. Nakamura
   
4. set(gcf, 'NumberTitle','off','Name', 'Figure B.3; plane_ ')
   
5. 
   
6. 
   
7. clc
   
8. clear, clf
   
9. 
   
10. 
    
11. mov1 = avifile('airplane.avi');
    
12. mov2 = avifile('airplane_xz.avi');
    
13. mov3 = avifile('airplane_xy,avi');
    
14. mov4 = avifile('airplane_yz,avi');
    
15. 
    
16. for(hhh = 0:-1:-15)
    
17.     clf
    
18.    
    
19.    
    
20. dth=pi/16;
    
21. fuselen=6;
    
22. thf=pi:-dth:pi/2;
    
23. xa = 0:0.5:fuselen;
    
24. 
    
25. hyq_x = hhh;
    
26. hyq_y = hhh;
    
27. hyq_z = -hhh;
    
28. 
    
29. xt=fuselen+0.25:0.25:fuselen+2;
    
30. dxt = 1.4/(length(xt)-0) ;
    
31. yt = -1+dxt:dxt  :0.4;
    
32. length(yt);
    
33. 
    
34. 
    
35. xft=[cos(thf),xa, xt] ;
    
36. 
    
37. yft=[sin(thf)- 0.3*sin(2*thf).^4, ones(size(xa)), ones(size(yt))] ;
    
38. 
    
39. 
    
40. xfb=[cos(thf),xa,xt] ;
    
41. 
    
42. yfb=[-sin(thf),-ones(size(xa)),yt ] ;
    
43. k=length([thf, xa]);
    
44. yfb(k)=( yfb(k-1)+yfb(k+1))/2;
    
45. xc =(xfb+xft)/2 ;
    
46. yc = (yfb+yft)/2;
    
47. L=length(xc);
    
48. for i=1:L
    
49. if xc(i)<0
    
50.    yc(i)=0;
    
51. end
    
52. end
    
53. 
    
54. 
    
55. 
    
56. 
    
57. % plot(xfb,yfb,'y', xft,yft,'y', xc, yc,':')
    
58. %
    
59. % axis([-2 8 -6 6])
    
60. % pause
    
61. 
    
62. 
    
63. a=0.5;
    
64. b=0.5;
    
65. 
    
66. dth=pi/8;
    
67. 
    
68. th=0:dth:2*pi;
    
69. jmax=length(th);
    
70. xr=cos(th);
    
71. 
    
72. yr=sin(th);
    
73. L=length(xc);
    
74. 
    
75. for i=1:L
    
76. %   xr=cos(th);
    
77. %   yr=sin(th);
    
78.   a = (yft(i)    -yc(i))/(-yfb(i) + yc(i));
    
79.   b = (-yfb(i)    +yc(i));
    
80.   for j=1:jmax
    
81.        y(i,j)=yr(j)*b+yc(i);
    
82.        if th(j)<pi    y(i,j)=yr(j)*b*a + yc(i); end
    
83.     x(i,j)=xr(j)*b;
    
84.     z(i,j)=xc(i);
    
85.   end
    
86. end
    
87. opengl neverselect
    
88. surface(z+hyq_x,x+hyq_y,y+hyq_z)
    
89. axis([-8 8 -16 16 -16 16])
    
90.   %mesh(z(1:30,1:11),x(1:30,1:11),y(1:30,1:11))
    
91.    hold on
    
92.    surfc(z+hyq_x,x+hyq_y,y+hyq_z)
    
93. [xw,yw,zw] = wing_;
    
94. F = 1.7;
    
95. xw=F*xw; yw=F*yw; zw=F*zw;
    
96. [x1,y1,z1] = rotz_(xw,yw,zw,-90);
    
97. [x2,y2,z2] = rotx_(xw,yw,zw,180);
    
98. [x2,y2,z2] = rotz_(x2,y2,z2,-270);
    
99. surface(x1+2+hyq_x,y1-0.5+hyq_y,z1+ 0.7+hyq_z);
    
100. surface(x2+2+hyq_x,y2+0.5+hyq_y,z2+ 0.7+hyq_z);
     
101. surface(0.8*x1+6.6+hyq_x,0.5*z1-0+hyq_y,-0.3*y1+1.2+hyq_z)
     
102. surface(0.7*x1+6.6+hyq_x,0.3*y1-0.7+hyq_y,0.9*z1+ 0.7+hyq_z);
     
103. surface(0.7*x2+6.6+hyq_x,0.3*y2+0.7+hyq_y,0.9*z2+ 0.7+hyq_z);
     
104. caxis([-3,1])
     
105. grid
     
106. %[x3,y3,z3] = rotx_(x1,y1,z1,+80)
     
107. 
     
108. %axis([-2 8 -6 6 -6 6])
     
109.     for(jjjj = 1:10)
     
110.     view(3)
     
111.     F = getframe(gca);
     
112.     mov1 = addframe(mov1,F);
     
113.     view([0 0])
     
114.     F = getframe(gca);
     
115.     mov2 = addframe(mov2,F);
     
116.     view([0 90])
     
117.     F = getframe(gca);
     
118.     mov3 = addframe(mov3,F);
     
119.     view([0 90])
     
120.     mov4 = addframe(mov4,F);
     
121.     end
     
122. 
     
123. title('Commuter Airplane')
     
124. caxis([-2, 2])
     
125. colormap(hsv)
     
126. end
     
127. 
     
128. mov1 = close(mov1);
     
129. mov2 = close(mov2);
     
130. mov3 = close(mov3);
     
131. mov4 = close(mov4);

复制代码

[ 本帖最后由 lxq 于 2007-1-1 17:08 编辑 ]

MVH

PLANE1_.m

1. % plane_ same as List_B4.
   
2. % Plots a wire frame airplane.
   
3. % Copyright S. Nakamura
   
4. set(gcf, 'NumberTitle','off','Name', 'Figure B.3; plane_ ')
   
5. %  scrsz = get(0,'ScreenSize');
   
6. %  set(gcf,'Position',[30 80 scrsz(3)-280 scrsz(4)-280])
   
7. 
   
8. clc
   
9. clear all, clf
   
10. 
    
11. mov1 = avifile('airplanh.avi');
    
12. for hhh = 0:0.2:50
    
13. dth=pi/16;
    
14. fuselen=6;
    
15. thf=pi:-dth:pi/2;
    
16. xa = 0:0.5:fuselen;
    
17. 
    
18. hyq_x = 0;
    
19. hyq_y = 0;
    
20. hyq_z = 0;
    
21. 
    
22. xt=fuselen+0.25:0.25:fuselen+2;
    
23. dxt = 1.4/(length(xt)-0) ;
    
24. yt = -1+dxt:dxt  :0.4;
    
25. length(yt);
    
26. 
    
27. 
    
28. xft=[cos(thf),xa, xt] ;
    
29. 
    
30. yft=[sin(thf)- 0.3*sin(2*thf).^4, ones(size(xa)), ones(size(yt))] ;
    
31. 
    
32. 
    
33. xfb=[cos(thf),xa,xt] ;
    
34. 
    
35. yfb=[-sin(thf),-ones(size(xa)),yt ] ;
    
36. k=length([thf, xa]);
    
37. yfb(k)=( yfb(k-1)+yfb(k+1))/2;
    
38. xc =(xfb+xft)/2 ;
    
39. yc = (yfb+yft)/2;
    
40. L=length(xc);
    
41. for i=1:L
    
42. if xc(i)<0
    
43.    yc(i)=0;
    
44. end
    
45. end
    
46. 
    
47. 
    
48. 
    
49. 
    
50. % plot(xfb,yfb,'y', xft,yft,'y', xc, yc,':')
    
51. %
    
52. % axis([-2 8 -6 6])
    
53. % pause
    
54. 
    
55. 
    
56. a=0.5;
    
57. b=0.5;
    
58. 
    
59. dth=pi/8;
    
60. 
    
61. th=0:dth:2*pi;
    
62. jmax=length(th);
    
63. xr=cos(th);
    
64. 
    
65. yr=sin(th);
    
66. L=length(xc);
    
67. 
    
68. for i=1:L
    
69. 
    
70.   a = (yft(i)    -yc(i))/(-yfb(i) + yc(i));
    
71.   b = (-yfb(i)    +yc(i));
    
72.   for j=1:jmax
    
73.        y(i,j)=yr(j)*b+yc(i);
    
74.        if th(j)<pi    y(i,j)=yr(j)*b*a + yc(i); end
    
75.     x(i,j)=xr(j)*b;
    
76.     z(i,j)=xc(i);
    
77.   end
    
78. end
    
79. opengl neverselect
    
80. 
    
81. 
    
82. 
    
83. [xw,yw,zw] = wing_;
    
84. F = 1.7;
    
85. xw=F*xw; yw=F*yw; zw=F*zw;
    
86. [x1,y1,z1] = rotz_(xw,yw,zw,-90);
    
87. [x2,y2,z2] = rotx_(xw,yw,zw,180);
    
88. [x2,y2,z2] = rotz_(x2,y2,z2,-270);
    
89. 
    
90. aa = -2.8;
    
91. z = z+aa;
    
92. x3 = x1+2+aa;
    
93. y3 = y1-0.5;
    
94. z3 = z1+0.7;
    
95. x4 = x2+2+aa;
    
96. y4 = y2+0.5;
    
97. z4 = z2+0.7;
    
98. x5 = 0.8*x1+6.6+aa;
    
99. z5 = 0.5*z1;
    
100. y5 = -0.3*y1+1.2;
     
101. x6 = 0.7*x1+6.6+aa;
     
102. y6 = 0.3*y1-0.7;
     
103. z6 = 0.9*z1+0.7;
     
104. x7 = 0.7*x2+6.6+aa;
     
105. y7 = 0.3*y2+0.7;
     
106. z7 = 0.9*z2+0.7;
     
107. 
     
108. % 侧倾角变换
     
109. [x3,y3,z3]=rotx_(x3,y3,z3,hhh);
     
110. [x4,y4,z4]=rotx_(x4,y4,z4,hhh);
     
111. [x5,z5,y5]=rotx_(x5,z5,y5,hhh);
     
112. [x6,y6,z6]=rotx_(x6,y6,z6,hhh);
     
113. [x7,y7,z7]=rotx_(x7,y7,z7,hhh);
     
114. [z,x,y] = rotx_(z,x,y,hhh);
     
115. 
     
116. % 俯仰角变换
     
117. [x3,y3,z3]=roty_(x3,y3,z3,hhh);
     
118. [x4,y4,z4]=roty_(x4,y4,z4,hhh);
     
119. [x5,z5,y5]=roty_(x5,z5,y5,hhh);
     
120. [x6,y6,z6]=roty_(x6,y6,z6,hhh);
     
121. [x7,y7,z7]=roty_(x7,y7,z7,hhh);
     
122. [z,x,y] = roty_(z,x,y,hhh);
     
123. 
     
124. % 航向角变换
     
125. [x3,y3,z3]=rotz_(x3,y3,z3,hhh);
     
126. [x4,y4,z4]=rotz_(x4,y4,z4,hhh);
     
127. [x5,z5,y5]=rotz_(x5,z5,y5,hhh);
     
128. [x6,y6,z6]=rotz_(x6,y6,z6,hhh);
     
129. [x7,y7,z7]=rotz_(x7,y7,z7,hhh);
     
130. [z,x,y] = rotz_(z,x,y,hhh);
     
131. 
     
132. clf
     
133. subplot(221)
     
134. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
135. axis([-16 16 -16 16 -16 16])
     
136. 
     
137. hold on
     
138. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
139. 
     
140. 
     
141. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
142. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
143. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
144. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
145. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
146. grid
     
147. 
     
148. 
     
149. subplot(222)
     
150. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
151. axis([-16 16 -16 16 -16 16])
     
152. 
     
153. hold on
     
154. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
155. 
     
156. 
     
157. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
158. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
159. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
160. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
161. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
162. view([0 0])
     
163. grid
     
164. 
     
165. 
     
166. subplot(223)
     
167. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
168. axis([-16 16 -16 16 -16 16])
     
169. 
     
170. hold on
     
171. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
172. 
     
173. 
     
174. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
175. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
176. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
177. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
178. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
179. view([0 90])
     
180. grid
     
181. 
     
182. 
     
183. subplot(224)
     
184. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
185. axis([-16 16 -16 16 -16 16])
     
186. 
     
187. hold on
     
188. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
189. 
     
190. 
     
191. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
192. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
193. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
194. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
195. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
196. view([90 0])
     
197. grid
     
198. 
     
199. 
     
200. caxis([-3,1])
     
201. 
     
202. title('Commuter Airplane')
     
203. caxis([-2, 2])
     
204. colormap(hsv)
     
205. F = getframe(gcf);
     
206. mov1 = addframe(mov1,F);
     
207. end
     
208. 
     
209. mov1 = close(mov1);

复制代码

MVH

PLANE2_.m

1. function PLANT2_(north, east, height, yaw, pitch, roll)
   
2. 
   
3. % plane_ same as List_B4.
   
4. % Plots a wire frame airplane.
   
5. % Copyright S. Nakamura
   
6. % Updated By: Yuing He
   
7. 
   
8. set(gcf, 'NumberTitle','off','Name', 'Figure B.3; plane_ ','Position',[10 40 1000 650])
   
9. 
   
10. clc
    
11. 
    
12. dth=pi/16;
    
13. fuselen=6;
    
14. thf=pi:-dth:pi/2;
    
15. xa = 0:0.5:fuselen;
    
16. 
    
17. hyq_x = north;
    
18. hyq_y = east;
    
19. hyq_z = height;
    
20. 
    
21. xt=fuselen+0.25:0.25:fuselen+2;
    
22. dxt = 1.4/(length(xt)-0) ;
    
23. yt = -1+dxt:dxt  :0.4;
    
24. length(yt);
    
25. 
    
26. 
    
27. xft=[cos(thf),xa, xt] ;
    
28. yft=[sin(thf)- 0.3*sin(2*thf).^4, ones(size(xa)), ones(size(yt))] ;
    
29. 
    
30. xfb=[cos(thf),xa,xt] ;
    
31. yfb=[-sin(thf),-ones(size(xa)),yt ] ;
    
32. k=length([thf, xa]);
    
33. yfb(k)=( yfb(k-1)+yfb(k+1))/2;
    
34. 
    
35. xc =(xfb+xft)/2 ;
    
36. yc = (yfb+yft)/2;
    
37. 
    
38. L=length(xc);
    
39. for i=1:L
    
40. if xc(i)<0
    
41.    yc(i)=0;
    
42. end
    
43. end
    
44. 
    
45. dth=pi/8;
    
46. 
    
47. th=0:dth:2*pi;
    
48. jmax=length(th);
    
49. xr=cos(th);
    
50. 
    
51. yr=sin(th);
    
52. L=length(xc);
    
53. 
    
54. for i=1:L
    
55.   a = (yft(i) - yc(i))/(-yfb(i) + yc(i));
    
56.   b = (-yfb(i) + yc(i));
    
57.   for j=1:jmax
    
58.     y(i,j)=yr(j)*b+yc(i);
    
59.     if th(j)<pi    y(i,j)=yr(j)*b*a + yc(i); end
    
60.     x(i,j)=xr(j)*b;
    
61.     z(i,j)=xc(i);
    
62.   end
    
63. end
    
64. 
    
65. [xw,yw,zw] = wing_;
    
66. F = 1.7;
    
67. xw=F*xw; yw=F*yw; zw=F*zw;
    
68. [x1,y1,z1] = rotz_(xw,yw,zw,-90);
    
69. [x2,y2,z2] = rotx_(xw,yw,zw,180);
    
70. [x2,y2,z2] = rotz_(x2,y2,z2,-270);
    
71. 
    
72. aa = -2.8;
    
73. z = z+aa;
    
74. x3 = x1+2+aa;
    
75. y3 = y1-0.5;
    
76. z3 = z1+0.7;
    
77. x4 = x2+2+aa;
    
78. y4 = y2+0.5;
    
79. z4 = z2+0.7;
    
80. x5 = 0.8*x1+6.6+aa;
    
81. z5 = 0.5*z1;
    
82. y5 = -0.3*y1+1.2;
    
83. x6 = 0.7*x1+6.6+aa;
    
84. y6 = 0.3*y1-0.7;
    
85. z6 = 0.9*z1+0.7;
    
86. x7 = 0.7*x2+6.6+aa;
    
87. y7 = 0.3*y2+0.7;
    
88. z7 = 0.9*z2+0.7;
    
89. 
    
90. % 横滚角变换
    
91. [x3,y3,z3]=rotx_(x3,y3,z3,roll);
    
92. [x4,y4,z4]=rotx_(x4,y4,z4,roll);
    
93. [x5,z5,y5]=rotx_(x5,z5,y5,roll);
    
94. [x6,y6,z6]=rotx_(x6,y6,z6,roll);
    
95. [x7,y7,z7]=rotx_(x7,y7,z7,roll);
    
96. [z,x,y] = rotx_(z,x,y,roll);
    
97. 
    
98. % 俯仰角变换
    
99. [x3,y3,z3]=roty_(x3,y3,z3,pitch);
    
100. [x4,y4,z4]=roty_(x4,y4,z4,pitch);
     
101. [x5,z5,y5]=roty_(x5,z5,y5,pitch);
     
102. [x6,y6,z6]=roty_(x6,y6,z6,pitch);
     
103. [x7,y7,z7]=roty_(x7,y7,z7,pitch);
     
104. [z,x,y] = roty_(z,x,y,pitch);
     
105. 
     
106. % 航向角变换
     
107. [x3,y3,z3]=rotz_(x3,y3,z3,yaw);
     
108. [x4,y4,z4]=rotz_(x4,y4,z4,yaw);
     
109. [x5,z5,y5]=rotz_(x5,z5,y5,yaw);
     
110. [x6,y6,z6]=rotz_(x6,y6,z6,yaw);
     
111. [x7,y7,z7]=rotz_(x7,y7,z7,yaw);
     
112. [z,x,y] = rotz_(z,x,y,yaw);
     
113. 
     
114. opengl neverselect
     
115. clf
     
116. subplot(221)
     
117. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
118. axis([north-20 north+20 east-20 east+20 height-20 height+20])
     
119. 
     
120. hold on
     
121. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
122. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
123. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
124. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
125. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
126. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
127. grid
     
128. title('3 dimention view');
     
129. 
     
130. subplot(222)
     
131. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
132. axis([north-20 north+20 east-20 east+20 height-20 height+20])
     
133. 
     
134. hold on
     
135. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
136. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
137. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
138. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
139. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
140. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
141. view([0 0])
     
142. grid
     
143. title('looking throgh negative x');
     
144. 
     
145. 
     
146. subplot(223)
     
147. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
148. axis([north-20 north+20 east-20 east+20 height-20 height+20])
     
149. 
     
150. hold on
     
151. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
152. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
153. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
154. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
155. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
156. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
157. view([0 90])
     
158. grid
     
159. title('looking throgh negative y');
     
160. 
     
161. 
     
162. subplot(224)
     
163. surface(z+hyq_x,x+hyq_y,y+hyq_z)
     
164. axis([north-20 north+20 east-20 east+20 height-20 height+20])
     
165. 
     
166. hold on
     
167. surfc(z+hyq_x,x+hyq_y,y+hyq_z)
     
168. surface(x3+hyq_x,y3+hyq_y,z3+hyq_z);
     
169. surface(x4+hyq_x,y4+hyq_y,z4+hyq_z);
     
170. surface(x5+hyq_x,z5+hyq_y,y5+hyq_z)
     
171. surface(x6+hyq_x,y6+hyq_y,z6+hyq_z);
     
172. surface(x7+hyq_x,y7+hyq_y,z7+hyq_z);
     
173. view([90 0])
     
174. grid
     
175. title('looking throgh negative z');
     
176. 
     
177. 
     
178. caxis([-3,1])
     
179. caxis([-2,2])
     
180. colormap(hsv)
     
181. 
     
182. end

复制代码

MVH

move.m

1. mov1 = avifile('airplan2.avi');
   
2. 
   
3. length1=length(height);
   
4. 
   
5. for hh = 1:1:length1/10
   
6.     hhh = hh*10;
   
7.     north1 =- north(hhh) ;
   
8.     east1 = east(hhh);
   
9.     height1 = height(hhh);
   
10.     yaw = psid(hhh);
    
11.     pitch = thetad(hhh);
    
12.     roll = phid(hhh);
    
13.     plane2_(north1,east1,height1,yaw,pitch,roll);
    
14.     F = getframe(gcf);
    
15.     mov1 = addframe(mov1,F);
    
16.     mov1 = addframe(mov1,F);
    
17.     mov1 = addframe(mov1,F);
    
18. end
    
19. 
    
20. mov1 = close(mov1);

复制代码

MVH

wing_.M

1. % wing_ generates a wing data used in plane_
   
2. % (originally named wing_2d)
   
3. % Example:   [x,y,z]=wing_; mesh(x,y,z)
   
4. % See Figure B.3
   
5. function [zw,xw,yw] = wing_
   
6. x=0:0.1:1;
   
7. n=length(x);
   
8. for k=1:30
   
9. 
   
10. for i=2:n-1
    
11. x(i)=0.5*x(i-1)+0.4*x(i+1);
    
12. end
    
13. end
    
14. for i=2:n
    
15. x(n+i-1)=x(n-i+1);
    
16. end
    
17. y=0.2969*sqrt(x) - 0.126*x - 0.3516*x.^2  + ...
    
18.       0.2843 * x.^3 - 0.1015*x.^4;
    
19. for i=n+1:length(y)
    
20. y(i)=-y(i);
    
21. end
    
22. %axis([-0.5, 1.5 ,-1 1])
    
23. %plot(x,y)
    
24. jmax=15;
    
25. for j=1:jmax
    
26. for i=1:2*n-1
    
27. 
    
28. xw(i,j)=x(i);
    
29. yw(i,j)=y(i);
    
30. zw(i,j)=0.3*(j-1);
    
31. end
    
32. 
    
33. end
    
34. yw(:,jmax)=zeros(size(yw(:,jmax)));
    
35. zw(:,jmax)=zw(:,jmax-1);
    
36. 
    
37. %mesh(zw,xw,yw)
    
38. %axis([0,2,0,2,-1,1])

复制代码

MVH

rotx_.M

1. % rotx_(x,y,z,th) rotates a vector [x,y,z]
   
2. % th degrees counter-clockwise about
   
3. % the x-axis.  See Appendix B; List B.1.
   
4. function [xd,yd,zd]=rotx_(x,y,z,th)
   
5. cosf=cos(th*pi/180);sinf=sin(th*pi/180);
   
6. xd =x;
   
7. yd =  cosf.*y - sinf.*z;
   
8. zd =  sinf.*y + cosf.*z;

复制代码

MVH

roty_.M

1. % roty_(x,y,z,th) rotates a vector [x,y,z]
   
2. % th degrees counter-clockwise about
   
3. % the y-axis.  See Appendix B; List B.1.
   
4. function [xd,yd,zd]=roty_(x,y,z,th)
   
5. cosf=cos(th*pi/180);sinf=sin(th*pi/180);
   
6. yd =y;
   
7. xd =  cosf.*x + sinf.*z;
   
8. zd = - sinf.*x + cosf.*z;

复制代码

MVH

rotz_.M

1. % rotz_(x,y,z,th) rotates a vector [x,y,z]
   
2. % th degrees counter-clockwise about
   
3. % the z-axis.  See Appendix B; List B.1.
   
4. function [xd,yd,zd]=rotz_(x,y,z,phi)
   
5. cosf=cos(phi*pi/180);sinf=sin(phi*pi/180);
   
6. xd =  cosf *x - sinf *y;
   
7. yd = sinf *x + cosf *y;
   
8. zd =z;

复制代码

xjtu211

suffer

原帖由 xjtu211 于 2006-11-15 20:39 发表

该程序还有很多缺陷，希望大家讨论改进

y-y8641

bjydly

晕，论坛的程序为啥都没法运行呢？

suffer

原帖由 y-y8641 于 2006-11-20 15:12 发表

这个我试过，可以运行，不要没搞清楚就随便下判断

suffer

原帖由 bjydly 于 2006-11-21 14:59 发表
晕，论坛的程序为啥都没法运行呢？

可以运行的，只是有一点点小问题

hqhuang

是有点问题，运行后不会自动跳出。
要强制才能跳出，不过我是菜鸟一个不能解决。