基于MATLAB的控制网平差程序设计--第六章源代码

近似坐标计算的函数-calcux0y0函数（126页）

function [x0,y0]=calcux0y0(x0,y0,e,d,sid,g,f,dir,s,t,az,pn,xyknow,xyunknow,point,aa,bb,cc)

%本函数的作用是计算待定点的近似坐标

format short; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% time=0;prelength=length(xyknow);non_orient=0;point_angle=0;

while length(xyunknow)>0

%考虑的计算方法有：1.极坐标;2.前方交会；3.测边交会；4.后方交会；

%5.无定向导线的两种情况：（1）已知两个点；（2）分离的已知点与方位角；基本思路：%采用循环的方法逐一对每一个未知点进行以上各种方法条件的搜索，满足后即解算。

aa0=[];bb0=[];cc0=[];%记录搜索到两条观测边但需用户给顺序的点，注意要放在while 里面。

time=time+1; % 用于统计循环次数。

way=0;

for i=xyunknow %依次循环向量中的各元素

%============================================================= =====

%方法1.极坐标条件搜索与计算-->way=1,基本思路：找到或求出一个方位角，找出一条边。

temp1=[]; temp2=[]; temp3=[]; temp4=[]; temp5=[]; temp6=[]; temp7=[]; temp8=[];

temp9=[]; temp10=[];A=[];B=[];P=[];

%第一步：寻找观测条件：两种情况：一是有已知方位角；二是由两个已知点及方向观测值推出方位角。

temp7=find(t==i);

if length(temp7)>0

temp8=find(xyknow==s(temp7(1)));

if length(temp8)>0

temp9=find(e==s(temp7(1))&d==t(temp7(1))|e==t(temp7(1))&d==s(temp7(1)));

if length(temp9)>0 %第一种情况：有已知方位角（一般适用于闭合导线）

S=mean(sid(temp9));

Alfa0=az(temp7(1));

x0(i)=x0(s(temp7(1)))+S*cos(Alfa0);

y0(i)=y0(s(temp7(1)))+S*sin(Alfa0);

way=1;method(i)=11;%----->由已知方位角算出

end

end

else

temp1=find(f==i);%第二种情况：由两个已知点及方向观测值推出方位角（适用于附合、支导线）

if length(temp1)>0

for j=xyknow

temp2=find(g(temp1)==j);

if temp2>0

temp3(end+1)=temp1(temp2);

end

end %找到所有以已知点为起点，以i为终点的方向值位置，并存放于temp3中

if length(temp3)>0

for j=temp3

temp4=find(e==g(j)&d==f(j)|e==f(j)&d==g(j));

if temp4>0 %找到一条观测边

temp5=g(j);

for k=xyknow

temp6=find(g==temp5&f==k);

if temp6>0 %找到另一个方向值位置

A=k;B=temp5;

S=mean(sid(temp4));

dir1=dir(temp6);dir2=dir(j);

way=1;

method(i)=12;%由两个已知点及方向观测值算得

break;

end

end

end

end

end

end

%第二步：按极坐标公式计算坐标

if way==1

deltax=x0(B)-x0(A);

deltay=y0(B)-y0(A);

a1=alfa(deltax,deltay);

a2=a1+dir2-dir1;

if a2<0

a2=a2+2*pi;

end

if a2>2*pi

a2=a2-2*pi;

end

x0(i)=x0(B)+S*cos(a2);

y0(i)=y0(B)+S*sin(a2);

end

end

%%================================================================ ======

%方法2：前方搜索与计算way=2

%思路：找出合适的A,B点，然后对A,B,P进行逆时针排序，照公式计算坐标

if way==0

%第一步：搜索合适的A,B点

temp1=[];temp2=[];temp3=[];temp4=[];temp5=[];temp6=[];temp7=[];temp8=[];

A=[];B=[];P=[];

temp1=find(f==i); %找出终点为i的方向值位置

for j=xyknow

temp2=find(g(temp1)==j);

if temp2>0

temp3(end+1)=temp1(temp2);

end

end %找出起点为已知点，终点为目标点未知点的方向值位置，记录在temp3中

if length(temp3)>1

for j=g(temp3)

for k=g(temp3)

temp4=find(g==j&f==k);

if temp4>0

temp5(end+1)=temp4;%temp5中记录了所有与目标未知点有观测的已知点之间方向观测值的位置

end

end

if length(temp5)>0

for j=temp5

temp6=find(g(temp5)==f(j)&f(temp5)==g(j));

if temp6>0

way=2;method(i)=2;

P=i;

A=g(i);

B=f(j);

break;

end

end

end %找出一组合适的A,B即可，退出循环

end

if way==2 %第二步：对A,B,P进行逆时针排序

[A,B,P]=reorderunknow(A,B,P,1,g,f,dir);

bet1=beta(P,A,B,g,f,dir);

bet2=beta(P,B,A,g,f,dir);

%%第三步：代入公式计算坐标

if(bet1~=pi/2)&(bet2~=pi/2)

x0(P)=(x0(A)*tan(bet1)+x0(B)*tan(bet2)+(y0(B)-y0(A))*tan(bet1)*tan(bet2))/(tan(bet1)+tan(bet2 ));