测绘常用程序C语言

测量平差程序设计

1. 角度（度分秒）到弧度AngleToRadian

#define PI 3.14159265

double AngleToRadian（double angle）

{

       int D，M；

       double S,radian,degree, angle,MS；

       D=int（angle+0.3）；

       MS=angle-D；

       M=int（（MS）*100+0.3）；

       S=（MS*100-M）*100;

       degree=D+M/60.0+S/3600.0;

       radian=degree*PI/180.0;

       return radian;

}

注意:防止数据溢出,要加个微小量,例如0.3.

1. 弧度换角度(度分秒) RadianToAngle

#define PI 3.14159265

double RadianToAngle(double radian)

{

     int D,M;

     double S,radian,degree,MS,angle;

     degree=radian*180/PI;

     D=int(degree);

     MS=degree-D;

     M=int(MS*60);

     S=(MS*60-M)*60;

     angle=D+M/100.0+S/10000.0;

     return angle;

}

1. 已知两点求坐标方位角Azimuth

#include <math.h>

double Azimuth(double xi,double yi,double xj,double yj)

{

       double Dx,Dy,S,T;

       Dx=xj-xi;

       Dy=yj-yi;

       S=sqrt(Dx*Dx+Dy*Dy);

       if(S<1e-10)  return 0;

       T=asin(Dy/S);

       if(Dx<0)  T=PI-T;

       if(Dx>0&&(Dy<0)||T<0)  T=2*PI+T;

return T;

}

4.开辟二维数组的动态空间的宏

#include <malloc.h>

#define NewArray2D(type,A,i,n,m){A=(type**)malloc(n*sizeof(type*));\

                                                  for(i=0;i<m;i++)\

                                                  A[i]=(type*)malloc(m*sizeof(type));\

                                                  }

5.释放开辟的二维数组的空间

#define FreeSpace(A,i,m){for(i=0;i<m;i++)\

                                     free(A[       i]);\

                                     free(A);\

                                    }

注意:释放空间与开辟空间相反,释放空间是先释放列,后释放行.

6.矩阵求转置transformmatrix

void transformmatrix（double **A，double **B，int i,int j）

{

     int m,n;

     for(m=0;m<=i;m++)

            for(n=0;n<=j;n++)

            {

                   B[n][m]=A[m][n]:

            }

}

7.矩阵相乘(mulmatrix)

void mulmatrix(double **A,double **B,double **C,int i,int j,int k)

{

     int m,n,p;

     for(m=0;m<i;m++)

            for(n=0;n<j;n++)

            {

                   C[m][n]=0;

for(p=0;p<k;p++)

{

       C[m][n]+=A[m][p]*B[p][n]:

}

              }

}

8.矩阵求逆(countermatrix)

#include <math.h>

void countermatrix(double **T, double **s, double **r, double **Q,double **N, double **rt,int n)

{

     for(i=0;i<n;i++)

     {

            s=N[i][i];

            for(k=0;k<i;k++)

            {

                   s-=T[k][i]*T[k][i];

            }

            T[i][i]=sqrt(s)

            for(j=i+1;j<n;j++)

            {

                   s=N[i][j];

                   for(k=0;k<i;k++)

                   {

s-=T[k][i]*T[k][j];

                     }

                     T[i][j]=s/T[i][i];

              }

       }

for(i=0;i<n;i++)

       for(j=0;j<n;j++)

       {

              T[i][j]=0;

       }

for(i=n-1;i>=0;i++)

{

       r[i][i]=1/T[i][i];

       for(j=i+1;j<n;j++)

       {

              s=0;

              for(k=i;k<j-1;k++)

              {

                     s-=r[i][k]*T[k][j];

              }

              r[i][j]=s/T[i][i];

       }

}

for(i=0;i<n;i++)

       for(j=0;j<n;j++)

       {

              r[i][j]=0;

       }

transformmatrix(r,rt,n,n)

mulmatrix(r,rt,Q,n,n)

}

9.平差主程序之读入数据

typedef struct POINT

{

       char name[8];

       double x,y;

       int type;

}POINT;

typedef struct READVALUE

{

       POINT *begin;

       POINT *end;

       double value;

}READVALUE;

POINT *GETPOINT(char *name,POINT *pPoint,int nPoint)

{

       int i;

       for(i=0;i<nPoint;i++)

       {

              if (strcmp(pPoint[i].name,name)==0)

              return (pPoint+i)

       }

      for(i=0;i<nPoint;i++)

       {

              if(pPoint[i]=NULL)

              strcmp(pPoint[i].name,name);

              pPoint[i].type=0;

              return(pPoint+i);

       }

}

double AngleToRadian（double angle）

{

       int D，M；

       double S,radian,degree, angle,MS；

       D=int（angle+0.3）；

       MS=angle-D；

       M=int（（MS）*100+0.3）；

       S=（MS*100-M）*100;

       degree=D+M/60.0+S/3600.0;

       radian=degree*PI/180.0;

       return radian;

}

main()

{

       POINT *pPoint=NULL;

       READVALUE *pDirect=NULL;

       READVALUE *pDistance=NULL;

       int nPoint,nKnownPoint,nDirect,nDistance,i;

       double mo,mf,ms;

       char begin[8],end[8];

       FILE *fp=0;

       fp=fopen(“c:\\dat\\t1.txt”,”r”)

       fscanf(fp,”%d,%d,%d,%d\n”,&nPoint,&nKnowPoint,&nDirect,&nDistance)

       if(nPoint>0)

       pPoint=(POINT*)malloc(nDirect*sizeof(POINT));

       if(nDirect>0)

       pDirect=(READVALUE*)malloc(nDirect*sizeof(READVALUE));

       if(nDistance>0)

       pDistance=(READVALUE*)malloc(nDistance*sizeof(RAADVALUE));

       fscanf(fp,”%lf,%lf,%lf\n”,&mo,&mf,&ms);

       for(i=0;i<nKnownPoint;i++)

       {

              fscanf(fp,”%s,%lf,%lf\n”,pPoint[i].name,&pPoint[i].x,&pPoint[i].y);

              type=1;  

       }

       for( ;i<nPoint;i++)

       {

              pPoint[i].name=NULL;

             pPoint[i].x=0;

              pPoint[i].y=0;

              pPoint[i].type=0;

       }

      for(i=0;i<nDirect;i++)

       {

              fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDirect[i].value);

              pDirect[i].begin=GetPoint(begin,pPoint,nPoint);

              pDirect[i].end=GetPoint(end,pPoint,nPoint);

       }

       for(i=0;i<nDistance;i++)

       {

              fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDistance[i].value);

              pDistance[i].begin=GetPoint(begin,pPoint,nPoint);

              pDistance[i].end=GetPoint(end,pPoint,nPoint);

       }

       fclose(fp);

}

10.角度检验（checkangle）

#include <math.h>

int checkangle(double angle)

{

int M,S;

double MS;

if(angle>=0&&angle<360)

{

       MS=angle-(int)(angle);

       if(M<6)

       {

              S=(int)(MS*1000);

              if(S%10<6)

              {

                     return 1;

              }

       }

}

return 0;

}

11.前方交会

#define PI=3014159265

/***此处调用程序角度换弧度AngleToRadian***/

Qianfang(double XE, double YE, double XF, double YF, doubleDEG, double DEF, double DFG, double DFE, double *DFE, double *DFG)

{

       double C,A,B;

       C=DGE-DGF;

       A=DEF-DEG;

       B=DFG-DFE;

       if((C<-PI&&C>-2*PI)||(C>0&&C<PI)

       {

              XG=(XE/tan(B)+XF/tan(A)-YE+YF)/(1/tan(A)+ 1/tan(B);

              YG=(YE/tan(B)+YF/tan(A)+XE-XF)/ (1/tan(A)+ 1/tan(B);

       }

       if((C>-PI&&C<0)||(C>PI&&C<2*PI))

       {

              XG=(XE/tan(B)+XF/tan(A)+YE-YF)/(1/tan(A)+ 1/tan(B);

              YG=(YE/tan(B)+YF/tan(A)-XE+XF)/ (1/tan(A)+ 1/tan(B);

       }

}

12.坐标概算全方向法

子函数取出观测方向GetAllDirect

int GetAllDirect(char *name,int nDirect,READVALUE *pDirect, READVALUE *pStation)

{

       int i,nCount=0;

       for(i=0;i<nDirect;i++)

       if(strcmp(pDirect[i].begin->name,name)==0))

       {    

              pStation[nCount].begin=p(pDirect[nCount].begin;

              pStation[nCount].end=p(pDirect[nCount].end;

              pStation[nCount].value=p(pDirect[nCount].value;

             nCount++;

       }

return nCount;

}

坐标概算全方向法子程序实现流程(coordinate)

coordinate (入口参数设置)

{

       READVALUE pStation[50],pObject[50];

       int nCount,i,j,k,m,n,p,nobject;

       for(i=0;i<nPoint;i++)

       {

              nCount=GetAllDirect(pPoint[i].name,nDirect,pStation)

              if((nCount>1)||( nCount=1))

              {

                     for(j=0;j<nCount;j++)

                     {

                            if(pStation[j].end->type==1)

                            {

                                   for(k=0;k<nCount;k++)

                                   {

                                          if(pStation[k].end->type==0)

                                                                                     nobject=GetAllDirect(pStation[j].end->name,nDirect,pDirect,pobject)

                                                   m=-1;

                                                   n=-1;

                                                   for(p=0;p<nobject;p++)

                                                   {

                                                        if(strcmp(pobject[p].end->name,pPoint[i].name)==0)

                                                        {

                                                               m=p;

                                                        }

                                                    if(strcmp(pobject[p].end->name,pStation[k].end->name)==0)

                                                        {

                                                               n=p;

                                                        }

                                                   if(m>=0&&n>=0)

                                                   {

                                                   pPoint[i]=pStation[k].end-pStation[j].end;

                                                   pStation[j].end=pObject[m].value-pObject[n].value;

                                                   {

Xe=pPoint[i].x;

                                         Ye=pPoint[i].y;

                                          Xf=pStation[j].end->x;

                                         Yf=pStation[j].end->y;

                                         Lef=pStation[j].value;

                                         Leg=pStation[k].value;

                                         Lfe=pObject[m].value;

                                         Lfg=pObject[n].value;

                                         Qianfang(Xe,Xf,Ye,Yf,Lef,Leg,Lfe,Lfg,*Xg,*Yg;)

                                         pStation[k].end->x=*xg;

                                         pStation[k].end->y=*yg;

                                         pStation[k].end.type=2;

                                    }

                                }

                          }

                     }

               }

          }

    }

 }

                                                  

13.坐标增量法(calcoordinate)

子函数由端点名称得边长值的函数GetDistance

double GetDistance(char *begin,char *end,int nDistance,READVALUE *pDistance)

{

       int i;

       for(i=0;i<nDistance;i++)

       {                   if((strcmp(pDistance[i].begin->name,begin)==0&&strcmp(pDistance[i].end->name,end==0)||(strcmp(pDistance[i].begin->name,end)==0&&strcmp(pDistance[i].end,begin)==0))

       return pDistance[i].value;

       }

return -1;

}

/***函数取出观测方向GetAllDirect***/

void calcoordinate(int nDirect,READVALUE *pDirect,int nDistace,READVALUE *pDistance,int nPoint,POINT *pPoint)

{

int nPoint,nCount,nDirect,nDistance;

      int m=-1,i,j,k;

      double x1,y1,x2,y2,A0,A,S,dx,dy;

      READVALUE*pDirect=NULL;

      READVALUE pStation[50];

      for(i=0;i<nPoint;i++)

      {

 if(pPoint[i].type>0)

            {

nCount=GetAllDirect(pPoint[i].name,nDirect,pDirect,pStation[50]);

                 for(j=0;j<nCount;j++)

                 {

if(pStation[j].end->type>0)m=j;

               if(m!=-1)

               {

for(k=0;k<nCount;k++)

              {

 if(pStation[k].end->type==0)

                     {

x1=pPoint[i].x;

                         y1=pPoint[i].y;

                         x2=pStation[j].end->x;

                         y2=pStation[j].end->y;

                         A0=Bearing(x1,y1,x2,y2);

                    A=A0-(DMSToRAD(pStation[m].value)-DMSToRAD(pStation[k].value));

                         if(A<0)A=A+2*PI;

                     if(A>2*PI)A=A-2*PI;

                         S=GetDistance(pPoint[i],pStation[k].end,nDistance,pDistance);

                         if(S<0)continue;

                         else

                         {

dx=S*cos(A);

                              dy=S*sin(A);

                              pStation[k].end->x=pPoint[i].x+dx;

                               pStation[k].end->y=pPoint[i].y+dy;               

                              pStation[k].end->type=2;}

                          }

                     }

               }

          }

    }

  }

}

14.高斯正反算

高斯正算：

#include <math.h>

#include <stdio.h>

#define PI 3.14159265

double DMSToRAD（double dDMS）

{

       int L1,L2；

       double T,L3；

       L1=(int)（dDMS+0.3）；

       L2=(int)((dDMS-L1)*100+0.3)；

       L3=((dDMS-L1)*100-L2)*100；

       T=(L1+L2/60.0+L3/3600.0)*PI/180.0;

       return T;

}

void PreGausePositive(double B,double L,double L0, double a, double b, double *N, double *l, double *c, double *t, double *X,double *B1)

{

   double a0,a2,a4,a6,a8,m0,m2,m4,m6,m8;

   double e,e1;

   e=(sqrt(a*a-b*b))/a;

   e1=(sqrt(a*a-b*b))/b;

   B1=DMSToRAD(B);

   t=tanB1;

   c=sqrt(e1*e1*cosB1*cos*B1);

   l=L-L0;

   N=a/(sqrt(1-e*e*sinB1*sinB1));

   m0=a*(1-e*e);

   m2=3/2*e*e*m0;

   m4=5/4*e*e*m2;

   m6=7/6*e*e*m4;

   m8=9/8*e*e*m6;

   a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;

   a2=m2/2+m4/2+15*m6/32+7/16*m8;

   a4=m4/8+3*m6/16+7*m8/32;

   a6=m6/32+m8/16;

   a8=m8/128;

   X=a0*B1-a2*(sin(2*B1))/2+a4*(sin(4*B1))/4-a6*(sin(6*B1))/6+a8*(sin(8*B1))/8;

}

Void BLToXY(double *x,double *y,double N,double l,double c,double t,double B1,double X)

  {

x=X+N*l*l*t*cosB1*cosB1*((3+l*l*cosB1*cosB1*(5-t*t+9*c*c+4*c*c*c*c)/4+l*l*\

cosB1*cosB1*(61-58*t*t+t*t*t*t)/30))/6;

y=N*l*cosB1(1+l*l*cosB1*((1+c*c-t*t)+l*l*cosB1*cosB1(5-18*t*t+t*t*t*t+14*c*c\

-58*t*t*c*c)));

}

高斯反算

void XYToBL(double x,double y,double L0,double a,double b,double q,double *B,\

double *L)

{

             double Bf,c,t,y,N,e1,e

             e=(sqrt(a*a-b*b))/(a*a);

             e1=(sqrt(a*a-b*b))/(b*b);

             for(Bf=0;;)

           {

                 t=tanBf;

                 c=e1*e1*cosBf;

                 N=a/(sqrt(1-e*e*sinBf*sinBf));

B=Bf-(1+c*c)*t*y*y/(2*N*N)*(1-y*y)/(12*N*N)*(15+3*t*t+c*c-9*t*t*c*c)-y*y/(30*N*N)*(61+90*t*t+45*t*t*t*t);

                 if(fabs(B-Bf)<q)

              break;

              Bf=B;

           }

L=L0+y/(N*cosBf)*(1-y*y/(6*N*N))*((1+2*t*t+c*c)-y*y/(20*N*N)*(5+6*c*c)+28*t*t+24*t*t*t*t+8*t*t*c*c)-y*y/(42*N*N)*(61+662*t*t+1320*t*t*t*t+720*t*t*t*t*t*t));

             B=RADToDMS(B);

             L=RADToDMS(L);

}