声明：本程序未经过工程实际测试，仅用于测绘程序的学习交流，望您留下对本程序宝贵的意见，谢谢！

适用于：水准网平差  三角网平差  水准网和三角网联合平差

//************************************************************************************************
void matout(double A[][1],int n,ofstream out)         // 向文件输出矩阵
{
for(int i=0;i<n;i++)
      out<<"     "<<A[i][0]<<endl;
}
//************************************************************************************************
void matout(double A[][MAX],int n,int m,ofstream out) // 向文件输出矩阵
{//1.set B[][] I;
   for(int i=0;i<n;i++)
   {  out<<"     ";
   for(int j=0;j<m;j++)
   out<<A[i][j]<<"  ";
   out<<endl;
   }
}
//**********************************************************************************
//****************       高程网平差程序     ****************************************
struct Hpnt{                            //  高程点结构
 char name[20];
      double H;
   double H0;
      double mH;
   int fixed; // known=1;else =0;
   int i;     // point number
};

struct line{                            //  水准线路结构            
      Hpnt *startp;
      Hpnt *endp;
      double length;
   double h;
   int style;  // 水准测量=1;三角高程=2
};

struct Hnet{                            //  高程网结构
 char netname[40];                   //  网名
      int allpnum;                      //  总点数
   int fixpnum;                      //  控制点个数
   int obnum;                        //  观测值个数
      double m0;                        //  验前单位权中误差或水准单位权中误差
   double m1;                         //  三角网单位权中误差
   Hpnt Pt[MAX];                     //  高程点数组
   line L[MAX];                      //  水准线路数组
   int  style;                       //  高程网类型（=1，水准网；=2，三角高程网；=3，混合网）
   adj aa;                           //  平差结构
};
//************************************************************************************************
int finHnet(Hnet &a,char *fname)        //  文件输入高程网数据函数
{
 ifstream in(fname,ios::nocreate);   // 建立文件流,并与输入文件名建立关联
 if(!in) {cout<<fname<<" error: file does not exist!   "<<endl; return 0;}
 //  文件存在判断
 char name[20];
    in>>a.netname;  
 in>>a.style;          
 in>>a.obnum;
 in>>a.allpnum;
 in>>a.fixpnum;
 in>>a.m0;
 if(a.style==3)
 in>>a.m1;
 else
  ;
 int n(a.fixpnum);                   // n为已输入名字的点的个数
// 输入控制点信息 
 for(int i=0;i<a.fixpnum;i++)
 {
  in>>a.Pt[i].name>>a.Pt[i].H;
        a.Pt[i].fixed=1;                // 控制点标记
  a.Pt[i].H0=a.Pt[i].H;
  a.Pt[i].mH=0;
  a.Pt[i].i=i;                    // 控制点编号，从0到a.fixpnum-1
 }
 for(i=a.fixpnum;i<a.allpnum;i++)    // 输入未知点相关信息（名字在后面输入）   
 {
     a.Pt[i].fixed=0;                // 未知点标记
  a.Pt[i].H=0;a.Pt[i].H0=-PI;
  a.Pt[i].mH=0;
  *(a.Pt[i].name)=0;
  a.Pt[i].i=i;                    // 为未知点编号,从a.fixpnum到a.allpnum-1
 }
// 输入观测值    
 for(i=0;i<a.obnum;i++)
 {int t=0;                           // 点名比较标志
   in>>name;                         // 输入起点名
     for(int k=0;k<n;k++)
    if(strnicmp(name,a.Pt[k].name,20)==0)
    {
     a.L[i].startp=&(a.Pt[k]); // 找到同名点，起点指针指向该点
     t++;                      // 找到标志
    }
  if(t==0) {strcpy(a.Pt[n].name,name);
   a.L[i].startp=&(a.Pt[n]);   // 找不到同名点，该名输给新点
      n++;}
   in>>name; t=0;                    // 输入终点名，操作过程同上   
       for(k=0;k<n;k++)
    if(strnicmp(name,a.Pt[k].name,20)==0)
    {
      a.L[i].endp=&(a.Pt[k]);
       t++;}
  if(t==0) {strcpy(a.Pt[n].name,name);
    a.L[i].endp=&(a.Pt[n]);
    n++;}
      if(a.style==3)
   in>>a.L[i].style;                   //输入线路类型
      in>>a.L[i].h;                     // 输入线路高差
   in>>a.L[i].length;                // 输入线路长度
 }
  if(n!=a.allpnum)
  {cout<<fname<<" error: file provide not correct point number !  "<<endl;
  return 0;}
   //  文件正确性判断
in.close();                             // 关闭输入流及关联文件

// 向屏幕输出原始平差文件
cout<<"     平差文件 "<<fname<<" 数据输入结果："<<endl;
cout<<"     "<<a.netname<<"  "<<a.style<<"  "<<a.obnum<<"  "<<a.allpnum<<"   "
<<a.fixpnum<<"   "<<a.m0<<"     "<<a.m1<<endl<<endl;
for(i=0;i<a.fixpnum;i++)                // 控制点数据
cout<<"     "<<a.Pt[i].name<<"   "<<a.Pt[i].H<<endl;
cout<<endl;
for(i=0;i<a.obnum;i++)                  // 水准线路数据
cout<<"     "<<a.L[i].startp->name<<"  "<<a.L[i].endp->name<<"  "<<a.L[i].style<<"     "<<a.L[i].h
<<"  "<<a.L[i].length<<endl; 
return 1;
}

//************************************************************************************************
int kinHnet(Hnet &a)                    //  键盘输入高程网数据函数
{

 cout<<"......输入高程网总体信息......"<<endl;
 char name[20];
 cout<<"    输入高程网网名：  ";
    cin>>a.netname;
 cout<<"    选择高程网类型：  "<<endl;
 cout<<"    1，水准网；  "<<endl;
 cout<<"    2，三角高程网;  "<<endl;
 cout<<"    3、混合网平差;   "<<endl;
 cin>>a.style;
    cout<<"    输入高程路线个数：  ";
 cin>>a.obnum;
 cout<<"    输入高程网总点数：  ";
 cin>>a.allpnum;
 cout<<"    输入已知高程点数：  ";
 cin>>a.fixpnum;
 if(a.style==3)
    {
 cout<<"    输入水准高程网验前单位权中误差：  ";
 cin>>a.m0;
 cout<<"    输入三角高程网验前单位权中误差：  ";
 cin>>a.m1;
 }
 else
 {cout<<"    输入高程网验前单位权中误差：  ";
 cin>>a.m0;
 }
 
 int n(a.fixpnum);                   // n为已输入名字的点的个数
// 输入控制点信息
 cout<<"......输入已知高程点数据（格式：点名  高程）：......"<<endl;
 for(int i=0;i<a.fixpnum;i++)
 {
  cin>>a.Pt[i].name>>a.Pt[i].H;
        a.Pt[i].fixed=1;                // 控制点标记
  a.Pt[i].H0=a.Pt[i].H;
  a.Pt[i].mH=0;
  a.Pt[i].i=i;                    // 控制点编号，从0到a.fixpnum-1
 }
 for(i=a.fixpnum;i<a.allpnum;i++)    // 设置未知点相关信息（名字在后面输入）   
 {
     a.Pt[i].fixed=0;                // 未知点标记
  a.Pt[i].H=0;a.Pt[i].H0=-PI;
  a.Pt[i].mH=0;
  *(a.Pt[i].name)=0;
  a.Pt[i].i=i;                    // 为未知点编号,从a.fixpnum到a.allpnum-1
 }
// 输入观测值
  cout<<"......输入高程线路信息（起点名  终点名  【混合网时输入】线路类型（1水准线路2三角线路）  高差   线路长度）：......"<<endl;

 for(i=0;i<a.obnum;i++)
 {int t=0;                           // 点名比较标志
   cin>>name;                         // 输入起点名
     for(int k=0;k<n;k++)
    if(strnicmp(name,a.Pt[k].name,20)==0)
    {
     a.L[i].startp=&(a.Pt[k]); // 找到同名点，起点指针指向该点
     t++;                      // 找到标志
    }
  if(t==0) {strcpy(a.Pt[n].name,name);
   a.L[i].startp=&(a.Pt[n]);   // 找不到同名点，该名输给新点
      n++;}
   cin>>name; t=0;                    // 输入终点名，操作过程同上   
       for(k=0;k<n;k++)
    if(strnicmp(name,a.Pt[k].name,20)==0)
    {
      a.L[i].endp=&(a.Pt[k]);
       t++;}
  if(t==0) {strcpy(a.Pt[n].name,name);
    a.L[i].endp=&(a.Pt[n]);
    n++;}
 if(a.style==3)
  cin>>a.L[i].style;               //输入线路类型
      cin>>a.L[i].h;                     // 输入线路高差
   cin>>a.L[i].length;                // 输入线路长度
 }
  if(n!=a.allpnum) {cout<<" input data error:  providing not correct point number !  "<<endl;
  return 0;
  }

//////////数据保存///////////////  
  cout<<"   是否保存数据文件："<<endl<<"     1，保存；"<<endl<<"     2，不保存。"<<endl;
  int bc;
  cin>>bc;
  if(1==bc)
  {
   cout<<"    输入保存文件名："<<endl;
   char fname[100];
   cin>>fname;
   ofstream out(fname);
   out<<"     "<<a.netname<<"  "<<a.style<<"  "<<a.obnum<<"  "<<a.allpnum<<"   "
         <<a.fixpnum<<"   "<<a.m0<<"      "<<a.m1<<endl<<endl;
       for(int i=0;i<a.fixpnum;i++)                // 控制点数据
        out<<"     "<<a.Pt[i].name<<"   "<<a.Pt[i].H<<endl;
        out<<endl;
       for(i=0;i<a.obnum;i++)                  // 水准线路数据
       out<<"     "<<a.L[i].startp->name<<"  "<<a.L[i].endp->name<<"  "<<a.L[i].style<<"        "<<a.L[i].h
        <<"  "<<a.L[i].length<<endl;
    out.close();
 
  }

// 向屏幕输出原始平差文件
cout<<"     平差数据输入结果："<<endl;
cout<<"     "<<a.netname<<"  "<<a.style<<"  "<<a.obnum<<"  "<<a.allpnum<<"   "
<<a.fixpnum<<"   "<<a.m0<<endl<<endl;
for(i=0;i<a.fixpnum;i++)                // 控制点数据
cout<<"     "<<a.Pt[i].name<<"   "<<a.Pt[i].H<<endl;
cout<<endl;
for(i=0;i<a.obnum;i++)                  // 水准线路数据
cout<<"     "<<a.L[i].startp->name<<"  "<<a.L[i].endp->name<<"  "<<a.L[i].h
<<"  "<<a.L[i].length<<endl; 
return 1;
}

//************************************************************************************************
int Hnetadj(Hnet &a,char *outfile)       // 高程网平差函数
{
  // 1.  近似高程计算
 int n;
   do{n=0;                     // 近似高程标志，=1表示仍有点近似高程未知
   for(int i=0;i<a.obnum; i++)          
    if( a.L[i].startp->H0==-PI || a.L[i].endp->H0==-PI) { n=1; break;}     
    for(i=0;i<a.obnum; i++)          // 由水准线路计算近似高程
    { 
      if( a.L[i].startp->H0!=-PI && a.L[i].endp->H0==-PI)
    a.L[i].endp->H0=a.L[i].startp->H0+a.L[i].h;   
   if(a.L[i].startp->H0==-PI && a.L[i].endp->H0!=-PI)
          a.L[i].startp->H0=a.L[i].endp->H0-a.L[i].h;
    }                               
  }while(n==1);
cout<<endl<<endl<<"     平差数据准备，请稍等......    "<<endl<<endl;
//2. 极大权法平差数据 
  a.aa.m=a.obnum;                       // 观测值个数
  a.aa.n=a.allpnum;                     // 未知点个数（极大权处理，含控制点）
//3. 观测值权阵计算
  if(1==a.style)    // 水准网观测值定权
  for(int i=0;i<a.aa.m;i++)
  for(int j=0;j<a.aa.m;j++)
   if(i!=j) a.aa.P[i][j]=0;
    else a.aa.P[i][j]=1/a.L[i].length;
  if(2==a.style)    // 三角高程网观测值定权
  for(int i=0;i<a.aa.m;i++)
  for(int j=0;j<a.aa.m;j++)
   if(i!=j) a.aa.P[i][j]=0;
   else a.aa.P[i][j]=1/a.L[i].length/a.L[i].length;
  if(3==a.style)
  {
     for(int i=0;i<a.aa.m;i++)
  { for(int j=0;j<a.aa.m;j++)
     if(a.L[i].style==1)
     a.aa.P[i][j]=1/a.L[i].length;
     else
      a.aa.P[i][j]=pow(a.m0,2)/(pow(a.m1,2)*pow(a.L[i].length,2));}
  for( i=0;i<a.aa.m;i++)
  for( int j=0;j<a.aa.m;j++)
            if(i!=j) a.aa.P[i][j]=0;
   
   
  }

            cout<<"      观测值权阵计算结果 :"<<endl;
   matdis(a.aa.P,a.obnum,a.obnum);
   cout<<endl;
//4. 误差方程系数阵计算
   for(int i=0;i<a.aa.m;i++)                // 根据水准线路号及起、终点号确定
    for(int j=0;j<a.aa.n;j++)
     a.aa.A[i][j]=0;
   for(i=0;i<a.aa.m;i++)
   {   a.aa.A[i][a.L[i].endp->i]=1;
         a.aa.A[i][a.L[i].startp->i]=-1;
   }
   cout<<"     "<<" 误差方程系数阵计算结果:  "<<endl;  
   matdis(a.aa.A,a.aa.m,a.aa.n);
   cout<<endl;
//5. 误差方程常数项计算
  double AX0;
  for(i=0;i<a.aa.m;i++)                   // V=AX-L
  {                                       // L=h-AX0
   AX0=0;
      for(int j=0;j<a.aa.n;j++)
    AX0+=a.aa.A[i][j]*a.Pt[j].H0;
      a.aa.l[i][0]=a.L[i].h-AX0;
  }
  cout<<endl;
  cout<<"     "<<" 误差方程常数项计算结果 ："<<endl;
  matdis(a.aa.l,a.aa.m);
  cout<<endl;

//6.极大权法平差
 cout<<endl<<"     正在平差计算，请稍等......    "<<endl<<endl;
  doadj(a.aa,a.fixpnum,0);               // 极大权法最小二乘平差：
  for( i=0;i<a.allpnum;i++)           // 未知点高程及精度计算
  {
   a.Pt[i].H=a.Pt[i].H0+a.aa.X[i][0];
      a.Pt[i].mH=a.aa.m0*sqrt(a.aa.QXX[i][i]);
  } 
//7. 结果屏幕输出：
 cout<<endl<<"     "<<" 验后单位权中误差：+-"<<a.aa.m0<<endl<<endl<<" 未知数改正数dX:  "<<endl;
  matdis(a.aa.X,a.aa.n);
cout<<endl<<"     "<<" 观测值改正数V:"<<endl;
matdis(a.aa.V,a.aa.m);
      cout<<endl<<"     "<<" 未知点高程及精度 ："<<endl;
     for(i=a.fixpnum;i<a.allpnum;i++)
        cout<<"     "<<a.Pt[i].name<<"高程="<<a.Pt[i].H<<"  +-"<<a.Pt[i].mH<<endl;
  cout<<endl;
//8.  平差结果文件保存
   ofstream out(outfile);
   if(!out) cout<<"can not open save file!"<<endl;
   out<<"         高程网"<<a.netname<<"平差计算   "<<endl<<endl;
   out<<"1. 原始数据：   "<<endl<<endl;
   out<<"    "<<a.netname<<"   "<<a.style<<"    "<<a.obnum<<"   "<<a.allpnum<<"  "<<a.fixpnum<<"  "<<a.m0<<endl;
   for(i=0;i<a.fixpnum;i++)
    out<<"    "<<a.Pt[i].name<<"   "<<a.Pt[i].H<<endl;
   out<<endl;
   for(i=0;i<a.obnum;i++)
    out<<"    "<<a.L[i].startp->name<<"   "<<a.L[i].endp->name<<"   "<<a.L[i].h
    <<"   "<<a.L[i].length<<endl;
   out<<endl;
   out<<"2. 平差数据："<<endl<<endl;
   out<<"   2.1  误差方程系数阵： "<<endl;
   matout(a.aa.A,a.aa.m,a.aa.n,out);
   out<<endl<<endl;
   out<<"   2.2  误差方程权阵： "<<endl;
   matout(a.aa.P,a.aa.m,a.aa.m,out);
   out<<endl<<endl;
   out<<"   2.3  未知点近似高程： "<<endl;
   for(i=a.fixpnum;i<a.allpnum;i++)
   {   out<<"     ";
    out<<a.Pt[i].name<<"   "<<a.Pt[i].H0<<"  ";
    out<<endl;
   }
   out<<endl<<endl;
   out<<"   2.4  误差方程常数项： "<<endl;
   matout(a.aa.l,a.aa.m,out);
   out<<endl;
   out<<"3. 平差结果 "<<endl<<endl;
   out<<"   3.1  观测值改正数V： "<<endl;
   matout(a.aa.V,a.aa.m,out);
   out<<endl<<endl;
   out<<"   3.2  单位权中误差m0：+-"<<a.aa.m0<<endl;
   out<<endl<<endl;
   out<<"   3.3  未知数改正数dH："<<endl;
   for(i=a.fixpnum;i<a.allpnum;i++)
   {   out<<"     ";
      out<<a.aa.X[i][0]<<"  ";
    out<<endl;
   }
   out<<endl<<endl;
   out<<"   3.4  未知点高程及精度： "<<endl;
   for(i=a.fixpnum;i<a.allpnum;i++)
   {   out<<"     ";
    out<<a.Pt[i].name<<"   "<<a.Pt[i].H<<"   +-"<<a.Pt[i].mH;
    out<<endl;
   }
      out<<endl<<endl;
   out<<"   3.5  未知参数协方差阵： "<<endl;
   for(i=a.fixpnum;i<a.allpnum;i++)
   {   out<<"     ";
    for(int j=a.fixpnum;j<a.allpnum;j++)
           out<<a.aa.m0*a.aa.m0*a.aa.QXX[i][j]<<"  ";
    out<<endl;
   }
   out.close();
  return 1;
}

//************************************************************************************************
int Hdoadj(char *infilename,char* outfilename)    // 文件输入数据，结果文件保存
{
 Hnet a;
    if(finHnet(a,infilename))
 {
      if(Hnetadj(a,outfilename))
      return 1;
 }
     return 0;
}

int Hdoadj(char* outfilename)        // 键盘输入高程网数据，结果文件保存
{
 Hnet a;
    if(kinHnet(a))
 {
     if(Hnetadj(a,outfilename))
     return 1;
 }
     return 0;
}
int Hadj()
{
  Hnet a;
  char infilename[100];
  char outfilename[100];
  cout<<endl<<"     ********************高程网平差Version 1.0********************* "<<endl<<endl;
   int style;
   // 1. 数据输入
   do{
    cout<<"    请选择数据输入方式：\n      1. 文件输入高程网数据； \n      2. 键盘输入高程网数据； \n      3. 退出。 "<<endl;
       cin>>style;
      switch(style)
   {
        case 1:
        cout<<"  输入文件名：";
  cin>>infilename;
        if(!finHnet(a,infilename))
  return 0;
     break;
        case 2:
   if(!kinHnet(a))
    return 0;
     break;
  case 3:
  return 0;
     break;
        default :
      cout<<"  错误：非法类型选择！   请重新选择 "<<endl;
   style=0;
   };
   }while(style==0);
    // 2.平差处理
         cout<<"\n\n  平差结果保存文件名：";
  cin>>outfilename;
  if(!Hnetadj(a,outfilename)) return 0;
 
  return 1;

}

//***********************************************************************************************
//************************************************************************************************

void main()
{
  Hadj();
}