|
进程被枪毙了.心里不痛快.把做了一半的matlab仿真程序拿出来供大家拍拍砖吧.写得比较乱.有不明白的可以回复或是发邮件给我.我的邮件是 phoenix8848@gmail.com
1 clear; 2 A=[0 0 0;0.65 0.25 0.376;1 0 0];R=0.25*1.414; 3 %calculate line vectors 4 v=ones(1,3);%initial line vectors 5 for i=1:2 6 v(i,:)=A(i+1,:)-A(i,:); 7 end 8 %calculate the interpolation angle 9 ang=acos(abs(v(1,1)*v(2,1)+v(1,2)*v(2,2)+v(1,3)*v(2,3))/(sqrt(v(1,1)^2+v(1,2)^2+v(1,3)^2)*sqrt(v(2,1)^2+v(2,2)^2+v(2,3)^2))); 10 %the dis pB to pAi 11 d=R*tan(ang/2); 12 %direction cosine of V(Ai-1,Ai) 13 cosa1=v(1,1)/sqrt(v(1,1)^2+v(1,2)^2+v(1,3)^2); 14 cosb1=v(1,2)/sqrt(v(1,1)^2+v(1,2)^2+v(1,3)^2); 15 cosc1=v(1,3)/sqrt(v(1,1)^2+v(1,2)^2+v(1,3)^2); 16 %interpolation start point 17 b1=[A(2,1)-d*cosa1,A(2,2)-d*cosb1,A(2,3)-d*cosc1]; 18 %direction cosine of V(Ai,Ai+1) 19 cosa2=v(2,1)/sqrt(v(2,1)^2+v(2,2)^2+v(2,3)^2); 20 cosb2=v(2,2)/sqrt(v(2,1)^2+v(2,2)^2+v(2,3)^2); 21 cosc2=v(2,3)/sqrt(v(2,1)^2+v(2,2)^2+v(2,3)^2); 22 %interpolation end point 23 b2=[A(2,1)+d*cosa2,A(2,2)+d*cosb2,A(2,3)+d*cosc2]; 24 %vector of angular bisector 25 v0=sqrt(v(1,1)^2+v(1,2)^2+v(1,3)^2); 26 v1=[v(1,1) v(1,2) v(1,3)]/v0; 27 v0=sqrt(v(2,1)^2+v(2,2)^2+v(2,3)^2); 28 v2=[v(2,1) v(2,2) v(2,3)]/v0; 29 L0=[(-v1(1)+v2(1))/2,(-v1(2)+v2(2))/2,(-v1(3)+v2(3))/2]; 30 %direction cosine of angular bisector 31 cosa3=L0(1)/sqrt(L0(1)^2+L0(2)^2+L0(3)^2); 32 cosb3=L0(2)/sqrt(L0(1)^2+L0(2)^2+L0(3)^2); 33 cosc3=L0(3)/sqrt(L0(1)^2+L0(2)^2+L0(3)^2); 34 %calculate the center point of interpolation circle 35 t=R*sec(ang/2); 36 p=[A(2,1)+cosa3*t,A(2,2)+cosb3*t,A(2,3)+cosc3*t]; 37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 38 %calculate the normal vector of interpolation circle 39 dx=[b1(1)-p(1) b1(2)-p(2) b1(3)-p(3)]; 40 dx=dx/sqrt(dx(1)^2+dx(2)^2+dx(3)^2); 41 dy=[b2(1)-p(1) b2(2)-p(2) b2(3)-p(3)]; 42 dy=dy/sqrt(dy(1)^2+dy(2)^2+dy(3)^2); 43 dz=[v(1,2)*v(2,3)-v(2,2)*v(1,3),v(1,1)*v(2,3)-v(2,1)*v(1,3),v(1,1)*v(2,2)-v(2,1)*v(1,2)]; 44 dz=dz/sqrt(dz(1)^2+dz(2)^2+dz(3)^2); 45![](http://www.cppblog.com/Images/OutliningIndicators/None.gif) 46 T=[dx;dy;dz]; 47 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 48 P=ones(3,11); 49 delta=pi*0.5/10; 50 for i=0:10 51 m=[R*cos(i*delta) R*sin(i*delta) 0]; 52 P(:,i+1)=T'*m'+p'; 53 end 54 plot3(P(1,:),P(2,:),P(3,:),'r+-',A(:,1),A(:,2),A(:,3),'b+-'); 55 grid on; 56 axis on; 57 xlabel('X');ylabel('Y');zlabel('z');
附上以上程序在matlab6.5环境下的运行结果: ![](http://www.cppblog.com/images/cppblog_com/phoenix8848cn/figure2.JPG)
|