近似坐标计算的函数-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)>0temp8=find(xyknow==s(temp7(1)));if length(temp8)>0temp9=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;%----->由已知方位角算出endendelsetemp1=find(f==i);%第二种情况：由两个已知点及方向观测值推出方位角（适用于附合、支导线）if length(temp1)>0for j=xyknowtemp2=find(g(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);endend %找到所有以已知点为起点，以i为终点的方向值位置，并存放于temp3中if length(temp3)>0for j=temp3temp4=find(e==g(j)&d==f(j)|e==f(j)&d==g(j));if temp4>0 %找到一条观测边temp5=g(j);for k=xyknowtemp6=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;endendendendendend%第二步：按极坐标公式计算坐标if way==1deltax=x0(B)-x0(A);deltay=y0(B)-y0(A);a1=alfa(deltax,deltay);a2=a1+dir2-dir1;if a2<0a2=a2+2*pi;endif a2>2*pia2=a2-2*pi;endx0(i)=x0(B)+S*cos(a2);y0(i)=y0(B)+S*sin(a2);endend%%================================================================ ======%方法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=xyknowtemp2=find(g(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);endend %找出起点为已知点，终点为目标点未知点的方向值位置，记录在temp3中if length(temp3)>1for j=g(temp3)for k=g(temp3)temp4=find(g==j&f==k);if temp4>0temp5(end+1)=temp4;%temp5中记录了所有与目标未知点有观测的已知点之间方向观测值的位置endendif length(temp5)>0for j=temp5temp6=find(g(temp5)==f(j)&f(temp5)==g(j));if temp6>0way=2;method(i)=2;P=i;A=g(i);B=f(j);break;endendend %找出一组合适的A,B即可，退出循环endif 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 ));y0(P)=(y0(A)*tan(bet1)+y0(B)*tan(bet2)+(x0(A)-x0(B))*tan(bet1)*tan(bet2))/(tan(bet1)+tan(bet2 ));elsesidab=sqrt((x0(A)-x0(B))^2+(y0(A)-y0(B))^2);alfaab=alfa(x0(B)-x0(A),y0(B)-y0(A));alfaap=alfaab-bet1;if alfaap<0alfaap=alfaap+2*pi;endif bet1==pi/2sidap=tan(bet2)*sidab;elsesidap=sidab/cos(bet1);endx0(P)=x0(A)+sidap*cos(alfaap);y0(P)=y0(A)+sidap*sin(alfaap);endend %for if wey==2end % for if way==0%%============================================================ ======%方法3 ：测边交会与计算way=3,思路：三种算法：一是找到一个端点为所求目标未知点的三条边%二是寻求方向观测信息以对A,B,P进行逆时针排序，三是给出A,B,P的排列顺序。

if way==0temp1=[];temp2=[];temp3=[];temp4=[];temp5=[];temp6=[];temp7=[];temp8=[];i1=0;A=[];B=[];C=[];P=[];%第一步：搜索合适的A,B,C点或带有方向观测值的A，B点temp1=find(e==i|d==i);if length(temp1)>0for j=xyknowtemp2=find(e(temp1)==j|d(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);endendif length(temp3)>2 %第一种: 找到至少带有已知点的3条边，不需排序即可进行计算way=3;i1=1;P=i;if e(temp3(1))==iA=d(temp3(1)); %AelseA=e(temp3(1));endif e(temp3(2))==iB=d(temp3(2)); %BelseB=e(temp3(2));endif e(temp3(3))==iC=d(temp3(3)); %CelseC=e(temp3(3));ends1=sid(temp3(1));s2=sid(temp3(2));s3=sid(temp3(3));elseif length(temp3)==2 %只能找到两条带有已知点的边P=i;if e(temp3(1))==iA=d(temp3(1)); %AelseA=e(temp3(1));endif e(temp3(2))==iB=d(temp3(2)); %BelseB=e(temp3(2));end%第二种：需进一步查找方向观测值信息，以对A,B,P进行逆时针排序n1=find(g==A&f==P);n2=find(g==A&f==B);n3=find(g==B&f==P);n4=find(g==B&f==A);n5=find(g==P&f==A);n6=find(g==P&f==B);if n1>0&n2>0way=3;i1=2;[A,B,P]=reorderunknow(A,B,P,1,g,f,dir);elseif n3>0&n4>0way=3;i1=2;[B,A,P]=reorderunknow(B,A,P,1,g,f,dir);elseif n5>0&n6>0way=3;i1=2;[A,B,P]=reorderunknow(A,B,P,2,g,f,dir);else %第三种：需给出A,B,P的排列顺序temp4=find(cc==P);if length(temp4)>0A=aa(temp4(1));B=bb(temp4(1));way=3;i1=2;elseaa0(end+1)=A;bb0(end+1)=B,cc0(end+1)=P;endendendend%第二步：坐标计算if way==3&i1==1 %三个已知点，三条边长method(i)=31;deltx=x0(B)-x0(A);delty=y0(B)-y0(A);alfaAB=alfa(deltx,delty);ss=sqrt(deltx*deltx+delty*delty);ee=(s1*s1+ss*ss-s2*s2)/(2*ss);ff=sqrt(s1*s1-ee*ee);xp1=x0(A)+ee*cos(alfaAB)+ff*sin(alfaAB);yp1=y0(A)+ee*sin(alfaAB)-ff*cos(alfaAB);xp2=x0(A)+ee*cos(alfaAB)-ff*sin(alfaAB);yp2=y0(A)+ee*sin(alfaAB)+ff*cos(alfaAB);s4=sqrt((xp1-x0(C))^2+(yp1-y0(C))^2);s5=sqrt((xp2-x0(C))^2+(yp2-y0(C))^2);if abs(s4-s3)<abs(s5-s3)x0(P)=xp1;y0(P)=yp1;elsex0(P)=xp2;y0(P)=yp2;endelseif way==3&i1==2 %由方向观测值确定逆时针顺序后坐标计算q1=find(e==A&d==P|e==P&d==A);s1=sid(q1);q2=find(e==B&d==P|e==P&d==B);s2=sid(q2);method(i)=32;deltx=x0(B)-x0(A);delty=y0(B)-y0(A);alfaAB=alfa(deltx,delty);ss=sqrt(deltx*deltx+delty*delty);ee=(s1*s1+ss*ss-s2*s2)/(2*ss);ff=sqrt(s1*s1-ee*ee);x0(P)=x0(A)+ee*cos(alfaAB)+ff*sin(alfaAB);y0(P)=y0(A)+ee*sin(alfaAB)-ff*cos(alfaAB);endend %if way==0%=============================================================== =====%方法4：后方交会搜索与计算way=4%思路：1.找到三个已知点；2.判断是否是危险圆；3.坐标计算%第一步：找点if way==0temp1=[];temp2=[];temp3=[];temp4=[];temp5=[];temp6=[];A=[];B=[];C=[];P=[];temp1=find(g==i);%找出起点为P的方向值位置if length(temp1)>2for j=xyknowtemp2=find(f(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);%记录起点为P终点为已知点的方向值位置endendif length(temp3)>2 %找到三个或者三个以上的点A=f(temp3(1));B=f(temp3(2));P=i;for j=3:length(temp3)C=f(temp3(j));[A,B,C]=reorderknow(A,B,C,x0,y0);%将A,B,C排成顺时针方向sa=sqrt((x0(C)-x0(B))^2+(y0(C)-y0(B))^2);sb=sqrt((x0(C)-x0(A))^2+(y0(C)-y0(A))^2);sc=sqrt((x0(A)-x0(B))^2+(y0(A)-y0(B))^2);angleA=acos((sb*sb+sc*sc-sa*sa)/(2*sb*sc));angleB=acos((sa*sa+sc*sc-sb*sb)/(2*sa*sc));angleC=acos((sa*sa+sb*sb-sc*sc)/(2*sa*sb));temp4=find(g==i&f==A);RA=dir(temp4);temp5=find(g==i&f==B);RB=dir(temp5);temp6=find(g==i&f==C);RC=dir(temp6);alfa1=RC-RB;beta1=RA-RC;gama1=RB-RA;alfa2=beta(B,P,C,g,f,dir);beta2=beta(A,P,C,g,f,dir);gama2=beta(A,P,B,g,f,dir);%%第二步：判断四点是否共圆if(angleA==alfa2&angleB==beta2)|(angleA==alfa2&angleC==gama2)|(angleC==gama2&angleB ==beta2)continue %危险圆！！！！！！！！！！！！！！！！else%第三步：非危险圆，坐标计算way=4;method(i)=4;PA=tan(alfa1)*tan(angleA)/(tan(alfa1)-tan(angleA));PB=tan(beta1)*tan(angleB)/(tan(beta1)-tan(angleB));PC=tan(gama1)*tan(angleC)/(tan(gama1)-tan(angleC));x0(P)=(PA*x0(A)+PB*x0(B)+PC*x0(C))/(PA+PB+PC);y0(P)=(PA*y0(A)+PB*y0(B)+PC*y0(C))/(PA+PB+PC);break;endend % for j=3:length(temp3)end %if length(temp3)>2end %if length(temp1)>2end %if way==0%=============================================================== =====%求出一个未知点的近似坐标后终止循环，更新已知点点号和未知点点号数组if way>0temp=[];xyknow(end+1)=i;temp=find(xyunknow==i);xyunknow(temp:end-1)=xyunknow(temp+1:end);xyunknow=xyunknow(1:end-1);break; %%%%%%% for i=xyunknowend % if way>0end %for i=xyunknow%================================================================= =====%方法5.倘若对所有未知点循环一次的过程中一个也未求出，则考虑用无定向导线计算%way==5%第一种情况：已知条件为两个或者两个以上点的坐标if(length(xyknow)==prelength)&(way==0)&(length(xyknow)>=2)siddir=0;seeksid=0;seekdir=0;for j=xyknowfor k=xyunknowtemp1=[];temp1=find(e==j&d==k|e==k&d==j);if(length(temp1)>0) %找到与已知点联测的一条边的另一端点seeksid=1;siddir=1;break;endendif seeksid==1break;endendif seeksid==0 %未找到与已知点联测的一条边的另一端点，则找方向观测的另一个端点for j=xyknowfor k=xyunknowtemp1=[];temp1=find(g==j&f==k|g==k&f==j);if(length(temp1)>0) %找到方向观测的另一端点seekdir=1;seekdir=1;break;endendif seekdir==1break;endendendif siddir==1 % (length(temp1)>0)trueX0=x0;trueY0=y0;x0=[];y0=[];xyknow=[];xyunknow=[];x0(j)=0.0;y0(j)=0.0; % 能用此方法解决，建立新的坐标系if seeksid==1x0(k)=mean(sid(temp1));y0(k)=0;elseif seekdir==1x0(k)=1000;y0(k)=0;endnon_orient=1; %更新已知点数组和未知点数组xyknow=[j k];point(find(point==k|point==j))=[];xyunknow=point';end%============================================================== =====%第二种情况：已知条件为一个已知点和一个方位角，但分散开。