我做了一个程序,最后要求误差的分布,整个过程只有实数,没有复数产生,但运行 M(r)=acos(abs(dot(N(r,:),NL))/(norm(N(r,:))*norm(NL)))*180/(pi);有时会产生复数,头疼呀,查了很多遍了,但找不出问题,望各位高手出手相助,不胜感激,下面是程序:
%用蒙特卡洛方法测定调平装置的误差分布
%以平面2x+3y+6z-8=0为例,测定平面上点(c1,c2,c3)=(1,2,0)处的法向误差分布
N=[];%计算法向量
M=[];%误差值
angle=20*pi/180;%传感头发射的激光与主轴的夹角
NL=[2,3,6];%平面2x+3y+6z+8=0的法向量
point=0;%计数点
s=unifrnd(-30*cos(angle),0,1,200);
%产生服从均匀分布U(-30cos20,0)的随机数200个
A1=[1/tan(angle),0,1;
2,3,6;
0,1,0];
A2=[-2*3^(1/2),-2,0;
0,4/tan(angle),2*3^(1/2);
2,3,6];
A3=[2*3^(1/2),-2,0;
0,4/tan(angle),-2*3^(1/2);
2,3,6];
%传感器1,2,3发射激光的光线方程与平面方程的联立的系数矩阵
for t=1:200
b1=[s(t)+1/tan(angle);8;2];
b2=[-2*3^(1/2)-4;8/tan(angle)+2*3^(1/2)*s(t);8];
b3=[2*3^(1/2)-4;8/tan(angle)-2*3^(1/2)*s(t);8];
x1=A1\b1;
x2=A2\b2;
x3=A3\b3;
%三束激光线与平面的焦点坐标
if (norm(x1'-[1,2,s(t)])<=30)&(norm(x2'-[1,2,s(t)])<=30)&(norm(x3'-[1,2,s(t)])<=30)
%判断是否超出测量范围
point=point+1;
h1=30-norm(x1'-[1,2,s(t)]);
h2=30-norm(x2'-[1,2,s(t)]);
h3=30-norm(x3'-[1,2,s(t)]);
L2=65;
T=95*sin(angle);
N(point,:)=[1/2*3^(1/2)*(2*T-(2*L2+h2+h3)*sin(angle))*(h2-h1)*cos(angle)-1/2*(h2-h3)*cos(angle)*3^(1/2)*(T-(L2+h2)*sin(angle)),(h2-h3)*cos(angle)*(-3/2*T+(3/2*L2+h1+1/2*h2)*sin(angle))-(1/2*h2-1/2*h3)*sin(angle)*(h2-h1)*cos(angle), 1/2*(1/2*h2-1/2*h3)*sin(angle)*3^(1/2)*(T-(L2+h2)*sin(angle))-1/2*3^(1/2)*(2*T-(2*L2+h2+h3)*sin(angle))*(-3/2*T+(3/2*L2+h1+1/2*h2)*sin(angle))];
%计算近似向量
if point==100
break;
end
end
end
for r=1:100
M(r)=acos(abs(dot(N(r,:),NL))/(norm(N(r,:))*norm(NL)))*180/(pi);%本次测量的误差
end
[h,p,j,cv]=jbtest(M)%判断是否为正态分布
[n,xout]=hist(M,25)
bar(xout,n)%绘制频率直方图
E=mean(M)%计算样本均值
sigma=std(M)%计算标准差