水准网间接平差程序设计(C++)

//////////////////////////////////////////////////////visualC++6.0编译通过///////////////////////////////////////////////////////////////////////////////////////////////////////////参考资料////部分网络资料////宋力杰《测量平差程序设计》////姚连壁《基于matlab的控制网平差程序设计》/////////////////////////////////////////////////////#includeiostream#includefstream#includestdlib.h#includemath.h#includeiomanipusingnamespacestd;//////////////////////////////////////////////////////////////////////////classclassSZWPC{private:intgcz_zs;//高差总数intszd_zs;//总点数intyz_szd_zs;//已知点数doublem_pvv;//[pvv]int*qsd_dh;//高差起点号int*zd_dh;//高差终点号char**dm;//点名地址数组double*gcz;//观测值数组double*szd_gc;//高程值数组double*P;//观测值的权double*ATPA,*ATPL;//法方程系数矩阵与自由项double*dX;//高程改正数、平差值double*V;//残差doublem_mu;//单位权中误差public:SZWPC();~SZWPC();intij(inti,intj);//对称矩阵下标计算函数boolinverse(doublea[],intn);//对称正定矩阵求逆(仅存下三角元素)(参考他人）voidinputdata(char*datafile);//输入原始数据函数intdm_dh(char*name);//点名转点号voidca_H0();//近似高程计算函数voidca_ATPA();//法方程组成函数voidca_dX();//高程平差值计算函数voidprintresult(char*resultfile);//精度估计与平差值输出函数doubleca_V();//残差计算函数voidzxecpc(char*resultfile);//最小二乘平差函数};////////////////////////////////////////////////////////////////////////构造函数SZWPC::SZWPC(){gcz_zs=0;szd_zs=0;yz_szd_zs=0;}////////////////////////////////////////////////////////////////////////析构函数SZWPC::~SZWPC(){if(gcz_zs0){delete[]qsd_dh;delete[]zd_dh;delete[]gcz;delete[]P;delete[]V;}if(szd_zs0){delete[]szd_gc;delete[]ATPA;delete[]ATPL;delete[]dX;for(inti=0;iszd_zs;i++)if(dm[i]!=NULL)delete[](dm[i]);delete[]dm;}}////////////////////////////////////////////////////////////////////////////对称矩阵下标计算函数intSZWPC::ij(inti,intj){return(i=j)?i*(i+1)/2+j:j*(j+1)/2+i;}////////////////////////////////////////////////////////////////////////////对称正定矩阵求逆(仅存下三角元素)(参考他人）boolSZWPC::inverse(doublea[],intn){double*a0=newdouble[n];for(intk=0;kn;k++){doublea00=a[0];if(a00+1.0==1.0){delete[]a0;returnfalse;}for(inti=1;in;i++){doubleai0=a[i*(i+1)/2];if(i=n-k-1)a0[i]=-ai0/a00;elsea0[i]=ai0/a00;for(intj=1;j=i;j++){a[(i-1)*i/2+j-1]=a[i*(i+1)/2+j]+ai0*a0[j];}}for(i=1;in;i++){a[(n-1)*n/2+i-1]=a0[i];}a[n*(n+1)/2-1]=1.0/a00;}delete[]a0;returntrue;}///////////////////////////////////////////////////////////////////////原始数据输入函数voidSZWPC::inputdata(char*datafile){ifstreaminfile(datafile,ios::in);if(!infile){cerrOpenerror!endl;}infilegcz_zsszd_zsyz_szd_zs;intunPnumber=szd_zs-yz_szd_zs;szd_gc=newdouble[szd_zs];dX=newdouble[szd_zs];ATPA=newdouble[szd_zs*(szd_zs+1)/2];ATPL=newdouble[szd_zs];qsd_dh=newint[gcz_zs];zd_dh=newint[gcz_zs];gcz=newdouble[gcz_zs];V=newdouble[gcz_zs];P=newdouble[gcz_zs];dm=newchar*[szd_zs];for(inti=0;iszd_zs;i++){dm[i]=NULL;//dm_dh函数根据dm[i]是否为NULL确定dm[i]是否为点名地址}charbuffer[128];//临时数组，保存从文件中读到的点名for(i=0;i=yz_szd_zs-1;i++)//读取已知高程数据{infilebuffer;intc=dm_dh(buffer);infileszd_gc[i];}for(i=0;igcz_zs;i++)//读取观测数据{infilebuffer;//读取高程起点名qsd_dh[i]=dm_dh(buffer);infilebuffer;//读取高程终点zd_dh[i]=dm_dh(buffer);infilegcz[i]P[i];//读取高差值与路线长度P[i]=1.0/P[i];//线路长转化为观测值的权}infile.close();}////////////////////////////////////////////////////////////////////点名转点号，返回点名对应的点号intSZWPC::dm_dh(char*name){for(inti=0;iszd_zs;i++){if(dm[i]!=NULL){if(strcmp(name,dm[i])==0)returni;//将待查点名与已经存入点名数组的点名比较，若存在返回点号}else{intlen=strlen(name);//判断点名长度dm[i]=newchar[len+1];//为点名申请存储空间strcpy(dm[i],name);//待查点名是一个新的点名，将新点名的地址放到dm数组中returni;//返回点号}}return-1;//dm数组已经存满，且没有待查点名}////////////////////////////////////////////////////////////////////////////高程近似值计算voidSZWPC::ca_H0(){for(inti=yz_szd_zs;iszd_zs;i++)szd_gc[i]=-10000.9;//为计算机设置辨别未知高程点的标志for(intj=1;;j++){intk=0;//计算出近似高程的点数for(i=0;igcz_zs;i++){intk1=qsd_dh[i];//高差起点号intk2=zd_dh[i];//高差终点号if(szd_gc[k1]-10000.0&&szd_gc[k2]-10000.0)//k1点高程或高程近似值已知，k2点高程或高程近似值未知{szd_gc[k2]=szd_gc[k1]+gcz[i];//计算近似高程k++;}elseif(szd_gc[k1]-10000.0&&szd_gc[k2]-10000.0)//k2点高程或高程近似值已知，k1点高程或高程近似值未知{szd_gc[k1]=szd_gc[k2]-gcz[i];//计算近似高程k++;}}if(k==(szd_zs-yz_szd_zs))break;//所有的近似高程计算完成，退出}}////////////////////////////////////////////////////////////////////////////组成法方程voidSZWPC::ca_ATPA(){//intt=szd_zs;for(inti=0;iszd_zs*(szd_zs+1)/2;i++)ATPA[i]=0.0;//赋初值for(i=0;iszd_zs;i++)ATPL[i]=0.0;//赋初值for(intk=0;kgcz_zs;k++){inti=qsd_dh[k];//获取点号intj=zd_dh[k];//获取点号doublePk=P[k];//获取权值doublelk=gcz[k]-(szd_gc[j]-szd_gc[i]);//获得第k个自由项ATPL[i]-=Pk*lk;//获得法方程自由项ATPL[j]+=Pk*lk;ATPA[ij(i,i)]+=Pk;//获得法方程系数矩阵ATPA[ij(j,j)]+=Pk;ATPA[ij(i,j)]-=Pk;}}////////////////////////////////////////////////////////////////////////////高程平差值计算voidSZWPC::ca_dX(){for(inti=0;iyz_szd_zs;i++)ATPA[ij(i,i)]=1.0e30;//处理已知点if(!inverse(ATPA,szd_zs))//矩阵求逆{cerr法方程系数矩阵降秩!endl;//矩阵为奇异矩阵，无法求逆exit(0);//退出程序}for(i=0;iszd_zs;i++)//计算高程改正数{doublexi=0.0;for(intj=0;jszd_zs;j++){xi+=ATPA[ij(i,j)]*ATPL[j];}dX[i]=xi;szd_gc[i]+=xi;//计算高程平差值}}////////////////////////////////////////////////////////////////////////////残差计算doubleSZWPC::ca_V(){doublepvv=0.0;for(inti=0;i=gcz_zs-1;i++){intk1=qsd_dh[i];intk2=zd_dh[i];V[i]=szd_gc[k2]-szd_gc[k1]-gcz[i];pvv+=V[i]*V[i]*P[i];}re