function leveling()
clc;
clear;
% a： 已知点个数 b：未知点个数 c：观测值个数 D：点号
% HH：已知点高程 h：观测值高差 s：观测值距离
% Q： 起点点号 Z：终点点号
disp('观测点号全部用数字表示！');
disp('数据格式：');
disp('a b c (已知点个数 未知点个数 观测值个数)');
disp('1 12.123 (已知点点号 已知高程 共a行)');
disp('1 3 0.123 2.5 (观测值起点 终点 观测高差 观测距离 共c行)')
disp('%%%%%%%%%%%%%%%%');
%% 读入原始数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D=[]; HH=[]; h=[]; s=[]; Q=[]; Z=[];
[filename,filepath]=uigetfile('*.txt','选择平差文件:');
name=[filepath filename];
fid1=fopen(name,'rt');
if(fid1==-1)
disp('文件未打开，请重试！');
return;
end
a=fscanf(fid1,'%f',1); %读入已知点个数
b=fscanf(fid1,'%f',1); %读入未知点个数
ab=a+b; %总点数ab
c=fscanf(fid1,'%f',1); %读入观测值个数
for(i=1:a)
D(i)=fscanf(fid1,'%f',1); %点号
HH(i)=fscanf(fid1,'%f',1); %高程值
end
for(j=1:c)
Q(j)=fscanf(fid1,'%f',1); %起点
Z(j)=fscanf(fid1,'%f',1); %终点
h(j)=fscanf(fid1,'%f',1); %高程值
s(j)=fscanf(fid1,'%f',1); %距离值
end
fclose(fid1);
%查看变量
%% 不用输入所有点号 自动寻找所有点号%%%%%%%%%%%%%%%%%%%%%%%%%%%
jj=1; %把未知点点号存到D矩阵中！
for(i=1:c)
ii=0;
for(j=1:a)
if(Q(i)==D(j))
break;
else
ii=ii+1;
end
if(ii==a)
D(a+jj)=Q(i);
a=a+jj;
break;
end
end
end
for(i=1:c)
ii=0;
for(j=1:a)
if(Z(i)==D(j))
break;
else
ii=ii+1;
end
if(ii==a)
D(a+jj)=Z(i);
a=a+jj;
break;
end
end
end
%% 转变点号%%%%%%%%%%%%%%%%%
i=length(D);
pn=D;
[be,k1]=changepiont(ab,pn,Q);
[en,k2]=changepiont(ab,pn,Z);
%% 生成误差方程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
HHH=HH;
for(i=1:b)
HHH(ab-b+i)=0; % 给未知高程点添加高程零生成HHH
end
%% 生成误差方程常数项l %%
for(i=1:c)
e=be(i);
f=en(i);
l(i)=h(i)-(HHH(f)-HHH(e));
end
%%ll=zeros(c,1); % 解常数项(书上代码)
%%ll=HH(be)-HH(en)+h; % 解常数项(书上代码)
l=l';
%% 生成系数阵B %%
B=zeros(c,b);
for(i=1:c)
e=be(i);
f=en(i);
i1=e-(ab-b);
i2=f-(ab-b);
if(HHH(e)==0)
B(i,i1)=-1;
end
if(HHH(f)==0)
B(i,i2)=1;
end
end
%% 求解 观测值改正数V 高程改正数x 并计算高程HH %%
p=diag(1./s); %解权阵
x=inv(B'*p*B)*(B'*p*l); % 高程改正数x
for(i=1:b)
HH(i+ab-b)=x(i);
end
v=B*x-l; %观测值改正数
%HH %高程成果
L=h'+v; %高差观测值
%L %高差改正值
%% 计算中误差%%%%%%%%%
u0=sqrt(v'*p*v/(c-b));
Qxx=inv(B'*p*B);
u1=u0*sqrt(diag(Qxx)); %待定点高程平差值中误差
Qll=B*Qxx*B';
uu=u0*sqrt(diag(Qll));
for(i=1:ab-b)
x2(i)=0;
u2(i)=0;
end
for(i=1:b)
x2(i+ab-b)=x(i);
u2(i+ab-b)=u1(i);
end
%% 输出水准平差结果%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('水准平差结果如下：');
disp('具体见工作目录result.txt');
fprintf('\n');
disp('点号：');
D %点号
disp('高程改正数(零为已知点)：'); x2 %高程改正数
disp('高程改正数(m)：'); HHH %高程改正数
disp('高程中误差(零为已知点)：'); u2 %高程中误差
disp('高程平差值(m)：'); HH %高程平差值
fid2=fopen('result.txt','wt'); %% wt不同于w或是w+，\n不能识别
if(fid2==-1)
disp('创建result文件未成功，请重试！');
return;
end
fprintf(fid2,'1.水准网观测值平差值:\n 起点 终点 观测高差 高差改正数 平差高差(m)\n');
for(i=1:ab)
fprintf(fid2,'%5d %5d %10.4f %10.4f %10.4f\n',Q(i),Z(i),h(i),v(i),L(i));
end
fprintf(fid2,'2.水准网平差结果: \n 点号 原始高程 高程改正数 高程中误差 高程(m)\n');
for(i=1:ab)
fprintf(fid2,'%5d %10.4f %10.4f %10.4f %10.4f\n',D(i),HHH(i),x2(i),u2(i),HH(i));
end
fprintf(fid2,' 注：高程改正数为零即是 已知点!\n');
fprintf(fid2,'3.单位权中误差u0= %5.4fm\n',u0);
fclose(fid2);
%%%%%%%%%%%%%%%%%%%%%%%%% 水准平差完成 王建%%%%%%%%%%%%%%%%%
