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