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 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(S1e-10) return 0;T=asin(Dy/S);if(Dx0&(Dy0)|T0) T=2*PI+T;return T;4.开辟二维数组的动态空间的宏#include #define NewArray2D(type,A,i,n,m)A=(type*)malloc(n*sizeof(type*); for(i=0;im;i+) Ai=(type*)malloc(m*sizeof(type); 5.释放开辟的二维数组的空间#define FreeSpace(A,i,m)for(i=0;im;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;mi;m+)for(n=0;nj;n+)Cmn=0;for(p=0;pk;p+)Cmn+=Amp*Bpn:8.矩阵求逆(countermatrix)#include void countermatrix(double *T, double *s, double *r, double *Q,double *N, double *rt,int n)for(i=0;in;i+)s=Nii;for(k=0;ki;k+)s-=Tki*Tki;Tii=sqrt(s)for(j=i+1;jn;j+)s=Nij;for(k=0;ki;k+)s-=Tki*Tkj;Tij=s/Tii;for(i=0;in;i+)for(j=0;j=0;i+)rii=1/Tii;for(j=i+1;jn;j+)s=0;for(k=i;kj-1;k+)s-=rik*Tkj;rij=s/Tii;for(i=0;in;i+)for(j=0;jn;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;inPoint;i+)if (strcmp(pP,name)=0)return (pPoint+i) for(i=0;i0)pPoint=(POINT*)malloc(nDirect*sizeof(POINT);if(nDirect0)pDirect=(READVALUE*)malloc(nDirect*sizeof(READVALUE);if(nDistance0)pDistance=(READVALUE*)malloc(nDistance*sizeof(RAADVALUE);fscanf(fp,”%lf,%lf,%lfn”,&mo,&mf,&ms);for(i=0;inKnownPoint;i+)fscanf(fp,”%s,%lf,%lfn”,pP,&pPointi.x,&pPointi.y);type=1;for( ;inPoint;i+)pP=NULL; pPointi.x=0;pPointi.y=0;pPointi.type=0; for(i=0;inDirect;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;inDistance;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 int checkangle(double angle)int M,S;double MS;if(angle=0&angle360)MS=angle-(int)(angle);if(M6)S=(int)(MS*1000);if(S%106)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-2*PI)|(C0&C-PI&CPI&C2*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;iname,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;i1)|( nCount=1)for(j=0;jtype=1)for(k=0;ktype=0) nobject=GetAllDirect(pStationj.end-name,nDirect,pDirect,pobject)m=-1;n=-1;for(p=0;pname,pP)=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;iname,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;i0) nCount=GetAllDirect(pP,nDirect,pDirect,pStation50); for(j=0;jtype0)m=j; if(m!=-1) for(k=0;ktype=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(A2*PI)A=A-2*PI; S=GetDistance(pPointi,pStationk.end,nDistance,pDistance); if(Sx=pPointi.x+dx; pStationk.end-y=pPointi.y+dy; pStationk.end-type=2; 14.高斯正反算高斯正算：#include #include #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*