?? adjpro0504.cpp
字號:
cout<<endl<<endl;
cout<<endl<<endl<<" 未知數解："<<endl;
matdis(aa.X,aa.n);
cout<<endl<<endl<<" 未知數協因數陣："<<endl;
matdis(aa.QXX,aa.n,aa.n);
cout<<endl<<endl<<" 改正數："<<endl;
matdis(aa.V,aa.m);
cout<<endl<<endl<<" 單位權中誤差：+-"<<aa.m0<<endl;
}
//************************************************************************************************
int foutadj(adj &aa, char *name)
{
ofstream on(name);
if(!on) return 0;
// 輸出平差項目名
on<<aa.name<<endl<<endl;
// output observation number:
on<<aa.m<<endl;
// output unknown data number:
on<<aa.n<<endl<<endl;
//out unknown data results
for(int i=0;i<aa.n;i++)
on<<aa.X[i][0]<<endl;
on<<endl;
// 未知數協因數陣："<<endl;
for(i=0;i<aa.n;i++)
{
for(int j=0;j<aa.n;j++)
on<<aa.QXX[i][j]<<" ";
on<<endl;
}
on<<endl;
//output 改正數：;
for(i=0;i<aa.m;i++)
on<<aa.V[i][0]<<endl;
on<<endl;
on<<aa.m0<<endl;
on.close();
return 1;
}
//************************************************************************************************
int rubust(adj &a) // 定義抗差估計
{
// 1.最小二成平差
if(!doadj(a)) return 0;
double *xold,*xnew,small(0),P[MAX],n0(0);
xold=new double[a.n];
int num(0);
xnew=new double[a.n];
for(int i=0;i<a.m;i++)
P[i]=a.P[i][i];
do{
small=0;
for(i=0;i<a.n;i++) *(xold+i)=a.X[i][0];
// 2.等權處理
for(i=0;i<a.m;i++)
{
if(fabs(a.V[i][0])>2.5*a.m0) { a.P[i][i]=0;n0++;}
else if(fabs(a.V[i][0])<1.8*a.m0*sqrt(1/a.P[i][i])) ;
else a.P[i][i]=a.P[i][i]/1.8/a.m0;
}
cout<<"a.m0="<<a.m0<<endl;
matdis(a.P,a.m,a.m);
doadj(a);
for(i=0;i<a.n;i++)
{ *(xnew+i)=a.X[i][0]-*(xold+i);
if(small<fabs(*(xnew+i))) small=fabs(*(xnew+i));
}
num++;
}while(small>0.0000000001);
if(n0) a.m0*=sqrt((a.m-a.n)/(a.m-a.n-n0));
delete [] xold;
delete [] xnew;
return 1;
}
//************************************************************************************************
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; // 驗前單位權中誤差
Hpnt Pt[MAX]; // 高程點數組
line L[MAX]; // 水準線路數組
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.obnum;
in>>a.allpnum;
in>>a.fixpnum;
in>>a.m0;
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++;}
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.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;
//*********************** 平差數據準備 ***************************
// 近似高程計算
do{n=0; // 近似高程標志，=1表示仍有點近似高程未知
for(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);
// 極大權法平差數據
a.aa.m=a.obnum; // 觀測值個數
a.aa.n=a.allpnum; // 未知點個數（極大權處理，含控制點）
// 觀測值權陣計算
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;
else a.aa.P[i][j]=1/a.L[i].length;
cout<<endl<<endl<<" 平差數據準備，請稍等...... "<<endl<<endl;
cout<<" 觀測值權陣計算結果 :"<<endl;matdis(a.aa.P,a.obnum,a.obnum);
cout<<endl;
// 誤差方程系數陣計算
for(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;
// 誤差方程常數項計算
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;
// 平差數據保存在date.txt文件中
ofstream out("date.txt");
out<<a.netname<<endl;
out<<a.obnum<<endl;
out<<a.allpnum<<endl;
matout(a.aa.A,a.aa.m,a.aa.n,out);
out<<endl;
matout(a.aa.P,a.aa.m,a.aa.m,out);
out<<endl;
matout(a.aa.l,a.aa.m,out);
out.close();
return 1;
}
//************************************************************************************************
int Hnetadj(Hnet &a,char *outfile) // 高程網平差函數
{
fsetadj(a.aa,"date.txt");
cout<<endl<<" 正在平差計算，請稍等...... "<<endl<<endl;
doadj(a.aa,a.fixpnum,0); // 極大權法最小二乘平差：
for(int 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]);
}
// 結果屏幕輸出：
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;
// 平差結果文件保存
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.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( Hnet &a,char *infilename,char* outfilename)
{
if(finHnet(a,infilename)){
Hnetadj(a,outfilename);
return 1;
}
else return 0;
}
// ***********************************************************************************************
// *********************************** 平面網平差******** ****************************************
double d_h(double angle) //角度化弧度
{
double a,b;
angle=modf(angle,&a);
angle=modf(angle*100.0,&b);
return (a+b/60.0+angle/36.0)*(PI+3.0E-16)/180.0;
}
// ***********************************************************************************************
double h_d(double angle) //弧度化角度
{
double a,b,c;
angle=modf(angle*180.0/(PI-3.0E-16),&a);
angle=modf(angle*60.0,&b);
angle=modf(angle*60.0,&c);
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -