📄 水平网平差.cpp
字号:
Li=S12-SSS[i].L;
caATPAi(A,Ain,Pi,Li);
}
Pi=1.0/(mT*mT+1.0e-15);
for(i=0;i<=Tnumber-1;i++)
{
int k1=TTT[i].k1;
int k2=TTT[i].k2;
double T12=ca_ab(k1,k2,A,Ain);
Li=(T12-TTT[i].L)*ROU;
caATPAi(A,Ain,Pi,Li);
}
// 处理已知点
for(i=0;i<=knPnumber-1;i++)
{
ATPA[ij(2*i,2*i)]+=1.0e20;
ATPA[ij(2*i+1,2*i+1)]+=1.0e20;
}
}
/////////////////////////////////////////////
//
//
/////////////////////////////////////////////
void known()
{
}
/////////////////////////////////////////////
//
// 计算参数平差值
//
/////////////////////////////////////////////
void ca_dX(void)
{
int i,j,t;
double xi;
t=2*Pnumber+Lnumber;
if(inverse2(ATPA,t)<0)
{
printf("\n法方程系数阵不满秩，程序中断执行。");
exit(0);
}
for(i=0;i<t;i++)
{
xi=0.0;
for(j=0;j<t;j++)xi+=ATPA[ij(i,j)]*ATPL[j];
dX[i]=xi;
if(i<2*Pnumber)XY[i]+=xi;
}
}
//////////////////////////////////////////////////////////////////////////
// 计算最小二乘残差及单位权中误差
//////////////////////////////////////////////////////////////////////////
double caV(void)
{
double VTPV=0.0;
double A[5];
int Ain[5];
A[4]=-1.0;
double Pi=1.0/(ma*ma);
for(int i=0;i<=Lnumber-1;i++)
{
int k1=dir1[i];
Ain[4]=2*Pnumber+i;
for(int j=dir0[i];j<dir0[i+1];j++)
{
int k2=dir2[j];
double T=ca_ab(k1,k2,A,Ain);
double vj=V[j];
for(int s=0;s<5;s++)
vj+=A[s]*dX[Ain[s]];
V[j]=vj;
VTPV+=vj*vj*Pi;
}
}
for(i=0;i<=Snumber-1;i++)
{
int k1=SSS[i].k1;
int k2=SSS[i].k2;
double S12=ca_cd(k1,k2,A,Ain);
double vi=S12-SSS[i].L;
double m=mS1+mS2*SSS[i].L;
Pi=1.0/(m*m);
VTPV+=Pi*vi*vi;
}
Pi=1.0/(mT*mT);
for(i=0;i<=Tnumber-1;i++)
{
int k1=TTT[i].k1;
int k2=TTT[i].k2;
double T12=compute_T12(k1,k2);
double vi=T12-TTT[i].L;
VTPV+=Pi*vi*vi;
}
int n = Nnumber+Snumber+Tnumber;
int t = 2*(Pnumber-knPnumber)+Lnumber;
m0 = sqrt(VTPV/(n-t));
fprintf(resultfp,"\n\n\n 残差二次型： [pvv]=%lf",VTPV);
fprintf(resultfp,"\n 单位权中误差: m=%lf",m0);
return(VTPV);
}
/////////////////////////////////////////////
// 打印坐标平差值及其精度
void PrintCoord()
{
printf("\n 打印平差计算成果...");
fprintf(resultfp,"\n\n ==== 坐标平差值及其精度 ====\n");
fprintf(resultfp,"\nNo. P X Y"
" dX dY mx my M\n");
for(int i=0;i<=Pnumber-1;i++)
{
double xi=XY[2*i];
double yi=XY[2*i+1];
double dxi=dX[2*i];
double dyi=dX[2*i+1];
fprintf(resultfp,"\n%2d %3s ",i+1,Pname[i]);
fprintf(resultfp,"%14.3lf %12.3lf", xi,yi);
fprintf(resultfp,"%7.3lf%7.3lf ", dxi,dyi);
double mx=sqrt(ATPA[ij(2*i,2*i)])*m0;
double my=sqrt(ATPA[ij(2*i+1,2*i+1)])*m0;
double M=sqrt(mx*mx+my*my);
fprintf(resultfp,"%6.3lf",mx);
fprintf(resultfp," %6.3lf",my);
fprintf(resultfp," %6.3lf",M);
}
}
void PrintResult()
{
double A[5];
int Ain[5];
fprintf(resultfp,"\n\n\n ==== 最 后 成 果 表 ====\n");
fprintf(resultfp,"\n\n P1 P2 L V"
" T mT s ms");
for(int i=0;i<=Lnumber-1;i++)
{
int k1=dir1[i];
for(int j=dir0[i];j<dir0[i+1];j++)
{
int k2=dir2[j];
if(j==dir0[i])
{
fprintf(resultfp,"\n\n%6s",Pname[k1]);
}
else fprintf(resultfp,"\n ");
fprintf(resultfp,"%7s",Pname[k2]);
fprintf(resultfp,"%15s%7.2lf",dms(L[j],2),V[j]);
double T = ca_ab(k1,k2,A,Ain);
double mT = sqrt(qq(A,Ain))*m0;
fprintf(resultfp,"%15s %5.2lf",dms(T,2),mT);
double S = ca_cd(k1,k2,A,Ain);
double mS=sqrt(qq(A,Ain))*m0;
fprintf(resultfp,"%10.3lf %6.3lf",S,mS);
}
}
}
/////////////////////////////////////////////
//
// 平面方位角计算
//
/////////////////////////////////////////////
double compute_T12(int k1,int k2) // 返回值为弧度值,
{
if(XY[2*k1]>1.0e29 || XY[2*k2]>1.0e29) return 1.0e30;
double dx=XY[2*k2]-XY[2*k1];
double dy=XY[2*k2+1]-XY[2*k1+1];
double T=atan2(dy,dx);
if(T<0.0)T=T+2.0*PAI;
return T;
}
//////////////////////////////////////////////////////////////////////////
// 权 倒 数 计 算
double qq(double A[],int Ain[])
{
int i,j,k1,k2;
double q=0.0;
for(i=0;i<4;i++)
{
k1=Ain[i];
for(j=0;j<4;j++)
{
k2=Ain[j];
q+=ATPA[ij(k1,k2)]*A[i]*A[j];
}
}
return q;
}
/////////////////////////////////////////////
//
// 拷贝字符串
//
/////////////////////////////////////////////
void stpcpy(char ch1[],char ch2[])
{
int i=0;
while(1)
{
ch1[i]=ch2[i];
if(ch1[i]=='\0')return;
i++;
}
}
//////////////////////////////////////////////
//
// 度分秒字符串
//
/////////////////////////////////////////////
#include <strstrea.h>
#include <iomanip.h>
char *dms(double a,int n)
{
static int num=-1;
static char ch[10][90];
num++;
if(num==10)num=0;
int d,m,sign=1;
a=a*2.06264806247096363e+05;
if(a<0.0){ sign=-1; a=-a; }
d=(int)((a+0.00000001)/3600.0);
a=a-d*3600.0;
m=(int)((a+0.00000001)/60.0);
a=a-m*60.0;
if(a<0.000000001)a=0.0;
ostrstream out(ch[num],90);
out<<setw(4);
if(d==0)
{
if(sign==1)out<<"0";
else out<<"-0";
}
else{ d=d*sign; out<<d; }
out<<" "<<setw(2)<<setfill('0')<<m;
out.setf(ios::showpoint);
out.setf(ios::fixed);
out<<" ";
out<<setw(3+n)<<setprecision(n)<<a<<ends;
return ch[num];
}
/////////////////////////////////////////////
//
// 对 称 正 定 矩 阵 求 逆
// (a[]仅存下三角元素）
/////////////////////////////////////////////
int inverse2(double a[],int n)
{
int i,j,k,m;
double w,g,*b;
b=new double[n];
for(k=0;k<=n-1;k++)
{
w=a[0];
if(w+1.0==1.0)
{
delete []b; return (-2);
}
m=n-k-1;
for(i=1;i<=n-1;i++)
{
g=a[i*(i+1)/2];
b[i]=g/w;
if(i<=m)b[i]=-b[i];
for(j=1;j<=i;j++)
a[(i-1)*i/2+j-1]=a[i*(i+1)/2+j]+g*b[j];
}
a[n*(n+1)/2-1]=1.0/w;
for(i=1;i<=n-1;i++)
a[(n-1)*n/2+i-1]=b[i];
}
delete []b;
return 1;
}
////////////////////////////////////////////////////
//
// 主函数
//
////////////////////////////////////////////////////
void main()
{
puts("\n\n\n\n 欢迎使用水平网平差软件");
while(1)
{
printf("\n默认的数据文件和结果文件: %s %s",DATAFILE,RESULTFILE);
printf("\n开始计算？（y/n)");
char yn=getch();
if(yn=='y' || yn=='\r')break;
printf("\n请从键盘输入数据文件名:");
scanf("%s",DATAFILE);
printf("\n请从键盘输入结果文件名:");
scanf("%s",RESULTFILE);
}
if((resultfp=fopen(RESULTFILE,"w"))==NULL)
{
printf("结果文件打不开");
getch();
exit(0);
}
puts("\n 输入数据文件...");
inputdata();
puts("\n\n 打印原始数据...");
printdata();
puts("\n\n 近似坐标计算...");
//caxy0();
ca_x0y0();
puts("\n 最小二乘平差....");
caATPA();
ca_dX();
VTPV=caV();
PrintCoord();
PrintResult();
fprintf(resultfp,"\n\n\n");
Printwcty();
fcloseall();
printf("\n\n 计算全部结束. 结果存入:%s\n\n",RESULTFILE);
}
//////////////////////////////////////////////////////////////////////////
// 导线网近似坐标计算函数
int ca_x0y0()
{
//设置未知点标志
for(int i=knPnumber;i<Pnumber;i++)
{
XY[2*i]=1.0e30;
}
int unknow=Pnumber-knPnumber; //未知点数
int No=0;
while(1)
{
if(unknow==0)return 1;
No++;
if(No>(Pnumber-knPnumber)) return 0;
for( i=0;i<=Lnumber-1;i++)
{
int k1=dir1[i]; //测站点号
double x1=XY[2*k1];
double y1=XY[2*k1+1];
if(x1>1.0e29) continue; //测站点是未知点
int j0=dir0[i];
for(int j=j0; j<dir0[i+1];j++)
{
int k3=dir2[j];
double T13=compute_T12(k1,k3);
if(T13>1.0e29)T13=get_T12(k1,k3);
if(T13>1.0e29) continue;
for(int k=j0;k<dir0[i+1];k++)
{
int k2=dir2[k];
double x2=XY[2*k2];
if(x2<1.0e29)continue;
double T12=T13+L[k]-L[j];
double S12=get_S12(k1,k2);
if(S12>1.0e29) continue;
XY[2*k2]=x1+S12*cos(T12);
XY[2*k2+1]=y1+S12*sin(T12);
unknow--;
}
}
}
}
return 1;
}
//在方位角观测值中查找
double get_T12(int k1,int k2)
{
for(int i=0; i<Tnumber; i++)
{
if(k1==TTT[i].k1 && k2==TTT[i].k2)
{
return TTT[i].L;
}
if(k1==TTT[i].k2 && k2==TTT[i].k1)
{
return TTT[i].L+PAI;
}
}
return 1.0e30;
}
//////////////////////////////////////////////////////////////////////////
//在边长观测值中查找
double get_S12(int k1,int k2)
{
for(int i=0; i<Snumber; i++)
{
if(k1==SSS[i].k1 && k2==SSS[i].k2)
{
return SSS[i].L;
}
if(k1==SSS[i].k2 && k2==SSS[i].k1)
{
return SSS[i].L;
}
}
return 1.0e30;
}
//////////////////////////////////////////////////////////////////////////
//误差椭圆
void Printwcty()
{
printf("\n 误差椭圆...");
fprintf(resultfp,"\n\n ==== 误差椭圆====\n");
fprintf(resultfp,"\nNo. name E F a1 a2");
double K,t1,t2,E,F,a1,a2;
for(int i=0;i<=Pnumber-1;i++)
{
fprintf(resultfp,"\n%2d %3s ",i+1,Pname[i]);
QX=ATPA[ij(2*i,2*i)];
QY=ATPA[ij(2*i+1,2*i+1)];
QXY=ATPA[ij(2*i,2*i+1)];
K=sqrt((QX-QY)*(QX-QY)+4*QXY*QXY);
t1=(QX+QY+K)/2;
t2=(QX+QY-K)/2;
E=m0*sqrt(t1);
F=m0*sqrt(t2);
a1=atan((t1-QX)/QXY);
a2=atan((t2-QX)/QXY);
fprintf(resultfp,"%10.3lf",E);
fprintf(resultfp," %10.3lf",F);
fprintf(resultfp," %10.3lf",a1);
fprintf(resultfp,"%10.3lf",a2);
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -