for(j=0;j<3;j++)
for(i=0;i<4;i++)
sumXYZ[j]+=X[i][j];
for(i=0;i<2;i++)
X0[i]=sumXYZ[i]/4;//X0,Y0初始化
X0[i]=1/a*f+sumXYZ[2]/4.0;//Z0初始化
do{
R[0][0]=cos(X0[3])*cos(X0[5])-sin(X0[3])*sin(X0[4])*sin(X0[5]);
R[0][1]=-cos(X0[3])*sin(X0[5])-sin(X0[3])*sin(X0[4])*cos(X0[5]);
R[0][2]=-sin(X0[3])*cos(X0[4]);
R[1][0]=cos(X0[4])*sin(X0[5]);
R[1][1]=cos(X0[4])*cos(X0[5]);
R[1][2]=-sin(X0[4]);
R[2][0]=sin(X0[3])*cos(X0[5])+cos(X0[3])*sin(X0[4])*sin(X0[5]);
R[2][1]=-sin(X0[3])*sin(X0[5])+cos(X0[3])*sin(X0[4])*cos(X0[5]);
R[2][2]=cos(X0[3])*cos(X0[4]);
//第一个像点的估计值,其坐标位于X[0][0],X[0][1],X[0][2]
dayue_x[0]=-f*(R[0][0]*(X[0][0]-X0[0])+R[1][0]*(X[0][1]-X0[1])+R[2][0]*(X[0][2]-X0[2]))/(R[0][2]*(X[0][0]-X0[0])+R[1][2]*(X[0][1]-X0[1])+R[2][2]*(X[0][2]-X0[2]));
dayue_x[1]=-f*(R[0][1]*(X[0][0]-X0[0])+R[1][1]*(X[0][1]-X0[1])+R[2][1]*(X[0][2]-X0[2]))/(R[0][2]*(X[0][0]-X0[0])+R[1][2]*(X[0][1]-X0[1])+R[2][2]*(X[0][2]-X0[2]));
//第二个像点的估计值,其坐标位于X[1][0],X[1][1],X[1][2]
dayue_x[2]=-f*(R[0][0]*(X[1][0]-X0[0])+R[1][0]*(X[1][1]-X0[1])+R[2][0]*(X[1][2]-X0[2]))/(R[0][2]*(X[1][0]-X0[0])+R[1][2]*(X[1][1]-X0[1])+R[2][2]*(X[1][2]-X0[2]));
dayue_x[3]=-f*(R[0][1]*(X[1][0]-X0[0])+R[1][1]*(X[1][1]-X0[1])+R[2][1]*(X[1][2]-X0[2]))/(R[0][2]*(X[1][0]-X0[0])+R[1][2]*(X[1][1]-X0[1])+R[2][2]*(X[1][2]-X0[2]));
//第三个像点的估计值,其坐标位于X[2][0],X[2][1],X[2][2]
dayue_x[4]=-f*(R[0][0]*(X[2][0]-X0[0])+R[1][0]*(X[2][1]-X0[1])+R[2][0]*(X[2][2]-X0[2]))/(R[0][2]*(X[2][0]-X0[0])+R[1][2]*(X[2][1]-X0[1])+R[2][2]*(X[2][2]-X0[2]));
dayue_x[5]=-f*(R[0][1]*(X[2][0]-X0[0])+R[1][1]*(X[2][1]-X0[1])+R[2][1]*(X[2][2]-X0[2]))/(R[0][2]*(X[2][0]-X0[0])+R[1][2]*(X[2][1]-X0[1])+R[2][2]*(X[2][2]-X0[2]));
//第四个像点的估计值,其坐标位于X[3][0],X[3][1],X[3][2]
dayue_x[6]=-f*(R[0][0]*(X[3][0]-X0[0])+R[1][0]*(X[3][1]-X0[1])+R[2][0]*(X[3][2]-X0[2]))/(R[0][2]*(X[3][0]-X0[0])+R[1][2]*(X[3][1]-X0[1])+R[2][2]*(X[3][2]-X0[2]));
dayue_x[7]=-f*(R[0][1]*(X[3][0]-X0[0])+R[1][1]*(X[3][1]-X0[1])+R[2][1]*(X[3][2]-X0[2]))/(R[0][2]*(X[3][0]-X0[0])+R[1][2]*(X[3][1]-X0[1])+R[2][2]*(X[3][2]-X0[2]));
for(i=0;i<4;i++)
{
//第i个像点估计值放在dayue_x[2*(i-1)]
A[2*i][0]=(R[0][0]*f+R[0][2]*dayue_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));
A[2*i][1]=(R[1][0]*f+R[1][2]*dayue_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));
A[2*i][2]=(R[2][0]*f+R[2][2]*dayue_x[2*i])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));
A[2*i+1][0]=(R[0][1]*f+R[0][2]*dayue_x[2*i+1])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));
A[2*i+1][1]=(R[1][1]*f+R[1][2]*dayue_x[2*i+1])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));
A[2*i+1][2]=(R[2][1]*f+R[2][2]*dayue_x[2*i+1])/(R[0][2]*(X[i][0]-X0[0])+R[1][2]*(X[i][1]-X0[1])+R[2][2]*(X[i][2]-X0[2]));
A[2*i][3]=dayue_x[2*i+1]*sin(X0[4])-(dayue_x[2*i]/f*(dayue_x[2*i]*cos(X0[5])-dayue_x[2*i+1]*sin(X0[5]))+f*cos(X0[5]))*cos(X0[4]);
A[2*i][4]=-f*sin(X0[5])-dayue_x[2*i]/f*(dayue_x[2*i]*sin(X0[5])+dayue_x[2*i+1]*cos(X0[5]));
A[2*i][5]=dayue_x[2*i+1];
A[2*i+1][3]=-1*dayue_x[2*i]*sin(X0[4])-(dayue_x[2*i+1]/f*(dayue_x[2*i]*cos(X0[5])-dayue_x[2*i+1]*sin(X0[5]))-f*sin(X0[5]))*cos(X0[4]);
A[2*i+1][4]=-1*f*cos(X0[5])-dayue_x[2*i+1]/f*(dayue_x[2*i]*sin(X0[5])+dayue_x[2*i+1]*cos(X0[5]));
A[2*i+1][5]=-dayue_x[2*i];
}
//初始化常数项
for(i=0;i<4;i++)
{
L[2*i]=x[i][0]-dayue_x[2*i];
L[2*i+1]=x[i][1]-dayue_x[2*i+1];
}
//A的转置矩阵
for(i=0;i<8;i++)
for(j=0;j<6;j++)
{
AT[j][i]=A[i][j];
}
//实现A与AT相乘
intk=0;
for(i=0;i<6;i++)
for(j=0;j<6;j++)
Asum[i][j]=0;
for(i=0;i<6;i++)
for(k=0;k<6;k++)
for(j=0;j<8;j++)
Asum[i][k]+=AT[i][j]*A[j][k];
//得到AT*A的逆矩阵存放在inverseAsum[6][6]中
inverse(Asum);
//实现矩阵Asum[6][6]与AT[6][8]的相乘,结果存放在result1[6][8]中
for(i=0;i<6;i++)
for(j=0;j<8;j++)
jieguo1[i][j]=0;
for(i=0;i<6;i++)
for(k=0;k<8;k++)
for(j=0;j<6;j++)
jieguo1[i][k]+=Asum[i][j]*AT[j][k];
//实现result1[6][8]与l[8]的相乘,得到结果放在result2[6]中;
for(i=0;i<6;i++)
jieguo2[i]=0;
for(i=0;i<6;i++)
for(j=0;j<8;j++)
jieguo2[i]+=jieguo1[i][j]*L[j];
for(i=0;i<6;i++)
{
X0[i]=X0[i]+jieguo2[i];
}
ofstreamf7("d:
\\A.txt");
f7<:
fixed;
cout<<"进行第"<\n";
for(i=0;i<6;i++)
{
cout<f7<}
cout<f7.close();
getchar();
m+=1;
}while(abs(jieguo2[3]*206265.0)>6||abs(jieguo2[4]*206265.0)>6||abs(jieguo2[5]*206265.0)>6);
cout<<"\n满足条件的结果为\n";
cout<ofstreamf7("d:
\\A.txt");
f7<:
fixed;
cout.precision(4);
for(i=0;i<6;i++)
{
cout<f7<}
f7.close();
doubleXG[6][1];
for(i=0;i<6;i++)
XG[i][0]=jieguo2[i];
doubleAXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6];
multi(A,XG,AXG,8,6,1);
for(i=0;i<8;i++)//计算改正数
V[i][0]=AXG[i][0]-L[i];
transpose(V,VT,1,8);
multi(VT,V,VTV,1,8,1);
m0=VTV[0][0]/2;
cout<ofstreamf6("d:
\\what.txt");
cout<<"评定完精度的结果是"<for(i=0;i<6;i++)
{
for(intj=0;j<6;j++)
{
D[i][j]=m0*Asum[i][j];
cout<f6<}
cout<f6<}
cout<<"所得中误差为"<for(i=0;i<6;i++)
cout<f6.close();
getchar();
return0;
}
voidinverse(doublec[n][n])
{
inti,j,h,k;
doublep;
doubleq[n][12];
for(i=0;ifor(j=0;jq[i][j]=c[i][j];
for(i=0;ifor(j=n;j<12;j++)
{
if(i+6==j)
q[i][j]=1;
else
q[i][j]=0;
}
for(h=k=0;kfor(i=k+1;i{
if(q[i][h]==0)
continue;
p=q[k][h]/q[i][h];
for(j=0;j<12;j++)
{
q[i][j]*=p;
q[i][j]-=q[k][j];
}
}
for(h=k=n-1;k>0;k--,h--)
for(i=k-1;i>=0;i--)
{
if(q[i][h]==0)
continue;
p=q[k][h]/q[i