基于MATLAB的控制网平差程序设计--第五章源代码精品资料.doc

《基于MATLAB的控制网平差程序设计--第五章源代码精品资料.doc》由会员分享，可在线阅读，更多相关《基于MATLAB的控制网平差程序设计--第五章源代码精品资料.doc（24页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、观测数据读入程序，rddat1函数（85页）global net ed dd sd dd1 pn x0 y0 m1 m2 m3 ms pp e d sid md g f dir ni si ma s t az aa bb cc rt rr ttglobal pathname filenamex0=;y0=;e=;d=;sid=;g=;f=;dir=;si=;ni=;s=;t=;az=;pn=;filename,pathname=uigetfile(*.txt,请选择原始数据);fit1=fopen(strcat(pathname,filename),rt);if(fit1=-1) msgbox

2、(Input File or Path is not correct,Warning,warn); return;endnet=fscanf(fit1,%d,1);a=fscanf(fit1,%d,3);ed=a(1);dd=a(2);dd1=a(3);sd=ed+dd;pn=fscanf(fit1,%d,sd);a=fscanf(fit1,%f,2*ed);for i=1:ed x0(i)=a(2*i-1); y0(i)=a(2*i);enda=fscanf(fit1,%d,3);m1=a(1);m2=a(2);m3=a(3);isid=0; a=fscanf(fit1,%f,2); ms=

3、a(1); pp=a(2); a=fscanf(fit1,%d %d %f,3*m1); for i=1:m1 e(i)=a(3*i-2);d(i)=a(3*i-1);sid(i)=a(3*i); end e,i1=chkdat(sd,pn,e); d,i2=chkdat(sd,pn,d); i3=0; isid=i1+i2+i3;idir=0; md=fscanf(fit1,%f,1); a=fscanf(fit1,%d %d %f,3*m2); for i=1:m2 n1(i)=a(3*i-2); n2(i)=a(3*i-1); unk(i)=a(3*i); end n1,i1=chkda

4、t(sd,pn,n1); n2,i2=chkdat(sd,pn,n2); i3=0; ik=1; si(1)=1; for i=1:sd ii=0; for j=1:m2 if(n1(j)=j) ii=ii+1; g(ik)=n1(j); f(ik)=n2(j); dir(ik)=unk(j); ik=ik+1; end end ni(i)=ii; si(i+1)=si(i)+ni(i); end idir=i1+i2+i3;iaz=0;if(m30) ma=fscanf(fit1,%f,1); a=fscanf(fit1,%d %d %f,3*m3); for i=1:m3 s(i)=a(3

5、*i-2); t(i)=a(3*i-1); az(i)=a(3*i); end s,i1=chkdat(sd,pn,s); t,i2=chkdat(sd,pn,t); i3=0; iaz=i1+i2+i3;endkk=isid+idir+iaz;if(kk0) msgbox(Error by function rddat1,Warning,warn); return;endfclose(all);误差方程与法方程的组成函数-obnorm函数（90页） function obnormglobal ed dd dd1 ni si e d g f s t global m1 m2 m3 ms pp

6、md ma x0 y0 sid dir az c fit1 fit2global a q1 pa3 qls wlo=2062.648062470964;m=m1+m2+m3;n=2*dd;sum=n*(n+1)/2.0;sd=ed+dd;a(1:m,1:9)=0.0;for i=1:sd ii=4*(ni(i)+1); pa3(i,1:ii)=0.0;endc(1:sum)=0.0;w(1:n)=0.0;for i=1:m1 %边长观测误差方程 dx=x0(d(i)-x0(e(i); dy=y0(d(i)-y0(e(i); ss=sqrt(dx*dx+dy*dy); cosa=dx/ss; s

7、ina=dy/ss; a(i,1)=2*e(i)-1-2*ed+1.0e-9; a(i,2)=-cosa; a(i,3)=a(i,1)+1; a(i,4)=-sina; a(i,5)=2*d(i)-1-2*ed+1.0e-9; a(i,6)=cosa; a(i,7)=a(i,5)+1; a(i,8)=sina; a(i,9)=100.0*(ss-sid(i); q1(i)=(ms2+(ss*pp*0.0001)2);endq1(m1+1:m2+m1)=md*md;for i=1:sd if(ni(i)=0) continue; end jj=5; z0=0.0; zal=0; for j=si

8、(i):s(i)+ni(i)-1 dx=x0(f(j)-x0(g(j); dy=y0(f(j)-y0(g(j); a0=alfa(dx,dy); z1=a0-dir(j); if(z1=2.0*pi) ss=ss-2.0*pi; end a(ii,9)=(a0-ss)*lo*100.0; pa3(i,jj)=a(ii,5); pa3(i,jj+1)=a(ii,6)/q1(ii); pa3(i,jj+2)=a(ii,7); pa3(i,jj+3)=a(ii,8)/q1(ii); pa3(i,2)=pa3(i,2)+a(ii,2)/q1(ii); pa3(i,4)=pa3(i,4)+a(ii,4)

9、/q1(ii); jj=jj+4; end pa3(i,1)=a(ii,1); pa3(i,3)=a(ii,3); qls(i)=-zal;endfor i=1:m3 dx=x0(t(i)-x0(s(i); dy=y0(t(i)-y0(s(i); a0=alfa(dx,dy,a0); ss=dx*dx+dy*dy; ai=-dy/ss*lo; bi=dx/ss*lo; ii=m1+m2+i; a(ii,1)=2*s(i)-1-2*ed+1.0e-9; a(ii,2)=-ai; a(ii,3)=a(ii,1)+1; a(ii,4)=-bi; a(ii,5)=2*t(i)-1-2*ed+1.0e-

10、9; a(ii,6)=ai; a(ii,7)=a(ii,5)+1; a(ii,8)=bi; if(a0-az(i)pi) a(ii,9)=(a0-az(i)-2.0*p)*lo*100; else a(ii,9)=(a0-az(i)*lo*100.0; end q1(ii)=ma*ma;endfor i=1:m %形成法方程 for j=1:4 jj=fix(a(i,2*j-1); if(jj=0) continue; end w(jj)=w(jj)+a(i,2*j)*a(i,9)/q1(i); di=(jj-1)*(n-jj/2.0); for k=1:4 kk=fix(a(i,2*k-1)

11、; if(kkkk) continue; end c(di+kk)=c(di+kk)+a(i,2*k)*a(i,2*j)/q1(i); end endendif(m20) %和误差方程形成法方程 for i=1:sd if(ni(i)=0) continue; end for j=1:2*(ni(i)+1) jj=fix(pa3(i,2*j-1); if(jj=0) continue; end di=(jj-1)*(n-jj/2); for k=1:2*(ni(i)+1) kk=fix(pa3(i,2*k-1); if(kkkk) continue; end c(di+kk)=c(di+kk)

12、+pa3(i,2*k)*pa3(i,2*j)/qls(i); end end endendreturn 平差值与精度评定（94页）global net ed dd sd dd1 pn x0 y0 m1 m2 m3 ms pp e d sid md g f dir ni si ma s t azglobal aa bb cc rt rr tt global a q1 pa3 qls w c x y uw0global pathname net_name s_datafile a_datafile;fit2=fopen(strcat(pathname,net_name,a_datafile,wt)

13、;if(fit2=-1) msgbox(Input File or Path is not correct,Warning,warn); return;endk=1;while(k) m=m1+m2+m3; obnorm; c=invsqr(c,2*dd); uw0,k=adjxy(fit2);endellipse(uw0,fit2);n=2*dd;sum=n*(n+1)/2.0;n1=2*(ed+dd);sum1=n1*(n1+1)/2.0;for i=sum1:-1:sum1-sum+1 c(i)=c(i-(sum1-sum);endfor i=1:2*ed di=(i-1)*(n1-i/

14、2.0); for j=i:n1 if(j=i) c(di+j)=0.00000001; else c(di+j)=0.0; end endendif(m10) adjs(uw0,fit2);endif(m20) adjd(uw0,fit2);endif(m30) adja(uw0,fit2);endfclose(fit2);open(strcat(pathname,net_name,a_datafile); 坐标改正数计算及单位权中误差计算函数-adjxy函数（95页）function uw0,k=adjxy(fit2)global ed dd dd1 ni si e d g f s t p

15、n x y global m1 m2 m3 ms md ma x0 y0 sid dir az cglobal a q1 pa3 qls wsd=ed+dd;n=2*dd;k=0;for i=1:n dxy(i)=0.0; di=(i-1)*(n-i/2.0); for j=1:n dj=(j-1)*(n-j/2.0); if(j1.0) k=1; end dxy(i)=dxy(i)/100.0;endx(1:ed)=x0(1:ed);y(1:ed)=y0(1:ed);for i=1:dd x(ed+i)=x0(ed+i)+dxy(2*i-1); y(ed+i)=y0(ed+i)+dxy(2*

16、i);endx0(1:sd)=x(1:sd);y0(1:sd)=y(1:sd);for i=1:sd if(i0) ii=0; for i=1:sd if(ni(i)=0; ii=ii+1; end endenduw0=sqrt(pvv/(m-n-ii)*1.0e0);return 计算各点误差椭圆-ellipse函数（98页）function ellipse(uw0,fit2)glosbal ed dd pn c x y ai bi fifprintf(fit2, parameter of error ellipsen);fprintf(fit2, pn(i) mx my mm a b fi

17、n);n=2.0*dd;maxmm=0.0;smm=0.0;for i=1:dd ii=ed+i; di=(2*i-2)*(n-(2*i-1)/2.0); dj=(2*i-1)*(n-i); q1=c(di+2*i-1); q2=c(dj+2*i); q3=c(di+2*i); d1=sqrt(abs(q1+q2-q3); d=uw0*dl; xx=q1-q2; yy=2*q3; zz=q1+q2; mx1=sqrt(q1); my1=sqrt(q2); mm1=sqrt(zz); if(abs(xx)=0&yy=0) fi(i)=fi(i)/2.0; elseif(xx=0&yy=0) fi

18、(i)=(fi(i)+360)/2.0; elseif(xxmaxmm) maxmm=mm; i1=ii; end smm=smm+mm; ai(i)=uw0*a1; bi(i)=uw0*b1; ab=ai(i)-bi(i); fprintf(fit2, %10d %10.3f %10.3f %10.3f %10.3f %10.3f %10.3fn,pn(ii),mx,my,mm,ai(i),bi(i),fi(i);endsmm=smm/dd;fprintf(fit2, mse of unit weight= %9.6fn,uw0);fprintf(fit2, the maximum stat

19、ion error mm= %8.3f(cm) pn= %4dn,maxmm,pn(i1);fprintf(fit2, the average station error mm= %8.3fn,smm);return边长观测值平差值改正数及精度评定-adjs函数（101页）function adjs(uw,fit2)global ed dd sd pn m1 e d sid x y fprintf(fit2,adjusted sidesn);fprintf(fit2,i pn(e) pn(d) side vs(cm) side+vs msn);for i=1:m1 dx=x(d(i)-x(e(

20、i); dy=y(d(i)-y(e(i); if(e(i)d(i) maa,mss=trel(uw,e(i),d(i); else maa,mss=trel(uw,d(i),e(i); end ss=dx*dx+dy*dy; ss=sqrt(ss); vs=(ss-sid(i)*100; sid1=sid(i)+vs/100; fprintf(fit2, %3d %8d %8d %15.4f %10.4f %15.4f %8.2fn,i,pn(e(i),pn(d(i),sid(i),vs,sid1,mss);endreturn方向观测值的平差与精度评定-adjd函数（102页） functio

21、n adjd(uw,fit2)global ed dd dd1 sd pn ni si e d g f s t netglobal m1 m2 m3 ms pp md ma x0 y0 x y sid dir az cglobal a q1 pa3 qls wfprintf(fit2,adusted directions and their accuracyn);fprintf(fit2, i pn(g) pn(f) dir vd() dir+vdn);lo=206264.8062470964;for j=1:m2 q1(j+m1)=md*md;endfor i=1:sd zi=0.0; za

22、l=0.0; for j=si(i):si(i)+ni(i)-1 dx=x(f(j)-x(g(j); dy=y(f(j)-y(g(j); a0=alfa(dx,dy); z1=a0-dir(j); if(z1=2.0*pi) ss=ss-2*pi; end vd=(a0-ss)*lo; dir1=dir(j)+vd/lo; dir1=rad_dms(dir1); dir(j)=rad_dms(dir(j); fprintf(fit2, %3d %8d %8d %14.5f %9.2f %16.5fn,j,pn(g(j),pn(f(j),dir(j),vd,dir1); endendfprint

23、f(fit2,adusted directionsn);fprintf(fit2, i pn(g) pn(f) dir az man);for i=1:sd for j=si(i):si(i)+ni(i)-1 dx=x(f(j)-x(g(j); dy=y(f(j)-y(g(j); if(f(j)g(j) maa,mss=trel(uw,f(j),g(j); else maa,mss=trel(uw,g(j),f(j); end a0=alfa(dx,dy); if(j=si(i) a00=a0; end ss=a0-a00; if(ss0) ss=ss+2.0*pi; end ss=rad_d

24、ms(ss); a0=rad_dms(a0); fprintf(fit2, %3d %10d %10d %14.5d %14.5d %10.4fn,j,pn(g(j),pn(f(j),ss,a0,maa); endendreturn 方位角观测值改正数与精度评定-adja函数（104页） function adja(uw,fit2)global m3 x y s t az l0=206264.8062470964;fprintf(fit2,adjusted azimath n i pn(s) pn(t) az va() az+va ma);for i=1:m3 dx=x(t(i)-x(s(i)

25、; dy=y(t(i)-y(s(i); a0=alfa(dx,dy); if(t(i)pi) va=(a0-az(i)-2*pi)*lo; else va=(a0-az(i)*lo; end az1=az(i)+va/lo; az1=wg(az1); az(i)=wg(az(i); fprintf(fit2, %3d %8d %8d %13.5f %8.2f %12.5f %8.2f,i,pn(s(i),pn(t(i),az(i),va,az1,maa);endreturn 调用的的trel函数 function maa,mss=trel(uw,rr,tt)global ed dd x0 y0

26、 pn cn=2*(dd+ed);dx=x0(tt)-x0(rr);dy=y0(tt)-y0(rr);ss=sqrt(dx*dx+dy*dy);a0=alfa(dx,dy);aij=2062.648*sin(a0)/ss;bij=-2062.648*cos(a0)/ss;a=2*rr-1;b=a+1;f=2*tt-1;g=f+1;da=(a-1)*(n-a/2.0);db=(b-1)*(n-b/2.0);df=(f-1)*(n-f/2.0);dg=(g-1)*(n-g/2.0);q1=c(da+a)+c(df+f)-2.0*c(da+f);q2=c(db+b)+c(dg+g)-2.0*c(db

27、+g);q3=c(da+b)+c(df+g)-c(da+g)-c(db+f);x=q1-q2;y=2.0*q3;z=q1+q2;qs=q1*cos(a0)2+q2*sin(a0)2+q3*sin(2.0*a0);qa=q1*aij*aij+q2*bij*bij+2.0*q3*aij*bij;ms=sqrt(qs)*uw;maa=sqrt(qa)*uw;return附录资料：MATLAB Cell函数使用技巧谈谈MATLAB中cell函数如果p为一个数，那么h(1)=p,是没有问题的。如果p为一个向量，那么h(1,:)=p是没有问题的。如果p是一个矩阵的话，上面的两种赋值方法都是会有错误的。那么

28、要如何处理呢？这时就用到了cell数据类型了。cell的每个单元都可以存储任何数据，比如传递函数等。当然，存储矩阵更是没有问题的了。但是用cell数据类型之前，要先初始化。a=cell(n,m)那么就把a初始化为一个n行m列的空cell类型数据。如何赋值呢？a1,1=rand(5)那么a的1行1列的单元中存储的就是一个随机的55的方阵了。那么要用第一个单元中的方阵中的某个值呢？可以如下引用：a1,1(2,3)就可以了，引用cell单元时要用,再引用矩阵的某个数据就要用()了。cell单元中的每个单元都是独立的，可以分别存储不同大小的矩阵或不同类型的数据。下面举个例子：a=cell(2,2);%

29、预分配a1,1=cellclass;a1,2=1 2 2;a2,1=a,b,c;a2,2=9 5 6;a1,1ans =cellclassa1,2ans = 1 2 2a2,:ans =abcans = 9 5 6 b=a1,1b =cellclass元胞数组：元胞数组是MATLAB的一种特殊数据类型，可以将元胞数组看做一种无所不包的通用矩阵，或者叫做广义矩阵。组成元胞数组的元素可以是任何一种数据类型的常数或者常量，每一个元素也可以具有不同的尺寸和内存占用空间，每一个元素的内容也可以完全不同，所以元胞数组的元素叫做元胞（cell）。和一般的数值矩阵一样，元胞数组的内存空间也是动态分配的。（1）

30、元胞数组的创建 a=matlab,20;ones(2,3),1:10a = matlab 20 2x3 double 1x10 double b=matlab,20;ones(2,3),1:10b = matlab 20 2x3 double 1x10 double c=10c = 10c(1,2)=2c = 10 2c(2,2)=5c = 10 2 5isequal(a,b)ans = 1whosName Size Bytes Class Attributesa 2x2 388 cell ans 1x1 1 logical b 2x2 388 cell c 2x2 208 cell 用cell函数创建元胞数组，创建的数组为空元胞。cell函数创建空元胞数组的主要目的是为数组预先分配连续的存储空间，节约内存占用，提高执行效率。 a=cell(1)a = b=cell(1,2)b = c=cell(3,3)c = d=cell(2,2,2)d(:,:,1) =