📄 adjpro050404.cpp
字号:
for(i=0;i<a.obnum+a.fixafn+a.fixdisn;i++)
if(i!=obi)a.L[i].A0=-PI;
zheng(a.L[obi]); // 坐标正算
// 计算改正后近似方位角
do{flag=0;
statangc(a);
for(int i=0;i<a.obnum;i++)
if(a.L[i].A0==-PI){flag=1;break;}
}while(flag==1);
// 计算改正后各点的近似坐标
do{flag=0;
for(int i=0;i<a.obnum-1;i++)
for(int j=i+1;j<a.obnum;j++)
XY0ang(a.L[i],a.L[j]);
for(i=0;i<a.allpnum;i++)
if(a.Pt[i].X0==-PI)
{ flag=1;break; }
}while(flag==1);
// for(i=0;i<a.allpnum;i++)
// cout<<a.Pt[i].name<<" "<<a.Pt[i].X0<<" "<<a.Pt[i].Y0<<" "
// <<a.Pt[i].X<<" "<<a.Pt[i].Y<<" "<<endl;
return 1;
}
//************************计算控制网未知点的近似坐标**********************************************
int setx0y0(XYnet &a)
{
int n=a.obnum+a.fixdisn+a.fixafn;
// 1.计算近似坐标、近似边长确定的方位角与边长
int t(0);
do{
for(int i=0;i<n;i++)
{
if(a.L[i].startp->X0!=-PI && a.L[i].endp->X0!=-PI && a.L[i].A0==-PI)
{
//1.1 近似坐标确定的边的方位角
a.L[i].A0=h_d(afa(*(a.L[i].startp),*(a.L[i].endp)));
for(int k=0;k<n;k++)
{
if(a.L[i].startp==a.L[k].endp && a.L[i].endp==a.L[k].startp && a.L[k].A0==-PI)
{
if(d_h(a.L[i].A0)-PI>=0)
a.L[k].A0=h_d(d_h(a.L[i].A0)-PI);
else a.L[k].A0=h_d(d_h(a.L[i].A0)+PI);
}
if(a.L[i].startp==a.L[k].startp && a.L[i].endp==a.L[k].endp && a.L[k].A0==-PI)
a.L[k].A0=a.L[i].A0;
}
}
//1.2 边长观测值计算的近似边长
if(a.L[i].style==2 || a.L[i].style==4)
for(int k=0;k<n;k++)
if(a.L[i].startp==a.L[k].endp && a.L[i].endp==a.L[k].startp && a.L[k].dist0==-PI
|| a.L[i].startp==a.L[k].startp && a.L[i].endp==a.L[k].endp && a.L[k].dist0==-PI)
a.L[k].dist0=a.L[i].dist0;
//1.3 由已知近似方位角推算未知近似方位角
if(a.L[i].A0!=-PI)
for(int k=0;k<n;k++)
{
if(a.L[i].startp==a.L[k].endp && a.L[i].endp==a.L[k].startp && a.L[k].A0==-PI)
{
if(d_h(a.L[i].A0)-PI>=0)
a.L[k].A0=h_d(d_h(a.L[i].A0)-PI);
else a.L[k].A0=h_d(d_h(a.L[i].A0)+PI);
}
if(a.L[i].startp==a.L[k].startp && a.L[i].endp==a.L[k].endp && a.L[k].A0==-PI )
a.L[k].A0=a.L[i].A0;
}
}
// 2. 以测站计算观测值近似方位角
statangc(a);
// 3.近似坐标计算
// 3.1 两方位角交会法
for(i=0;i<n-1;i++)
for(int j=i+1;j<n;j++)
XY0ang(a.L[i],a.L[j]);
// 3.2 坐标正算法
for(i=0;i<a.obnum+a.fixafn+a.fixdisn;i++)
zheng(a.L[i]);
// 3.3 角度后方交会法
for(i=0;i<a.obnum-2;i++)
for(int j=i+1;j<a.obnum-1;j++)
for(int k=j+1;k<a.obnum;k++)
houj(a.L[i],a.L[j],a.L[k]);
// 3.4边长后方交会法
for(i=0;i<a.obnum+a.fixafn+a.fixdisn-2;i++)
for(int j=i+1;j<a.obnum+a.fixafn+a.fixdisn-1;j++)
for(int k=j+1;k<a.obnum+a.fixafn+a.fixdisn;k++)
XY0dist(a.L[i],a.L[j],a.L[k]);
// 3.5 无定向导线
Udxdsetx0y0(a);
// 4. 迭代条件的判断
t=0;
for( i=0;i<a.allpnum;i++)
{
// cout<<a.Pt[i].name<<" "<<a.Pt[i].X0<<" "<<a.Pt[i].Y0<<endl;
if(a.Pt[i].X0==-PI || a.Pt[i].Y0==-PI)
t=1;
}
}while(t==1);
cout<<"---------------------------------------------------------------------------"<<endl;
cout<<" 未知点近似坐标 "<<endl;
cout<<"---------------------------------------------------------------------------"<<endl;
for(int ii=a.fixpnum; ii<a.allpnum;ii++)
cout<<" "<<a.Pt[ii].name<<" "<<setf(a.Pt[ii].X0,3)<<" "<<setf(a.Pt[ii].Y0,3)<<endl;
cout<<"---------------------------------------------------------------------------"<<endl;
return 1;
}
//********************************设置平面网平差的A,P,L*******************************************
int setXYadj(XYnet &a)
{
// 1 确定平差网基本信息
a.aa.m=a.obnum+a.fixafn+a.fixdisn;
a.aa.n=2*a.allpnum+a.statnum;
// 2 确定观测值的权镇P
for(int i=0;i<a.aa.m;i++)
{
for(int j=0;j<a.aa.m;j++)
a.aa.P[i][j]=0;
if(a.L[i].style==1) // 设置角度观测值的权
{
a.aa.P[i][i]=1;
}
if(a.L[i].style==2) // 设置距离观测值的权
{
double mlml=a.msa+a.msb/1000000*a.L[i].dist;
mlml*=mlml;
a.aa.P[i][i]=a.mangle/rou*a.mangle/rou/mlml;
}
if(a.L[i].style==3 || a.L[i].style==4)// 设置固定距离与方位角观测值的权
a.aa.P[i][i]=1000000;
}
// matdis(a.aa.P,a.aa.m,a.aa.m);
// 3 确定误差方程系数阵A
for(i=0;i<a.aa.m;i++)
{
double s=dist(*(a.L[i].startp),*(a.L[i].endp));
// cout<<s<<endl;
double ss=s*s;
double dy=a.L[i].endp->Y0-a.L[i].startp->Y0;
double dx=a.L[i].endp->X0-a.L[i].startp->X0;
double aki=dy/ss; double bki=-1*dx/ss;
for(int j=0;j<a.aa.n;j++)
{a.aa.A[i][j]=0;
if(a.L[i].style==1) // 方向观测值误差方程系数
{
a.aa.A[i][2*a.allpnum+a.L[i].sti]=-1; // 测站定向角未知数系数
a.aa.A[i][2*(a.L[i].startp->i)]=aki;
a.aa.A[i][2*(a.L[i].startp->i)+1]=bki;
a.aa.A[i][2*(a.L[i].endp->i)]=-aki;
a.aa.A[i][2*(a.L[i].endp->i)+1]=-bki;
}
if(a.L[i].style==2 || a.L[i].style==4) // 边长观测值、固定边长误差方程系数
{
a.aa.A[i][2*(a.L[i].startp->i)]=bki*s;
a.aa.A[i][2*(a.L[i].startp->i)+1]=-aki*s;
a.aa.A[i][2*(a.L[i].endp->i)]=-bki*s;
a.aa.A[i][2*(a.L[i].endp->i)+1]=aki*s;
}
if(a.L[i].style==3) // 固定方位角误差方程系数
{
a.aa.A[i][2*a.L[i].startp->i]=aki;
a.aa.A[i][2*a.L[i].startp->i+1]=bki;
a.aa.A[i][2*a.L[i].endp->i]=-aki;
a.aa.A[i][2*a.L[i].endp->i+1]=-bki;
}
}
}
// matdis(a.aa.A,a.aa.m,a.aa.n);
// 4 确定误差方程常数项L
for(i=0;i<a.aa.m;i++)
{
a.L[i].A0=h_d(afa(*(a.L[i].startp),*(a.L[i].endp)));
a.L[i].dist0=dist(*(a.L[i].startp),*(a.L[i].endp));
if(a.L[i].style==3)
{
a.aa.l[i][0]=d_h(a.L[i].A-a.L[i].A0);
}
if(a.L[i].style==2 || a.L[i].style==4)
{
a.aa.l[i][0]=a.L[i].dist-a.L[i].dist0;
}
}
//方向观测值误差方程常数项
int n=0;
for(int j=0;j<a.statnum;j++)
{
for(int i=n;i<n+a.st[j].aglnum;i++)
{
double scd=d_h(a.L[i].A0)-d_h(a.L[n].A0);
if(scd<0) scd+=2*PI;
a.aa.l[i][0]=d_h(a.L[i].angle)-scd;
}
n+=a.st[j].aglnum+a.st[j].disnum;
}
//matdis(a.aa.l,a.aa.m);
return 1;
}
//********************************控制网平差计算**************************************************
int doXYadj(XYnet &a)
{
// 6 平差计算
int r0=a.fixafn+a.fixdisn;
if(!doadj(a.aa,2*a.fixpnum,r0))return 0; // 其中第二项为控制点数与测站数之和
cout<<"---------------------------------------------------------------------------"<<endl;
cout<<" 未知数改正数 "<<endl;
cout<<"---------------------------------------------------------------------------"<<endl;
matdis(a.aa.X,a.aa.n);
cout<<"---------------------------------------------------------------------------"<<endl;
cout<<" 验后单位权中误差："<<a.aa.m0*206264.8<<" 秒 "<<endl;
cout<<"---------------------------------------------------------------------------"<<endl;
cout<<" 平差后未知点坐标: "<<endl<<endl;
cout<<"---------------------------------------------------------------------------"<<endl;
for(int i=a.fixpnum;i<a.allpnum;i++)
{
a.Pt[i].X=a.Pt[i].X0+a.aa.X[2*i][0];
a.Pt[i].Y=a.Pt[i].Y0+a.aa.X[2*i+1][0];
a.Pt[i].X0+=a.aa.X[2*i][0];
a.Pt[i].Y0+=a.aa.X[2*i+1][0];
cout.precision(11);cout.width(10);
cout<<" "<<a.Pt[i].name<<" "<<setf(a.Pt[i].X,4)<<" "<<setf(a.Pt[i].Y,4)<<endl;
}
cout<<"---------------------------------------------------------------------------"<<endl;
cout<<" 未知点坐标的协方差阵:"<<endl<<endl;
cout<<"---------------------------------------------------------------------------"<<endl;
for(i=2*a.fixpnum;i<2*a.allpnum;i++)
{ cout<<" ";
for(int j=2*a.fixpnum;j<2*a.allpnum;j++)
{
cout<<a.aa.m0*a.aa.m0*a.aa.QXX[i][j]<<" ";
}cout<<endl;
}
cout<<"---------------------------------------------------------------------------"<<endl;
cout<<" 观测值改正数:"<<endl<<endl;
cout<<"---------------------------------------------------------------------------"<<endl;
matdis(a.aa.V,a.aa.m);
cout<<"---------------------------------------------------------------------------"<<endl;
//7 点位中误差与误差椭圆元素计算
for(i=a.fixpnum;i<a.allpnum;i++)
{
// 坐标中误差计算
a.Pt[i].mX=a.aa.m0*sqrt(a.aa.QXX[2*i][2*i]);
a.Pt[i].mY=a.aa.m0*sqrt(a.aa.QXX[2*i+1][2*i+1]);
// 误差椭圆参数计算
double K=sqrt((a.aa.QXX[2*i][2*i]-a.aa.QXX[2*i+1][2*i+1])
*(a.aa.QXX[2*i][2*i]-a.aa.QXX[2*i+1][2*i+1])
+4*a.aa.QXX[2*i][2*i+1]*a.aa.QXX[2*i][2*i+1]);
a.Pt[i].mp=sqrt(a.Pt[i].mX*a.Pt[i].mX+a.Pt[i].mY*a.Pt[i].mY);
a.Pt[i].E=sqrt(1/2.0*a.aa.m0*a.aa.m0*(a.aa.QXX[2*i][2*i]+a.aa.QXX[2*i+1][2*i+1]+K));
a.Pt[i].F=sqrt(1/2.0*a.aa.m0*a.aa.m0*(a.aa.QXX[2*i][2*i]+a.aa.QXX[2*i+1][2*i+1]-K));
double F2=atan2(2*a.aa.QXX[2*i][2*i+1],a.aa.QXX[2*i][2*i]-a.aa.QXX[2*i+1][2*i+1]);
if(F2<=0) F2+=2*PI;F2/=2.0;
if(a.aa.QXX[2*i][2*i+1]>0 && F2>PI/2.0 && F2<PI
|| a.aa.QXX[2*i][2*i+1]>0 && F2>PI*3/2.0 && F2<2*PI)
F2+=PI/2.0;
if(a.aa.QXX[2*i][2*i+1]<0 && F2<PI/2.0 && F2>0
|| a.aa.QXX[2*i][2*i+1]<0 && F2>PI && F2<PI*3/2.0)
F2+=PI/2.0;
if(F2>=2*PI) F2-=2*PI;
a.Pt[i].T=h_d(F2);
}
return 1;
}
//*******************************平面网的文件输入*************************************************
int finXYnet(XYnet &a,char *fname) // 文件输入平面网函数
{
ifstream in(fname,ios::nocreate); // 建立文件流,并与输入文件名建立关联
if(!in) {cout<<fname<<" error: file does not exist! "<<endl; return 0;}
// 文件现实性判断
// 1 输入网的基本信息
char name[20];
in>>a.netname;
in>>a.allpnum;
in>>a.fixpnum;
in>>a.statnum;
// in>>a.obnum;
a.obnum=0;
in>>a.fixdisn;
in>>a.fixafn;
in>>a.mangle;
in>>a.msa;
in>>a.msb;
int n(a.fixpnum); // n为已输入名字的点的个数
// 2输入控制点信息
for(int i=0;i<a.fixpnum;i++)
{
in>>a.Pt[i].name>>a.Pt[i].X>>a.Pt[i].Y;
a.Pt[i].fixed=1; // 控制点标记
a.Pt[i].X0=a.Pt[i].X;
a.Pt[i].Y0=a.Pt[i].Y;
a.Pt[i].mX=a.Pt[i].mY=0;
a.Pt[i].i=i; // 控制点编号，从0到a.fixpnum-1
}
// 3 输入未知点相关信息（名字在后面输入）
for(i=a.fixpnum;i<a.allpnum;i++)
{
a.Pt[i].fixed=0; // 未知点标记
a.Pt[i].X=a.Pt[i].Y=-PI;
a.Pt[i].X0=a.Pt[i].Y0=-PI;
a.Pt[i].mX=a.Pt[i].mY=0;
*(a.Pt[i].name)=0;
a.Pt[i].i=i; // 为未知点编号,从a.fixpnum到a.allpnum-1
}
a.anglenum=0;
a.distnum=0;
// 4 输入测站及观测值
int obsnum(0);
for(i=0;i<a.statnum;i++)
{
// 4-1 输入测站信息
int t=0; // 点名比较标志
in>>name; // 输入测站名
for(int k=0;k<n;k++)
if(strnicmp(name,a.Pt[k].name,20)==0)
{
a.st[i].stp=&(a.Pt[k]); // 找到同名点，起点指针指向该点
t++; // 找到标志
}
if(t==0) {strcpy(a.Pt[n].name,name);
a.st[i].stp=&(a.Pt[n]); // 找不到同名点，该名输给新点
n++;}
in>>a.st[i].obnum;
a.st[i].aglnum=0; // 测站角度观测个数
a.st[i].disnum=0; // 测站距离观测个数
a.st[i].i=i;
// 4-2 输入本站观测值
for(int j=obsnum;j<obsnum+a.st[i].obnum;j++)
{
in>>name; t=0; // 输入照准点名，操作过程同上
a.L[j].startp=a.st[i].stp;
for(k=0;k<n;k++)
if(strnicmp(name,a.Pt[k].name,20)==0)
{
a.L[j].endp=&(a.Pt[k]);
t++;
}
if(t==0)
{strcpy(a.Pt[n].name,name);
a.L[j].endp=&(a.Pt[n]);
n++;
}
in>>a.L[j].style;
a.L[j].A=a.L[j].A0=-PI;
if(a.L[j].style==1)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -