C语言版测绘程序.doc

测量平差程序设计实习报告 中国地质大学（武汉）郑杨琳测量平差程序设计1. 角度（度分秒）到弧度AngleToRadian#define PI 3.14159265double 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.14159265double 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+) Ai=(type*)malloc(m*sizeof(type); 5.释放开辟的二维数组的空间#define FreeSpace(A,i,m)for(i=0;i<m;i+) free(Ai); free(A); 注意:释放空间与开辟空间相反,释放空间是先释放列,后释放行.6.矩阵求转置transformmatrixvoid transformmatrix（double *A，double *B，int i,int j）int m,n;for(m=0;m<=i;m+)for(n=0;n<=j;n+)Bnm=Amn: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+)Cmn=0;for(p=0;p<k;p+)Cmn+=Amp*Bpn: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=Nii;for(k=0;k<i;k+)s-=Tki*Tki;Tii=sqrt(s)for(j=i+1;j<n;j+)s=Nij;for(k=0;k<i;k+)s-=Tki*Tkj;Tij=s/Tii;for(i=0;i<n;i+)for(j=0;j<n;j+)Tij=0;for(i=n-1;i>=0;i+)rii=1/Tii;for(j=i+1;j<n;j+)s=0;for(k=i;k<j-1;k+)s-=rik*Tkj;rij=s/Tii;for(i=0;i<n;i+)for(j=0;j<n;j+)rij=0;transformmatrix(r,rt,n,n)mulmatrix(r,rt,Q,n,n)9.平差主程序之读入数据typedef struct POINTchar name8;double x,y;int type;POINT;typedef struct READVALUEPOINT *begin;POINT *end;double value;READVALUE;POINT *GETPOINT(char *name,POINT *pPoint,int nPoint)int i;for(i=0;i<nPoint;i+)if (strcmp(pPointi.name,name)=0)return (pPoint+i) for(i=0;i<nPoint;i+)if(pPointi=NULL)strcmp(pPointi.name,name);pPointi.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 begin8,end8;FILE *fp=0;fp=fopen(“c:datt1.txt”,”r”)fscanf(fp,”%d,%d,%d,%dn”,&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,%lfn”,&mo,&mf,&ms);for(i=0;i<nKnownPoint;i+)fscanf(fp,”%s,%lf,%lfn”,pPointi.name,&pPointi.x,&pPointi.y);type=1;for( ;i<nPoint;i+)pPointi.name=NULL; pPointi.x=0;pPointi.y=0;pPointi.type=0; for(i=0;i<nDirect;i+)fscanf(fp,”%s,%s,%lfn”,begin,end,&pDirecti.value);pDirecti.begin=GetPoint(begin,pPoint,nPoint);pDirecti.end=GetPoint(end,pPoint,nPoint);for(i=0;i<nDistance;i+)fscanf(fp,”%s,%s,%lfn”,begin,end,&pDistancei.value);pDistancei.begin=GetPoint(begin,pPoint,nPoint);pDistancei.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.坐标概算全方向法子函数取出观测方向GetAllDirectint GetAllDirect(char *name,int nDirect,READVALUE *pDirect, READVALUE *pStation)int i,nCount=0;for(i=0;i<nDirect;i+)if(strcmp(pDirecti.begin->name,name)=0)pStationnCount.begin=p(pDirectnCount.begin;pStationnCount.end=p(pDirectnCount.end;pStationnCount.value=p(pDirectnCount.value; nCount+;return nCount;坐标概算全方向法子程序实现流程(coordinate)coordinate (入口参数设置)READVALUE pStation50,pObject50;int nCount,i,j,k,m,n,p,nobject;for(i=0;i<nPoint;i+)nCount=GetAllDirect(pPointi.name,nDirect,pStation)if(nCount>1)|( nCount=1)for(j=0;j<nCount;j+)if(pStationj.end->type=1)for(k=0;k<nCount;k+)if(pStationk.end->type=0) nobject=GetAllDirect(pStationj.end->name,nDirect,pDirect,pobject)m=-1;n=-1;for(p=0;p<nobject;p+)if(strcmp(pobjectp.end->name,pPointi.name)=0)m=p; if(strcmp(pobjectp.end->name,pStationk.end->name)=0)n=p;if(m>=0&&n>=0)pPointi=pStationk.end-pStationj.end;pStationj.end=pObjectm.value-pObjectn.value; Xe=pPointi.x; Ye=pPointi.y; Xf=pStationj.end->x; Yf=pStationj.end->y; Lef=pStationj.value; Leg=pStationk.value; Lfe=pObjectm.value; Lfg=pObjectn.value; Qianfang(Xe,Xf,Ye,Yf,Lef,Leg,Lfe,Lfg,*Xg,*Yg;) pStationk.end->x=*xg; pStationk.end->y=*yg; pStationk.end.type=2; 13.坐标增量法(calcoordinate)子函数由端点名称得边长值的函数GetDistancedouble GetDistance(char *begin,char *end,int nDistance,READVALUE *pDistance)int i;for(i=0;i<nDistance;i+) if(strcmp(pDistancei.begin->name,begin)=0&&strcmp(pDistancei.end->name,end=0)|(strcmp(pDistancei.begin->name,end)=0&&strcmp(pDistancei.end,begin)=0)return pDistancei.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 pStation50; for(i=0;i<nPoint;i+) if(pPointi.type>0) nCount=GetAllDirect(pPointi.name,nDirect,pDirect,pStation50); for(j=0;j<nCount;j+) if(pStationj.end->type>0)m=j; if(m!=-1) for(k=0;k<nCount;k+) if(pStationk.end->type=0) x1=pPointi.x; y1=pPointi.y; x2=pStationj.end->x; y2=pStationj.end->y; A0=Bearing(x1,y1,x2,y2); A=A0-(DMSToRAD(pStationm.value)-DMSToRAD(pStationk.value); if(A<0)A=A+2*PI; if(A>2*PI)A=A-2*PI; S=GetDistance(pPointi,pStationk.end,nDistance,pDistance); if(S<0)continue; else dx=S*cos(A); dy=S*sin(A); pStationk.end->x=pPointi.x+dx; pStationk.end->y=pPointi.y+dy; pStationk.end->type=2; 14.高斯正反算高斯正算：#include <math.h>#include <stdio.h>#define PI 3.14159265double 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 -