测绘常用程序C语言.doc

测量平差程序设计
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.

2. 弧度换角度(度分秒) 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;
}

3. 已知两点求坐标方位角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);
}

- 11 -