//#define MAX 100
#define PI 3.14159265358979312
#define rou (180.0*60*60/PI)
#include"MyClassLibrary0504.h"
#include<iostream.h>
//************************************************************************************************
class adj : public comadj{
// adjustment model :
// V=AX-L P
// A[m][n], P[m][m], L[m][1]
public:
adj(){}
adj(int k,int w):comadj(k,w) {}
void ksetadj();
int fsetadj(char *name) ; //file_data input
int doadj(int known,int r); // 极大权法最小二乘平差,known--控制点已知数据个数;
// r－固定数据个数＋测站数（平面网）
int doadj();
int rubust();
void adjdis(); // 屏幕输出
int foutadj(char *name); // 文件输出
private:
};
//*******************************************************************************************
int adj::doadj() // 普通最小二乘平差
{
MAT APA,AT;
AT=A.T();
APA=AT*P*A;
N_1=APA.inverse1();
int flag;
if(APA.R()==APA.GetRow())
flag=1;
else flag=0;
if(flag!=1)
{
this->flag=0;
return 0;
}
MAT AX=A.T()*P*l;
X=N_1*AX;
AX=A*X;
for(int i=0;i<m;i++)
V=AX-l;
MAT VPV=V.T()*P*V;
double cc=VPV.GetElem(0,0);
this->m0=sqrt(cc/(m-n));
this->flag=1;
return 1;
}
//************************************************************************************************
void adj::ksetadj() //keyboard and screen input
{
cout<<" input the object name: ";
cin>>this->name;cout<<endl;
cout<<" input XYobservation number: ";
cin>>this->m;cout<<endl;
cout<<" input unknown data number: ";
cin>>n;cout<<endl;
t=m-n;
A.SetRow(m); A.SetRank(n);
P.SetRow(m); P.SetRank(m);
N_1.SetRow(n); N_1.SetRank(n);
X.SetRow(n); X.SetRank(1);
l.SetRow(m); l.SetRank(1);
V.SetRow(m); V.SetRank(1);
mX.SetRow(n); mX.SetRank(1);
N.SetRow(n); N.SetRank(n);
cout<<" input data of mat A["<<m<<"]["<<n<<"]"<<endl;
for(int i=0;i<this->m;i++)
for(int j=0;j<this->n;j++)
{double cc; cin>>cc;A.SetElem(i,j,cc);}
cout<<" input data of mat P["<<this->m<<"]["<<this->m<<"]"<<endl;
for(i=0;i<this->m;i++)
for(int j=0;j<this->m;j++)
{double cc; cin>>cc;P.SetElem(i,j,cc);}
cout<<" input data of mat L["<<this->m<<"][1]"<<endl;
for(i=0;i<this->m;i++)
{double cc; cin>>cc;l.SetElem(i,0,cc);}
}
//************************************************************************************************
int adj::fsetadj(char *name) //keyboard and screen input
{
ifstream in(name,ios::nocreate);
if(!in) return 0;
// 输入平差项目名
in>>this->name;
// input XYobservation number:
in>>this->m;
// input unknown data number:
in>>this->n;
t=m-n;
A.SetRow(m); A.SetRank(n);
P.SetRow(m); P.SetRank(m);
N_1.SetRow(n); N_1.SetRank(n);
X.SetRow(n); X.SetRank(1);
l.SetRow(m); l.SetRank(1);
V.SetRow(m); V.SetRank(1);
mX.SetRow(n); mX.SetRank(1);
N.SetRow(n); N.SetRank(n);
// input data of mat A[m][n]:
for(int i=0;i<this->m;i++)
for(int j=0;j<this->n;j++)
{double cc; in>>cc;A.SetElem(i,j,cc);}
// input data of mat P[m][m]:
for(i=0;i<this->m;i++)
for(int j=0;j<this->m;j++)
{double cc; in>>cc;P.SetElem(i,j,cc);}
// input data of mat L[m][1]:
for(i=0;i<this->m;i++)
{double cc; in>>cc;l.SetElem(i,0,cc);}
in.close();
return 1;
}
//*******************************************************************************************
int adj::doadj(int known,int r) // 极大权法最小二乘平差,known--控制点已知数据个数;
{ // r－固定数据个数＋测站数（平面网）
MAT APA,AT;
AT=A.T();
APA=AT*P*A;
double add(0);
for(int i=0;i<this->n;i++)
if(add<=APA.GetElem(i,i)) add=APA.GetElem(i,i);
add*=1000;
for(i=0;i<known;i++) // 对已知点方法程系数阵的处理：
APA.SetElem(i,i,APA.GetElem(i,i)+add); // 控制点点号对应矩阵主元加平均权的10000倍
N_1=APA.inverse1();
if(N_1.R()==N_1.GetRow()) flag=1;
else flag=0;
if(flag!=1) // 不可求逆的判断及处理
{ this->flag=0; return 0;}
MAT AX; // 极大权平差过程
AX=AT*P*l;
X=N_1*AX;
AX=A*X;
V=AX-l;
MAT VPV=V.T()*P*V;
double cc=VPV.GetElem(0,0); // VTPV计算
this->m0=sqrt(cc/(this->m-this->n+known-r)); // 单位权中误差计算，必要观测数为（this->n-knownp）
this->flag=1; // 平差完毕
return 1;
}
//*******************************************************************************************
void adj::adjdis()
{
cout<<" 平差项目名 ："<<this->name<<endl<<endl;
cout<<" 误差方程阵："<<endl;
this->A.print();
cout<<endl<<endl<<" 观测值权矩阵："<<endl;
this->P.print();
cout<<endl<<endl<<" 常数项："<<endl;
l.print();
cout<<endl<<endl;
cout<<endl<<endl<<" 未知数解："<<endl;
X.print();
cout<<endl<<endl<<" 未知数协因数阵："<<endl;
N_1.print();
cout<<endl<<endl<<" 改正数："<<endl;
V.print();
cout<<endl<<endl<<" 单位权中误差：+-"<<this->m0<<endl;
}
//************************************************************************************************
int adj::foutadj(char *name)
{
ofstream on(name);
if(!on) return 0;
// 输出平差项目名
on<<this->name<<endl<<endl;
// output XYobservation number:
on<<this->m<<endl;
// output unknown data number:
on<<this->n<<endl<<endl;
//out unknown data results
for(int i=0;i<n;i++)
on<<X.GetElem(i,0)<<endl;
on<<endl;
// 未知数协因数阵："<<endl;
for(i=0;i<n;i++)
{
for(int j=0;j<n;j++)
on<<N_1.GetElem(i,j)<<" ";
on<<endl;
}
on<<endl;
//output 改正数：;
for(i=0;i<m;i++)
on<<V.GetElem(i,0)<<endl;
on<<endl;
on<<m0<<endl;
on.close();
return 1;
}
//************************************************************************************************
int adj :: rubust() // 定义抗差估计
{
// 1.最小二乘平差
if(!doadj()) return 0;
double *xold,*xnew,small(0),n0(0);
MAT P(m,1);
xold=new double[n];
int num(0);
xnew=new double[n];
for(int i=0;i<m;i++)
P.SetElem(i,0,this->P.GetElem(i,i));
do{
small=0;
for(i=0;i<n;i++) *(xold+i)=X.GetElem(i,0);
// 2.等权处理
for(i=0;i<m;i++)
{
if(fabs(V.GetElem(i,0))>2.5*m0) { this->P.SetElem(i,i,0);n0++;}
else if(fabs(this->V.GetElem(i,0))<1.8*this->m0*sqrt(1/this->P.GetElem(i,i))) ;
else this->P.SetElem(i,i,this->P.GetElem(i,i)/1.8/this->m0);
}
cout<<"this->m0="<<this->m0<<endl;
P.print();
doadj();
for(i=0;i<this->n;i++)
{ *(xnew+i)=this->X.GetElem(i,0)-*(xold+i);
if(small<fabs(*(xnew+i))) small=fabs(*(xnew+i));
}
num++;
}while(small>0.0000000001);
if(n0) this->m0*=sqrt((this->m-this->n)/(this->m-this->n-n0));
delete [] xold;
delete [] xnew;
return 1;
}
//**********************************************************************************
//**************** 高程网平差程序 ****************************************
class Hpnt{ // 高程点结构
public:
char name[20];
double H;
double H0;
double mH;
int fixed; // known=1;else =0;
int i; // point number
};
class line{ // 水准线路结构
public:
Hpnt *startp;
Hpnt *endp;
double length;
double h;
int style; // 水准测量=1;三角高程=2
};
class Hnet{ // 高程网结构
public:
Hnet(char *fname); // 文件输入高程网函数
Hnet(){obnum=allpnum=fixpnum=0;}
Hnet(int allp,int fixp,int obn):aa(obn,allp-fixp)
{
allpnum=allp;
fixpnum=fixp;
obnum=obn;
Pt=new Hpnt[allpnum];
L=new line[obnum];
}
评论0