|
|
進程被槍斃了.心里不痛快.把做了一半的matlab仿真程序拿出來供大家拍拍磚吧.寫得比較亂.有不明白的可以回復或是發郵件給我.我的郵件是
phoenix8848@gmail.com
1 clear;
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環境下的運行結果:  
|