clc;
clear;
disp（‘*********************************************測繪工程16-2楊帆*****************************************************‘）
%讀取txt文本
[filename1filepath]=uigetfile（‘*.txt‘‘請選擇GPS數據文件‘）;%uigetfile顯示檢索文件的對話框，filename1為文件名
yf=fopen（strcat（filepathfilename1）‘rt‘）;%strcat:橫向連接字符串
if?（yf==-1）
???msgbox（‘出錯！出錯！出錯！‘）;?%magbox為彈出警告的指令
???return;
end

yzds=fscanf（yf‘%f‘1）;%已知點數
if?yzds==0
????x01=fscanf（yf‘%f‘1）;%已知點x坐標
????y01=fscanf（yf‘%f‘1）;%已知點y坐標
????z01=fscanf（yf‘%f‘1）;%已知點z坐標
????wzds=fscanf（yf‘%f‘1）;%未知點數
????jxs=fscanf（yf‘%f‘1）;%基線數
????sigma0=fscanf（yf‘%f‘1）;%中誤差
????i=1;
??while?（i????qd（i）=fscanf（yf‘%f‘1）;%起點
????zd（i）=fscanf（yf‘%f‘1）;%終點
????dx（i）=fscanf（yf‘%f‘1）;%△x
????dy（i）=fscanf（yf‘%f‘1）;%△y
????dz（i）=fscanf（yf‘%f‘1）;%△z
????sigmaxx（i）=fscanf（yf‘%f‘1）;%σxx
????sigmayy（i）=fscanf（yf‘%f‘1）;%σyy
????sigmazz（i）=fscanf（yf‘%f‘1）;%σzz
????sigmayx（i）=fscanf（yf‘%f‘1）;%σyx
????sigmazy（i）=fscanf（yf‘%f‘1）;%σzy
????sigmazx（i）=fscanf（yf‘%f‘1）;%σzx
????i=i+1;
??end
disp（‘**************************************************************************************************************‘）
disp（‘????????????????????????????????????????????GPS無約束網平差結果‘）
??%權陣P的形成
??disp（‘**********************************************************P陣**********************************************************‘）
??T=zeros（jxs*3jxs*3）;%建立行數列數為基線向量條數3倍的零矩陣，用來存儲方差協方差陣
??for?i=1:jxs?%存協方差
????????T（3*i-13*i-2）=sigmayx（i）;%σyx
????????T（3*i3*i-1）=sigmazy（i）;%σzy
????????T（3*i3*i-2）=sigmazx（i）;%σzx
??end
????T;
????D=T+T‘;%方差陣
????disp（D）;
??for?i=1:jxs?%存方差
????????D（3*i-23*i-2）=sigmaxx（i）;%σxx
????????D（3*i-13*i-1）=sigmayy（i）;%σyy
????????D（3*i3*i）=sigmazz（i）;%σzz
??end
????D;
????P=inv（（D/sigma0^2））;%p=inv（Q）q=D/σ0^2
????disp（P）;
????%系數矩陣B的形成
????disp（‘**********************************************************B陣**********************************************************‘）
????B=zeros（jxs*3wzds*3）;%nxt
??for?i=1:jxs
????????B（3*i-2:3*i3*qd（i）-2:3*qd（i））=diag（[-1-1-1]）;%對角陣
??end
??for?i=1:jxs
????????B（3*i-2:3*i3*zd（i）-2:3*zd（i））=diag（[111]）;
??end????
????disp（B）;
??%常數項矩陣l的形成
??disp（‘**********************************************************l陣**********************************************************‘）
??l=zeros（jxs*31）;
??for?i=1:jxs
??????l（3*i-2:3*i）=[dx（i）;dy（i）;dz（i）];
??end
???disp（l）;
?%間接平差
????NBB=B‘*P*B+G0*G;
????w=B‘*P*l;
????xp=inv（NBB）*w;
????v=B*xp-l;%誤差方程
????delta0=sqrt（（v‘*P*v）/（wzds-1））;%σ0
????disp（‘-------------------------------------------單位權中誤差（單位：m）--------------------------------------------‘）
????disp（strcat（‘σ0=‘num2str（delta0）））?%num2str為把數值轉換成字符串，轉換后可以使用disp函數輸出
????Q0=inv（NBB）;%σBbb-1
????Q=Q0*B‘*P*B*Q0;
????%由假定已知點計算其他未知點平差后坐標
????jxxl=v+l;
??for?i=1:wzds-1
????????X（1）=x01;Y（1）=y01;Z（1）=z01;
????????X（i+1）=X（i）-jxxl（3*i-2）;
????????Y（i+1）=Y（i）-jxxl（3*i-1）;
????????Z（i+1）=Z（i）-jxxl（3*i）;
??end
????dh=ze