#include
#include
#include
#define M 20
#define BM 8
#define GS 13
double JH[M],l[M],v[M],B[M][M],BT[M][M],P[M][M],BTP[M][M],BTPB[M][M],BTPBN[M][M],BTPL[M],x[M],H[M],h[M],Qhh[M][M];
double danweiquanzhongwucha,gcz_zhongwucha[M],H_zhongwucha[M];
int fg_JH[M],fg[M][M];
struct
{
int qianshidian;
int houshidian;
double gaocha;
double changdu;
}gaochaguancezhi[M];
void shuzufuchushizhi()
{
int i,j;
for (i=0;i<M;i++)
{
for (j=0;j<M;j++)
{
fg_JH[i]=fg[i][j]=0;
l[i]=P[i][j]=BTP[i][j]=BTPB[i][j]=BTPBN[i][j]=BTPL[i]=0;
}
}
}
void shujuduqu()
{
int i,a,b,no;
double t1,t2;
FILE *fp;
fp=fopen("水准网原始数据.txt","r");
if (fp==NULL)
{
printf("打开文件失败!\n");
exit(0);
}
else
//读取高差观测值
for (i=1;i<=GS;i++)
{
fscanf(fp,"%d%d%d%lf%lf",&no,&a,&b,&t1,&t2);
gaochaguancezhi[i].houshidian=a;
gaochaguancezhi[i].qianshidian=b;
gaochaguancezhi[i].gaocha=t1;
gaochaguancezhi[i].changdu=t2;
}
//读取起始高程
fscanf(fp,"%d%lf",&a,&t1);
JH[a]=t1;
fg_JH[a]=1;
fclose(fp);
}
//计算近似高程
void jisuangaochengjinsizhi()
{
int i,a,b;
double t1;
for (i=1;i<=GS;i++)
{
a=gaochaguancezhi[i].houshidian;
b=gaochaguancezhi[i].qianshidian;
t1=gaochaguancezhi[i].gaocha;
if (fg_JH[a]!=0&&fg_JH[b]==0)
{
JH[b]=JH[a]+t1;
fg_JH[b]=1;
}
if(fg_JH[a]==0&&fg_JH[b]!=0)
{
JH[a]=JH[b]-t1;
fg_JH[a]=1;
}
}
}
void cal_xishu()
{
int i,hs,qs,x1,x2;
double xhs,xqs;
for (i=1;i<=GS;i++)
{
x1=x2=1;
hs=gaochaguancezhi[i].houshidian;
qs=gaochaguancezhi[i].qianshidian;
if (hs==8)
{
x1=0;
}
if (qs==8)
{
x2=0;
}
xhs=x1*(-1);
xqs=1*x2;
B[i][hs]=xhs;
B[i][qs]=xqs;
}
}
//计算常数矩阵
void cal_l()
{
int i,qs,hs;
double t;
for (i=1;i<=GS;i++)
{
qs=gaochaguancezhi[i].qianshidian;
hs=gaochaguancezhi[i].houshidian;
t=gaochaguancezhi[i].gaocha;
l[i]=t-(JH[qs]-JH[hs]);
l[i]=l[i]*1000;
}
}
//转置矩阵
void zhuanzhi()
{
int i,j;
for (i=1;i<=GS;i++)
{
for(j=1;j<BM;j++)
{
BT[j][i]=B[i][j];
}
}
}
//计算权阵
void cal_quanzhen()
{
double t;
int i;
for (i=1;i<=GS;i++)
{
t=gaochaguancezhi[i].changdu;
P[i][i]=1.0/t;
}
}
//矩阵乘法
void juzhenxiangcheng()
{
int i,j,k;
double t;
//计算BTP
for(i=1;i<BM;i++)
{
for(j=1;j<=GS;j++)
{
t=0;
for(k=1;k<=GS;k++)
{
t+=BT[i][k]*P[k][j];
}
BTP[i][j]=t;
}
}
//计算BTPB
for(i=1;i<BM;i++)
{
for(j=1;j<BM;j++)
{
t=0;
for(k=1;k<=GS;k++)
{
t+=BTP[i][k]*B[k][j];
}
BTPB[i][j]=t;
}
}
//计算BTPL
for(i=1;i<BM;i++)
{
t=0;
for(k=1;k<=GS;k++)
{
t+=BTP[i][k]*l[k];
}
BTPL[i]=t;
}
}
//交换两个数的值
void swap(double *a,double *b)
{
double c;
c=*a;
*a=*b;
*b=c;
}
//BTPB矩阵求逆
int inverse()
{
int i,j,k,aa[M],bb[M];
double d;
for (i=0;i<=M;i++)
{
for (j=0;j<=M;j++)
{
BTPBN[i][j]=BTPB[i][j];
}
}
for (k=1;k
{
d=0;
for (i=k;i
{
for (j=k;j
{
if (fabs(BTPBN[i][j])>d)
{
d=fabs(BTPBN[i][j]);
aa[k]=i;
bb[k]=j;
}
}
}
if (d+1.0==1.0)
{
return 0;
}
if (aa[k]!=k)
{
for (j=1;j
{
swap(&BTPBN[k][j],&BTPBN[aa[k]][j]);
}
}
if (bb[k]!=k)
{
for (i=1;i
{
swap(&BTPBN[i][k],&BTPBN[i][bb[k]]);
}
}
BTPBN[k][k]=1/BTPBN[k][k];
for (j=1;j
{
if (j!=k)
{
BTPBN[k][j]=BTPBN[k][j]*BTPBN[k][k];
}
}
for (i=1;i
{
if (i!=k)
{
for (j=1;j
{
if (j!=k)
{
BTPBN[i][j]=BTPBN[i][j]-BTPBN[i][k]*BTPBN[k][j];
}
}
}
}
for (i=1;i
{
if (i!=k)
{
BTPBN[i][k]=-BTPBN[i][k]*BTPBN[k][k];
}
}
}
for (k=BM-1;k>=1;k--)
{
for (j=1;j<BM;j++)
{
if (bb[k]!=k)
{
swap(&BTPBN[k][j],&BTPBN[bb[k]][j]);
}
}
for (i=1;i<BM;i++)
{
if (aa[k]!=k)
{
swap(&BTPBN[i][k],&BTPBN[i][aa[k]]);
}
}
}
return 1;
}
//解法方程
void jiefafangcheng()
{
int i,k;
double t;
for(i=1;i<BM;i++)
{
t=0;
for(k=1;k<=GS;k++)
{
t+=BTPBN[i][k]*BTPL[k];
}
x[i]=t;
}
}
// 计算平差值
void cal_gaochengpingchazhi()
{
int i;
for(i=1;i<=BM;i++)
{
H[i]=JH[i]+x[i]/1000;
}
}
//计算观测值改正数
void guancezhigaizhengshu()
{
int i,j;
double t;
for(i=1;i<=GS;i++)
{
t=0;
for(j=1;j<=BM;j++) //
{
t+=B[i][j]*x[j];
}
v[i]=t-l[i];
}
}
//计算观测值平差值
void cal_guancezhipingchazhi()
{
int i;
for (i=1;i<=GS;i++)
{
h[i]=gaochaguancezhi[i].gaocha+(v[i]/1000);
}
}
// 计算观测值协因数
void cal_gczxieyinshu()
{
int i,j,k;
double BNBBN[M][M],t;
//计算BNBBN
for (i=1;i<=GS;i++)
{
for (j=1;j<BM;j++)
{
t=0;
for (k=1;k<BM;k++)
{
t+=B[i][k]*BTPBN[k][j];
}
BNBBN[i][j]=t;
}
}
//计算Qhh
for (i=1;i<=GS;i++)
{
for (j=1;j<=GS;j++)
{
t=0;
for (k=1;k<BM;k++)
{
t+=BNBBN[i][k]*BT[k][j];
}
Qhh[i][j]=t;
}
}
}
//精度评定
void jingdupingding()
{
int i;
double t=0,VTP[M];
//计算VTP
for(i=1;i<=GS;i++)
{
VTP[i]=P[i][i]*v[i];
}
//计算单位权中误差
for(i=1;i<=GS;i++)
{
t+=VTP[i]*v[i];
}
danweiquanzhongwucha=sqrt(t/(GS-BM+1));
//计算点位高程中误差
for(i=1;i<=BM;i++)
{
H_zhongwucha[i]=sqrt(BTPBN[i][i])*danweiquanzhongwucha;
}
//计算观测值平差值中误差
for (i=1;i<=GS;i++)
{
gcz_zhongwucha[i]=sqrt(Qhh[i][i])*danweiquanzhongwucha;
}
}
//数据输出
void shujushuchu()
{
int i,j;
FILE *fp;
fp=fopen("水准网数据处理成果.txt","w");
if (fp==NULL)
{
printf("打开文件失败!\n");
exit(0);
}
//输出观测数据
fprintf(fp,"原始观测数据为：\n");
fprintf(fp,"水准点1 水准点2 高差观测值(m) 两水准点间距离(km)\n");
for (i=1;i<=GS;i++)
{
fprintf(fp,"%10d%10d%20.5lfm%20.5lfkm\n",gaochaguancezhi[i].houshidian,gaochaguancezhi[i].qianshidian,gaochaguancezhi[i].gaocha,gaochaguancezhi[i].changdu);
}
//输出各点近似高程
fprintf(fp,"\n各点近似高程为：\n");
fprintf(fp,"点号 近似高程（m）\n");
for (i=1;i<=BM;i++)
{
fprintf(fp,"%5d%15.5lfm\n",i,JH[i]);
}
//输出系数矩阵（B）
fprintf(fp,"\n系数矩阵（B）为：\n");
for (i=1;i<BM;i++)
{
fprintf(fp,"x%-19d",i);
}
fprintf(fp,"\n");
for (i=1;i<=GS;i++)
{
for (j=1;j<BM;j++)
{
fprintf(fp,"%-15.5lf",B[i][j]);
}
fprintf(fp,"\n");
}
fprintf(fp,"\n");
// 输出常数项
fprintf(fp,"常数矩阵（l）为：\n");
for (i=1;i<=GS;i++)
{
fprintf(fp,"l%d=%-10.6lfmm\n",i,l[i]);
}
//输出BTPB
fprintf(fp,"\n法方程系数阵（BTPB）为:\n");
for (i=1;i<BM;i++)
{
for(j=1;j<BM;j++)
{
fprintf(fp,"%-15.5lf ",BTPB[i][j]);
}
fprintf(fp,"\n");
}
fprintf(fp,"\n");
//输出法方程常数向量（BTPL）
fprintf(fp,"法方程常数向量（BTPL）为:\n");
for (i=1;i<BM;i++)
{
fprintf(fp,"%15.5lfmm\n",BTPL[i]);
}
fprintf(fp,"\n");
//输出BTPBN
fprintf(fp,"BTPBN:\n");
for (i=1;i<BM;i++)
{
for(j=1;j<BM;j++)
{
fprintf(fp,"%-15.6lf ",BTPBN[i][j]);
}
fprintf(fp,"\n");
}
fprintf(fp,"\n");
//输出Qhh
fprintf(fp,"Qhh:\n");
for (i=1;i<=GS;i++)
{
for(j=1;j<=GS;j++)
{
fprintf(fp,"%-15.6lf ",Qhh[i][j]);
}
fprintf(fp,"\n");
}
//输出权阵
fprintf(fp,"\n观测值的权阵为：\n");
for (i=1;i<=GS;i++)
{
for(j=1;j<=GS;j++)
{
fprintf(fp,"%-15.6lf ",P[i][j]);
}
fprintf(fp,"\n");
}
//输出法方程的解
fprintf(fp,"\n法方程参数解为:\n");
for (i=1;i<BM;i++)
{
fprintf(fp,"x%d=%-10.5lfmm\n",i,x[i]);
}
fprintf(fp,"\n");
//输出观测值改正数的解
fprintf(fp,"观测值改正数的解:\n");
for (i=1;i<=GS;i++)
{
fprintf(fp,"v%d=%lfmm\n",i,v[i]);
}
fprintf(fp,"\n");
//输出观测值平差值
fprintf(fp,"观测值平差值:\n");
for (i=1;i<=GS;i++)
{
fprintf(fp,"h%d=%lfm\n",i,h[i]);
}
fprintf(fp,"\n");
//输出单位权中误差
fprintf(fp,"单位权中误差为：\n");
fprintf(fp,"σ0=%lfmm\n",danweiquanzhongwucha);
fprintf(fp,"\n");
//输出观测值协因数
fprintf(fp,"观测值平差值的协因数为:\n");
for (i=1;i<=GS;i++)
{
fprintf(fp,"Qh%dh%d=%lf\n",i,i,Qhh[i][i]);
}
fprintf(fp,"\n");
//输出观测值平差值中误差
fprintf(fp,"观测值平差值的中误差为:\n");
for (i=1;i<=GS;i++)
{
fprintf(fp,"σh%d=%lfmm\n",i,gcz_zhongwucha[i]);
}
fprintf(fp,"\n");
//输出高程平差值中误差
fprintf(fp,"各点高程平差值的中误差为:\n");
for (i=1;i<=BM;i++)
{
fprintf(fp,"σH%d=%lfmm\n",i,H_zhongwucha[i]);
}
fprintf(fp,"\n");
//输出最终成果
fprintf(fp,"最终成果\n观测值平差值及精度\n");
fprintf(fp,"起止点 高差观测值平差值 高差观测值平差值中误差\n");
for (i=1;i<=GS;i++)
{
fprintf(fp,"%5d%5d%20.5lfm%25.5lfmm\n",gaochaguancezhi[i].houshidian,gaochaguancezhi[i].qianshidian,h[i],gcz_zhongwucha[i]);
}
fprintf(fp,"\n各点高程及精度\n");
fprintf(fp,"点号 高程(m) 高程中误差(mm)\n");
for (i=1;i<=BM;i++)
{
fprintf(fp,"%5d%15.5lfm%15.5lfmm\n",i,H[i],H_zhongwucha[i]);
}
fclose(fp);
}
void main()
{
shuzufuchushizhi(); // 给数组赋初始值
shujuduqu(); // 读取数据
jisuangaochengjinsizhi(); // 计算高程近似值
cal_xishu(); // 计算系数矩阵
cal_l(); // 计算常数项
zhuanzhi(); // 计算BT
cal_quanzhen(); // 计算权阵
juzhenxiangcheng(); // 矩阵相乘，计算BTP、BTPB、BTPL
inverse(); // 计算BTPBN
jiefafangcheng(); // 解法方程
cal_gaochengpingchazhi(); // 计算高程平差值
guancezhigaizhengshu(); // 计算观测值改正数
cal_guancezhipingchazhi(); // 计算观测值平差值
cal_gczxieyinshu(); // 计算观测值协因数
jingdupingding(); // 精度评定
shujushuchu(); // 数据输出
}
这是水准网平差原始数据.txt文件，老是无法读取文件，无法输出正确的处理后的数值。每次执行只显示“press any key to continue”
