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

观测数据读入程序，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('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=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=chkdat(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(m3>0) ma=fscanf(fit1,'%f',1); a=fscanf(fit1,'%d %d %f',3*m3); for i=1:m3 s(i)=a(3*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(kk>0) 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 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; sina=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(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<=0.0) z1=z1+2.0*pi; end zal=zal+1./ql(m1+j); z0=z0+z1/q1(m1+j); end z0=z0/zal; for j=si(i):si(i)+ni(i)-1 dx=x0(f(j)-x0(g(j); dy=y0(f(j)-y0(g(j); ss=dx*dx+dy*dy; a0=alfa(dx,dy); ai=-dy/ss*lo; bi=dx/ss*lo; ii=m1+j; a(ii,1)=2*g(j)-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*f(j)-1-2*ed+1.0e-9; a(ii,6)=ai; a(ii,7)=a(ii,5)+1; a(ii,8)=bi; ss=dir(j)+z0; if(ss>=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)/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-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); if(kk<=0|jj>kk) continue; end c(di+kk)=c(di+kk)+a(i,2*k)*a(i,2*j)/q1(i); end endendif(m2>0) %和误差方程形成法方程 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(kk<=0|jj>kk) continue; end c(di+kk)=c(di+kk)+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');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/2.0); for j=i:n1 if(j=i) c(di+j)=0.00000001; else c(di+j)=0.0; end endendif(m1>0) adjs(uw0,fit2);endif(m2>0) adjd(uw0,fit2);endif(m3>0) 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 pn 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(j<i) dxy(i)=dxy(i)-c(dj+i)*w(j); else dxy(i)=dxy(i)-c(di+j)*w(j); end end if(abs(dxy(i)>1.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*i);endx0(1:sd)=x(1:sd);y0(1:sd)=y(1:sd);for i=1:sd if(i<=ed) vx(i)=0.0; vy(i)=0.0; else vx(i)=dxy(2*(i-ed)-1); vy(i)=dxy(2*(i-ed); endendfprintf(fit2,' adjusted coordinatesn');fprintf(fit2,' pn vx x vy yn');for i=1:sd fprintf(fit2,' %6d %8.4f %14.4f %8.4f %14.4fn',pn(i),vx(i),x(i),vy(i),y(i);endpvv=0.0;for i=1:n pvv=pvv+w(i)*dxy(i)*100.0;endm=m1+m2+m3;for i=1:m pvv=pvv+a(i,9)*a(i,9)/ql(i);endif(m2>0) 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 fin');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)<1d-10) fi(i)=sign(90.0,q3); else fi(i)=atan(yy/xx)*57.2958; end if(xx>=0&yy>=0) fi(i)=fi(i)/2.0; elseif(xx>=0&yy<=0) fi(i)=(fi(i)+360)/2.0; elseif(xx<0) fi(i)=(fi(i)+180)/2.0; end ww=sqrt(xx*xx+yy*yy); a1=sqrt(zz+ww)/2.0); b1=sqrt(zz-ww)/2.0); ab1=a1-b1; mx=uw0*mx1; my=uw0*my1; mm=uw0*mm1; if(mm>maxmm) 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 station 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(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页） function 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; zal=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<0.0) z1=z1+2.0*pi; end zal=zal+1./ql(m1+j); zi=zi+z1/ql(m1+j); end if(ni(i)=0) zi=zi/zal; end 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); ss=dir(j)+zi; if(ss>=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); endendfprintf(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(ss<0) ss=ss+2.0*pi; end ss=rad_dms(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); dy=y(t(i)-y(s(i); a0=alfa(dx,dy); if(t(i)<s(i) maa,mss=trel(uw,t(i),s(i); else maa,mss=trel(uw,s(i),t(i); end if(a0-az(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 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+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是一个矩阵的话，上面的两种赋值方法都是会有错误的。那么要如何处理呢？这时就用到了cell数据类型了。cell的每个单元都可以存储任何数据，比如传递函数等。当然，存储矩阵更是没有问题的了。但是用cell数据类型之前，要先初始化。a=cell(n,m)那么就把a初始化为一个n行m列的空cell类型数据。如何赋值呢？a1,1=rand(5)那么a的1行1列的单元中存储的就是一个随机的5×5的方阵了。那么要用第一个单元中的方阵中的某个值呢？可以如下引用：a1,1(2,3)就可以了，引用cell单元时要用,再引用矩阵的某个数据就要用()了。cell单元中的每个单元都是独立的，可以分别存储不同大小的矩阵或不同类型的数据。下面举个例子：a=cell(2,2);%预分配a1,1='cellclass'a1,2=1 2 2;a2,1='a','b','c'a2,2=9 5 6;>>a1,1ans =cellclass>>a1,2ans = 1 2 2>>a2,:ans =abcans = 9 5 6>> b=a1,1b =cellclass元胞数组：元胞数组是MATLAB的一种特殊数据类型，可以将元胞数组看做一种无所不包的通用矩阵，或者叫做广义矩阵。组成元胞数组的元素可以是任何一种数据类型的常数或者常量，每一个元素也可以具有不同的尺寸和内存占用空间，每一个元素的内容也可以完全不同，所以元胞数组的元素叫做元胞（cell）。和一般的数值矩阵一样，元胞数组的内存空间也是动态分配的。（1）元胞数组的创建>> 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 = 10>>c(1,2)=2c = 10 2>>c(2,2)=5c = 10 2 5>>isequal(a,b)ans = 1>>whosName 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) =