code

#include"stdio.h"
 #include"math.h"
 #include "mem.h"
 #include "stdlib.h"

 void Matr_muti(double *destine,double *source1,int row1,int col1,
double *source2,int row2,int col2);
 void Matr_rotation(double *destine,double *source,int row,int col);
 void Matr_converse(double *a,int w);

 main()
 {
int i,j,k;
int *knum;
FILE *in,*out;
double db[50][50],dx[50],dy[50];
int n,m;
double x,y,hr,hg,xg,yg,dx1,dy1;
double *pes0,*pes0k,*pes1,*v;
double vkmax,vcmax,muk,muc;
double *xx,*xk,*xt,*mid0,*mid1,*b;
char signal;
double dis,s[50],limit4[50],limit3[50];
     double *e,*mid;
      /*----------- edit data -------------*/
  in = fopen("d:\data.txt","r");
if(in==NULL) {printf("cannot open data file");
      exit(1);}
  out=fopen("d:\resu.c","a");
if(out==NULL) {printf("cannot open result2 file");
      exit(1);}
printf("\nInput the quantities of GPS-level points & known-points:\n");
scanf("\n%d %d",&n,&m);
printf("Input known-point numbers:");
knum=(int *)malloc(m*sizeof(int));
fprintf(out,"\nknown-point numbers:\n");
for(i=0;i<m;i++)
{
scanf("%d",&k);
*(knum+i)=k;
fprintf(out,"%d  ",k);
}
xx=(double *)malloc(n*6*sizeof(double));
xk=(double *)malloc(m*6*sizeof(double));
pes0=(double *)malloc(n*sizeof(double));
pes1=(double *)malloc(n*sizeof(double));
pes0k=(double *)malloc(m*sizeof(double));
v=(double *)malloc(n*sizeof(double));

 /*-----Compute the height error of GPS-level points-----*/
for(i=0;i<n;i++)
{
for(j=0;j<4;j++)
fscanf(in,"%lf",&db[j]);
}

 /*------compute the center-point's coordinates------*/
xg=0.0;yg=0.0;
for(i=0;i<n;i++)
   {  xg=xg+db[0];
      yg=yg+db[1];
   }
xg=xg/n;yg=yg/n;
 /*------compute the difference of coordinates-------*/
for(i=0;i<n;i++)
   {   dx=(db[0]-xg)*0.001;
dy=(db[1]-yg)*0.001;
   }

 /*-----Build all GPS-level points' matrix X------*/
for(i=0;i<n;i++)
{
*(pes0+i)=db[3]-db[2];
*(xx+i*6)=1;
*(xx+i*6+1)=dx;
*(xx+i*6+2)=dy;
*(xx+i*6+3)=dx*dx;
*(xx+i*6+4)=dy*dy;
*(xx+i*6+5)=dx*dy;
}
 /*-----Build known-point's matrix------*/
for(i=0;i<m;i++)
{
k=*(knum+i);
for(j=0;j<6;j++)
{
*(xk+i*6+j)=*(xx+k*6+j);
}
*(pes0k+i)=*(pes0+k);
}
 /*-----Compute the parametres------*/
xt=(double*)malloc(6*m*sizeof(double));
mid0=(double*)malloc(6*6*sizeof(double));
mid1=(double*)malloc(6*m*sizeof(double));
b=(double*)malloc(6*sizeof(double));

Matr_rotation(xt,xk,m,6);
Matr_muti(mid0,xt,6,m,xk,m,6);
Matr_converse(mid0,6);
Matr_muti(mid1,mid0,6,6,xt,6,m);
Matr_muti(b,mid1,6,m,pes0k,m,1);

 /*-----Compute the height-error of GPS-level points by the model-----*/
Matr_muti(pes1,xx,n,6,b,6,1);
for(i=0;i<n;i++)
*(v+i)=(*(pes1+i)-*(pes0+i))*1000.0;
muk=0.0;muc=0.0;vkmax=0.0;vcmax=0.0;

  for(j=0;j<m;j++)
{k=*(knum+j);
muk=muk+*(v+k)*(*(v+k));
if(fabs(*(v+k))>fabs(vkmax)) vkmax=*(v+k);
}
  for(i=0;i<n;i++)
{ muc=muc+*(v+i)*(*(v+i));
if(fabs(*(v+i))>fabs(vcmax)) vcmax=*(v+i);
}
muc=muc-muk;
muk=sqrt(muk/(m-1));
muc=sqrt(muc/(n-m-1));

 /*------compute limited errors---------*/
      for(i=0;i<n;i++)
{s=100.0;
for(k=0;k<m;k++)
      if(i!=*(knum+k))
{dis=(db[0]-db[*(knum+k)][0])*(db[0]-db[*(knum+k)][0])
+(db[1]-db[*(knum+k)][1])*(db[1]-db[*(knum+k)][1]);
dis=sqrt(dis)/1000.0;
if(dis<s) s=dis;
     }
 limit4=sqrt(s)*20.0;
 limit3=sqrt(s)*12.0;
}
  /*-----output the checked result------*/
       fprintf(out,"\npoint   v    limited3  limit4  distance");
       for(i=0;i<n;i++)
{if(*(v+i)>limit4) printf("point %d beyond",i);
fprintf(out,"\n%3d   %5.1f",i,*(v+i));
fprintf(out,"  %5.1f   %5.1f  %5.1f",limit3,limit4,s);
}
fprintf(out,"\nmu(known):%5.1f  vmax:%5.1f\n",muk,vkmax);
fprintf(out,"mu(unknow):%5.1f  vmax:%5.1f\n",muc,vcmax);
      printf("compute points' hr? (y/n)");
      signal=getch();
      while(signal=='y')
{
  printf("\ninput x y hg:\n");
  scanf("%lf %lf %lf",&x,&y,&hg);
  dx1=0.001*(x-xg);
  dy1=0.001*(y-yg);
  *pes1=*b+*(b+1)*dx1+*(b+2)*dy1+*(b+3)*dx1*dx1+*(b+4)*dy1*dy1+*(b+5)*dx1*dy1;
  hr=hg-*pes1;
  printf("\nx=%lf",x);
  printf("\ny=%lf",y);
  printf("\nhr=%lf",hr);
  printf("continute to another point?(y/n)");
  signal=getch();
 }

      }
 /*------two matrixs multipate each other-------*/
 void Matr_muti(double *destine,double *source1,int row1,int col1,
 double *source2,int row2,int col2)
 {
int i,j,k;
for(i=0;i<row1;i++)
  for(j=0;j<col2;j++)
  {
   *(destine+i*col2+j)=0;
   for(k=0;k<col1;k++)
   *(destine+i*col2+j)=*(destine+i*col2+j)+(*(source1+i*col1+k))*
     (*(source2+k*col2+j));
}
 }
 void Matr_rotation(double *destine,double *source,int row,int col)
 {
int i,j;
for(i=0;i<row;i++)
for(j=0;j<col;j++)
*(destine+j*row+i)=*(source+i*col+j);
 }
 void Matr_converse(double *a,int w)
 {
   int i,j,k;
   double p,q,*h;
   h=(double*)malloc(w*sizeof(double));
   for(k=w;k>=1;k--)
     { p=*a;              printf("p=%f k=%d ",p,k);
       if(p<=0)
{free(h);
 return(0); }
       for(i=2;i<=w;i++)
 {  q=*(a+(i-1)*w);
 if(i>k)
     h[i-1]=q/p;
 else h[i-1]=-q/p;
 for(j=2;j<=i;j++)
     *(a+(i-2)*w+j-2)=*(a+(i-1)*w+j-1)+q*h[j-1];
 }
*(a+(w-1)*w+w-1)=1/p;
for(i=2;i<=w;i++)
   *(a+(w-1)*w+i-2)=h[i-1];
     }
     for(i=0;i<w;i++)
for(j=0;j<w;j++)
if(j>i) *(a+i*w+j)=*(a+j*w+i);
     free(h);
     return(1);
 }



来点说明好不好

怎么用

就是啦，最好给点说明，大家都好迷茫啦

看看再说，只要贴就是好的。

kk

<img src="images/post/smile/dvbbs/em02.gif" />