📄 水平网平差.cpp
字号:
//
// 程 序 设 计 作 业
//
// 题目：水平网平差程序设计、调试及运行。
// 程序功能:(1)从键盘读入数据文件名和结果文件名。
// (2)输入已知坐标和观测方向值。
// (3)输出原始数据。
// (4)近似坐标计算。
// (5)最小二乘平差。
// (6)精度估计。
// 上交成果：源程序、运行结果文件。
//
// 本文件是教员在课堂上讲过的程序，但删除了部分函数的代码，
// 请将空函数的内容补充完整，编译、调试。
// 编译成功后，运行程序，从“导线网data.txt” 文件中读入已
// 知坐标和观测值，结果存入：“result.dat”文件。
//
// 注意：本作业的要求与教材的程序不同，所有点均在同一投影带，
// 没有换带问题，网中方位角和边长当作观测值处理。
//
//
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <string.h>
//自定义方位角、边长数据结构
struct ST{
int k1;
int k2;
double L;
};
///////////////////////////////////////////////////
// 全局变量
char netname[100]; //网名
int Pnumber; //总点数
int knPnumber;//已知点总数
int Lnumber; //方向组总数
int Nnumber; //方向值总数
int Snumber; //边长总数
int Tnumber; //方位角总数
double ma; //方向值中误差
double mS1; //边长固定误差
double mS2; //边长比例误差
double mT; //方位角中误差
double QX;
double QY;
double QXY;
char **Pname; //点名地址数组
double *XY; //坐标数组
double *L; //方向观测值
int *dir1; //测站点号
int *dir2; //观测方向的点号
int *dir0; //测站零方向在方向值
//数组L中的下标
ST *SSS; //边长观测值数组
ST *TTT; //方位角观测值数组
char DATAFILE[100]={"导线网Data.txt"}; //数据文件名
char RESULTFILE[100]={"result.txt"}; //结果文件名
FILE *resultfp; //输出文件指针
double *V;
double *ATPA,*ATPL;
double *dX;
double VTPV,m0;
double ROU=2.062648062470963630e+05;
double PAI=3.141592653589793120e+00;
////////////////////////////////////////////////////
int inverse2(double a[],int n);
void stpcpy(char ch1[],char ch2[]);
void known();
double compute_T12(int k1,int k2); // 返回值为弧度值,
double get_S12(int k1,int k2); //在边长观测值中查找
double get_T12(int k1,int k2); //方位角观测值中查找
double qq(double A[],int Ain[]);
char *dms(double a,int n);
int ca_x0y0();
char *dms(double a,int n);
void PrintResult();
void PrintCoord();
void Printwcty();
/////////////////////////////////////////////
//
// 编 点 号
//
/////////////////////////////////////////////
int cpoint(char *name)
{
int i=0,j,c;
static int pnum=0;
while(i<pnum)
{
if(strcmp(name,Pname[i])==0)return i;
i++;
}
if(i<Pnumber)
{
c=strlen(name)+1;
Pname[i]=new char[c];
stpcpy(Pname[i],name);
pnum++;
return i;
}
else
{
printf("\n\nError: Total points >%d ",Pnumber);
printf("\nPoint names:\n");
for(j=0;j<Pnumber;j++)printf("%10s",Pname[j]);
printf("%10s",name);
getch();
exit(0);
}
return i;
}
////////////////////////////////////////////////////
//
// 下标计算函数
//
///////////////////////////////////////////////////
int ij(int i,int j)
{
return j<=i ? i*(i+1)/2+j : j*(j+1)/2+i;
}
//////////////////////////////////////////////////////////////////////////
// 输入数据文件
void inputdata()
{
FILE *fp;
if((fp=fopen(DATAFILE,"r"))==NULL)
{
printf("\n 无法打开数据文件\"%s\"!",DATAFILE);
getch();
exit(0);
}
fscanf(fp,"%s", netname); //网名
fscanf(fp,"%d%d%d%d",&Pnumber, // 总点数
&knPnumber, // 已知点数
&Lnumber, // 测站总数
&Nnumber // 观测值总数
);
fscanf(fp,"%d%d",&Snumber, // 边长点数
&Tnumber // 方位角数
);
XY=new double [2*Pnumber];
SSS=new ST[Snumber];
TTT=new ST[Tnumber];
int t=2*Pnumber+Lnumber;
int tt=t*(t+1)/2;
ATPL=new double [t];
ATPA=new double [tt];
dX=new double [t];
Pname=new char* [Pnumber];
dir0=new int [Lnumber+1]; //各测站首方向在L的偏移量
dir1=new int [Lnumber];
int n=Nnumber;
dir2=new int [n];
L=new double [n];
V=new double [n]; //在计算V之前，放的是误差方程自由项l
int unPnumber=Pnumber-knPnumber;
//输入观测精度表
fscanf(fp,"%lf%lf%lf%lf",&ma,&mS1,&mS2,&mT);
char name[100];
for(int i=0;i<=knPnumber-1;i++)//输入已知点坐标
{
fscanf(fp,"%s",name);
int k=cpoint(name);
fscanf(fp,"%lf%lf",&XY[2*k],&XY[2*k+1]);
}
int d,m;
double s;
dir0[0]=0;
int ii=0;
for(i=0;i<=Lnumber-1;i++) //输入方向值，按测站
{
int sn;
fscanf(fp,"%s%d",name,&sn); // sn: 测站方向数
dir0[i+1]=dir0[i]+sn;
int k1=cpoint(name);
if(sn<=0){ printf("\n方向数等于零!"); exit(0); }
dir1[i]=k1;
for(int j=0;j<=sn-1;j++)
{
fscanf(fp,"%s%d%d%lf",name,&d,&m,&s);
dir2[ii]=cpoint(name);
L[ii]=(d*3600.0+m*60.0+s)/ROU;
ii++;
}
}
for(i=0;i<=Snumber-1;i++) //输入边长
{
fscanf(fp,"%s",name);
SSS[i].k1=cpoint(name);
fscanf(fp,"%s",name);
SSS[i].k2=cpoint(name);
fscanf(fp,"%lf",&SSS[i].L);
}
for(i=0;i<=Tnumber-1;i++) //输入方位角
{
fscanf(fp,"%s",name);
TTT[i].k1=cpoint(name);
fscanf(fp,"%s",name);
TTT[i].k2=cpoint(name);
fscanf(fp,"%d%d%lf",&d,&m,&s);
TTT[i].L=(d*3600.0+m*60.0+s)/ROU;
}
fclose(fp);
}
/////////////////////////////////////////////////////
//
// 打印原始数据
//
/////////////////////////////////////////////////////
void printdata(void)
{
fprintf(resultfp,"\n网名：%s " ,netname);
fprintf(resultfp,"\n总点数：%d",Pnumber);
fprintf(resultfp,"\n已知点数：%d",knPnumber);
fprintf(resultfp,"\n方向组总数：%d",Lnumber);
fprintf(resultfp,"\n方向值总数: %d",Nnumber);
fprintf(resultfp,"\n边长总数：%d",Snumber);
fprintf(resultfp,"\n方位角总数：%d",Tnumber);
fprintf(resultfp,"\n\n ==== 已知点坐标 ====\n");
for(int i=0;i<=knPnumber-1;i++)
{
fprintf(resultfp,"\n%4s ",Pname[i]);
fprintf(resultfp,"%12.3lf%14.3lf",
XY[2*i],XY[2*i+1]);
}
fprintf(resultfp,"\n\n ==== 方向观测值 ====\n");
for(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%4s %d \n",
Pname[k1],dir0[i+1]-dir0[i]);
fprintf(resultfp,"%4s %15s\n",
Pname[k2],dms(L[j],2));
}
}
for(i=0;i<=Snumber-1;i++)
{
if(i==0)fprintf(resultfp,
"\n\n ==== 边长观测值 ====\n\n");
int k1=SSS[i].k1;
int k2=SSS[i].k2;
fprintf(resultfp,"%4s %4s %10.3lf\n",
Pname[k1],Pname[k2],SSS[i].L);
}
for(i=0;i<=Tnumber-1;i++)
{
if(i==0)fprintf(resultfp,
"\n\n ==== 方位角观测值 ====\n\n");
int k1=TTT[i].k1;
int k2=TTT[i].k2;
fprintf(resultfp,"%4s %4s %s\n",
Pname[k1],Pname[k2],dms(TTT[i].L,3));
}
}
///////////////////////////////////////////////////
//
// 根据三点号返回以p为顶点的一个三角形内角
//
///////////////////////////////////////////////////
double getangle(int p,int p1,int p2)
{
int i,j,k,i1,i2;
double A;
for(i=0;i<=Lnumber-1;i++)
{
i1=i2=-1;
k=dir1[i];
for(j=dir0[i];j<dir0[i+1];j++)
{
if(k!=p)break;
if(p1==dir2[j])i1=j;
if(p2==dir2[j])i2=j;
}
if(i1>=0&&i2>=0)
{
A=L[i2]-L[i1];
if(A<0.0)A+=2.0*PAI;
return A;
}
}
return -1.0; //找不到符合条件的角，返回负值
}
//////////////////////////////////////////////////////////////////////////
// 计算近似坐标(三角网)
void caxy0()
{
int i,j,jj,k,k1,k2,s,f;
double A,B,C,xA,yA,xB,yB;
double sinA,sinB,cosA,cosB,sinC,xk,yk;
for(i=knPnumber;i<=Pnumber-1;i++)XY[2*i]=1.0e30; //设置坐标未知点标志
s=0;
while(1)
{
s++;
for(i=0;i<=Lnumber-1;i++)
{
k=dir1[i];
if(XY[2*k]<1.0e20)continue;
for(j=dir0[i];j<dir0[i+1];j++)
{
if(XY[2*k]<1.0e20)break;
k2=dir2[j];
xB=XY[2*k2];
yB=XY[2*k2+1];
if(xB>1.0e29)continue;
for(jj=j+1;jj<dir0[i+1];jj++)
{
k1=dir2[jj];
xA=XY[2*k1];
yA=XY[2*k1+1];
if(xA>1.0e29)continue;
A=getangle(k1,k,k2);
B=getangle(k2,k1,k);
C=L[jj]-L[j];
if(A<0.0||B<0.0)continue;
sinA=sin(A); sinB=sin(B); sinC=sin(C);
if(C>PAI)
{
cosB=cos(B);
xk=xB+((xA-xB)*cosB-(yA-yB)*sinB)*sinA/sinC;
yk=yB+((xA-xB)*sinB+(yA-yB)*cosB)*sinA/sinC;
}
else
{
cosA=cos(A);
xk=xA+((xB-xA)*cosA+(yB-yA)*sinA)*sinB/sinC;
yk=yA+((yB-yA)*cosA-(xB-xA)*sinA)*sinB/sinC;
}
XY[2*k]=xk;
XY[2*k+1]=yk;
break;
} //jj
} //j
}
f=1;
for(i=0;i<Pnumber;i++)if(XY[2*i]>1.0e29){ f=0; break; }
if(f)break;
int unPnumber=Pnumber-knPnumber;
if(s>unPnumber && !f)
{
printf("\n有下列点无法计算出坐标:\n");
for(i=0;i<Pnumber;i++)
if(XY[2*i]>1.0e29)printf("\n%s",Pname[i]);
exit(0);
getch();
}
}
fprintf(resultfp,"\n近似坐标\n");
for(i=0;i<=Pnumber-1;i++)
{
double xi=XY[2*i];
double yi=XY[2*i+1];
fprintf(resultfp,"\n%2d %3s ",i+1,Pname[i]);
fprintf(resultfp,"%14.3lf%12.3lf", xi,yi+500000.0);
}
fflush(resultfp);
}
///////////////////////////////////////////////
//
// 计算方向误差方程式ab系数
//
///////////////////////////////////////////////
double ca_ab(int k1,int k2,double A[],int Ain[])
{
double dx=XY[2*k2]-XY[2*k1];
double dy=XY[2*k2+1]-XY[2*k1+1];
double s2=dx*dx+dy*dy;
A[0]= dy/s2*ROU; Ain[0]=2*k1;
A[1]=-dx/s2*ROU; Ain[1]=2*k1+1;
A[2]=-dy/s2*ROU; Ain[2]=2*k2;
A[3]= dx/s2*ROU; Ain[3]=2*k2+1;
double T=atan2(dy,dx);
if(T<0.0)T=T+2.0*PAI;
return T;
}
///////////////////////////////////////////////
// 计算边长误差方程系数
///////////////////////////////////////////////
double ca_cd(int k1,int k2,double A[],int Ain[])
{
double dx=XY[2*k2]-XY[2*k1];
double dy=XY[2*k2+1]-XY[2*k1+1];
double s=sqrt(dx*dx+dy*dy);
A[0]=-dx/s; Ain[0]=2*k1;
A[1]=-dy/s; Ain[1]=2*k1+1;
A[2]= dx/s; Ain[2]=2*k2;
A[3]= dy/s; Ain[3]=2*k2+1;
return s;
}
/////////////////////////////////////////////
// 一个误差方程式的累加项计算
void caATPAi(double A[],int Ain[],double pk,double Lk)
{
for(int i=0;i<=4;i++)
{
int k1=Ain[i];
double ai=A[i];
ATPL[k1]-=pk*ai*Lk;
for(int j=0;j<=i;j++)
{
int k2=Ain[j];
double aj=A[j];
ATPA[ij(k1,k2)]+=pk*ai*aj;
}
}
}
/////////////////////////////////////////////
//
// 组 成 法 方 程 式
//
/////////////////////////////////////////////
void caATPA(void)
{
double Li,A[5];
int Ain[5];
int t=2*Pnumber+Lnumber;
int tt=t*(t+1)/2;
for(int i=0;i<=tt-1;i++) ATPA[i]=0.0;
for(i=0;i<=t-1;i++) ATPL[i]=0.0;
double Pi=1.0/(ma*ma);
A[4]=-1.0;
for( i=0;i<=Lnumber-1;i++)
{
int k1=dir1[i];
Ain[4]=2*Pnumber+i;
double z;
for(int j=dir0[i];j<dir0[i+1];j++)
{
int k2=dir2[j];
double T12=ca_ab(k1,k2,A,Ain);
if(j==dir0[i])z=T12; //定向角近似值
Li=T12-z;
if(Li<0.0)Li+=2.0*PAI;
Li=(Li-L[j])*ROU;
caATPAi(A,Ain,Pi,Li);
V[j]=Li; //自由项放在V数组里，计算残差时使用
}
}
A[4]=0.0;
Ain[4]=0;
for(i=0;i<=Snumber-1;i++)
{
int k1=SSS[i].k1;
int k2=SSS[i].k2;
double m=mS1+mS2*SSS[i].L;
Pi=1.0/(m*m+1.0e-15);
double S12=ca_cd(k1,k2,A,Ain);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -