obnorm.m

利用ＭＡＴＬＡＢ实现了对测量中的平面网

function obnorm
global 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 fit2
global a ql pa3 qls w
lo=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;
end
c(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));
     ql(i)=(ms^2+(ss*pp*0.0001)^2);
 end

ql(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):si(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/ql(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)/ql(ii);
       pa3(i,jj+2)=a(ii,7);
       pa3(i,jj+3)=a(ii,8)/ql(ii);
       pa3(i,2)=pa3(i,2)+a(ii,2)/ql(ii);
       pa3(i,4)=pa3(i,4)+a(ii,4)/ql(ii);
       jj=jj+4;
   end
   pa3(i,1)=a(ii,1);
   pa3(i,3)=a(ii,3);
   qls(i)=-zal;
end

for 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.0;
     else
        a(ii,9)=(a0-az(i))*lo*100.0;
     end
     ql(ii)=ma*ma;
 end

for 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)/ql(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)/ql(i);
        end
    end
end

if(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.0);
            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
    end
end
return