C语言版测绘程序

1、测量平差程序设计1. 角度(度分秒)到弧度 An gleToRadian#defi ne PI 3.14159265double AngleToRadian (double angle)int D , M ;double S,radia n, degree, an gle,MS;D=int (angle+0.3);MS=angle-D ;M=i nt (MS) *100+0.3 );S= (MS*100-M ) *100;degree=D+M/60.0+S/3600.0;radia n=degree*PI/180.0;return radia n;注意:防止数据溢出，要加个微小量，例如0.3.

2、2. 弧度换角度(度分秒)RadianToAngle#defi ne PI 3.14159265double Radia nToA ngle(double radia n)int D,M;double S,radia n, degree,MS,a ngle;degree=radia n*180/PI;D=in t(degree);MS=degree-D;M=i nt(MS*60);S=(MS*60-M)*60;an gle=D+M/100.0+S/10000.0;return an gle;3. 已知两点求坐标方位角Azimuth#in clude double Azimuth(double

3、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=asi n( Dy/S);if(Dx0&(Dy0)|T0)T=2*PI+T;return T;4开辟二维数组的动态空间的宏#i nclude #defi ne NewArray2D(type,A,i, n, m)A=(type*)malloc( n*sizeof(type*); for(i=0;im;i+)Ai=(type*)malloc(m*sizeof(type); 5.

4、释放开辟的二维数组的空间#defi ne FreeSpace(A,i,m)for(i=0;im;i+)free(A i);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+)Bn m=Am n:7. 矩阵相乘(mulmatrix)void mulmatrix(double *A,double *B,double *C,i nt i,

5、i nt j,i nt k)int m, n,p;for(m=0;mi;m+)for(n=0;nj;n+)Cm n =0;for(p=0;pk;p+)Cm n+=Amp*Bp n:8. 矩阵求逆(countermatrix)#in elude void coun termatrix(double *T, double *s, double *r, double *Q,double *N, double *rt,i nt n)for(i=0;i n;i+)s=Nii;for(k=0;ki;k+)s-=Tki*Tki;Tii=sqrt(s)for(j=i+1;j n;j+)s=Nij;for(k=

6、0;ki;k+)s-=Tki*Tkj;Tij=s/Tii;for(i=0;i n;i+)for(j=0;j=0;i+)rii=1/Tii;for(j=i+1;j n ;j+)s=0;for(k=i;kj-1;k+)s-=rik*Tkj;rij=s/Tii;for(i=0;i n;i+)for(j=0;j n;j+)rij=0;tran sformmatrix(r,rt, n,n) mulmatrix(r,rt,Q, n,n)9平差主程序之读入数据typedef struct POINTchar n ame8;double x,y;int type;POINT;typedef struct RE

7、ADV ALUEPOINT *begi n;POINT *en d;double value;READV ALUE;POINT *GETPOINT(char *n ame,POINT *pPoi nt,i nt nPoi nt) int i;for(i=0;i nPoi nt;i+)if (strcmp(pPo in ti .n ame ,n ame)=0)return (pPoi nt+i)for(i=0;i0)pPoi nt=(POINT*)malloc( nDirect*sizeof(POINT);if(n Direct0)pDirect=(READV ALUE*)malloc(nDir

8、ect*sizeof(READV ALUE);if(n Dista nce0)pDistance=(READV ALUE*)malloc(nDistance*sizeof(RAADVALUE);fscanf(fp, %lf,%lf,%lfn ,&mo,&mf,&ms);for(i=0;i nKnownPoin t;i+)fscanf(fp, %s,%lf,%lfn ,pP,&pPointi.x,&pPointi.y); type=1;for( ;i n Poi nt;i+)pPoi nti. name=NULL;pPoi nti.x=0;pPoi nti.y=O;pPo i

9、n ti.type=0;for(i=0;i nDirect;i+)fscanf(fp, %s,%s,%lfn ,begin,end,&pDirecti.value); pDirecti.begi n=GetPo in t(begi n, pPo in t, nPoin t);pDirecti.e nd=GetPo in t(e nd,pPo in t, nPoin t);for(i=0;i nDista nce;i+)fscanf(fp, %s,%s,%lfn ”,begin,end,&pDistancei.value);pDista ncei.beg in=GetPo in t(begi n

10、,pPo int,nPoin t);pDista ncei.e nd=GetPo in t(e nd,pPo in t, nPoin t);fclose(fp);10. 角度检验(checkangle)#in elude int checka ngle(double an gle)int M,S;double MS;if(a ngle=0&an gle360)MS=a ngle-(i nt)(a ngle);if(M6)S=(i nt)(MS*1000);if(S%106)return 1;return 0;11. 前方交会#defi ne Pl=3014159265/*此处调用程序角度换弧度

11、AngleToRadian*/,double DEF, double DFG,Qia nfan g(double XE, double YE, double XF, double YF, doubleDEGdouble 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/ta n(B)+XF/ta n(A)+YE-YF)/(1/ta n(A)+ 1/ta n(B);YG=(YE/ta n(B)+YF/ta n(A)-XE+XF

12、)/ (1/ta n( A)+ 1/ta n(B);12坐标概算全方向法子函数取出观测方向GetAIIDirectint GetAIIDirect(char *name,int nDirect,READVALUE *pDirect, READV ALUE *pStation)int i,nCoun t=0;for(i=0;in ame ,n ame)=0)pStati onnCoun t.begi n=p(pDirect nCoun t.begi n;pStati onnCoun t.e nd=p(pDirect nCoun t.e nd;pStati onnCoun t.value=p(pD

13、irect nCoun t.value;nCoun t+;retur n nCount;坐标概算全方向法子程序实现流程(coordi nate)coordinate (入口参数设置)READV ALUE pStatio n 50,pObject50;int nCoun t,i,j,k, m,n,p,no bject;for(i=0;i1)|( n Cou nt=1)for(j=0;jtype=1)for(k=0;ktype=O)no bject=GetAIIDirect(pStati on j.e nd- name, nDirect,pDirect, pobject) m=-1;n=_1;fo

14、r(p=0;p n ame,pPo in ti .n ame)=0)m=p;if(strcmp(pobjectp.e nd-n ame,pStatio n k.e nd-n ame)=O) n=p；if(m=0&n =0)pPo in ti=pStati on k.e nd-pStati on j.e nd;pStatio nj.e nd=pObjectm.value-pObject n .value;Xe=pPoi nti.x;Ye=pPoi nti.y;Xf=pStatio nj.e nd-x;Yf=pStatio nj.e nd-y;Lef=pStati on j.value;Leg=p

15、Statio n k.value;Lfe=pObjectm.value;Lfg=pObject n .value;Qia nfan g(Xe,Xf,Ye,Yf,Lef,Leg,Lfe ,L fg,*Xg,*Yg;)pStati on k.e nd-x=*xg;pStati on k.e nd-y=*yg;pStati on k.e nd.type=2;13.坐标增量法(calcoordinate)子函数由端点名称得边长值的函数GetDista needouble GetDistance(char *begin,char *end,int nDistance,READV ALUE *pDista

16、nee)int i;for(i=0;in ame,begi n)=0&strcmp(pDista ncei.e nd-n ame,e nd=0)| |(strcmp(pDista ncei.begi n-n ame,e nd)=0&strcmp(pDista ncei.e nd,begi n)=0) retur n pDista ncei.value;return -1;/*函数取出观测方向GetAIIDirect*/void calcoordinate(int nDirect,READV ALUE *pDirect,int nDistace,READV ALUE *pDistance,int

17、 nPoin t,POINT *pPoi nt)int nPoint,nCount,n Direct ,n Dista nee;int m=-1,i,j,k;double x1,y1,x2,y2,A0,A,S,dx,dy;READV ALUE*pDirect=NULL;READVALUE pStatio n 50;for(i=0;i0)n Cou nt=GetAIIDirect(pPoi nti. name, nDirect,pDirect,pStatio n50);for(j=0;jtype0)m=j;if(m!=-1)for(k=0;ktype=0)x1= pPoi nti. x;y1=p

18、Poi nti.y;x2=pStati on j.e nd-x;y2=pStati on j.e nd-y;A0=Beari ng(x1,y1,x2,y2);A=A0-(DMSToRAD(pStatio nm.value)-DMSToRAD(pStati on k.value); if(A2*PI)A=A-2*PI;S=GetDista nce(pPo in ti,pStati on k.e nd,n Dista nce,pDista nee);if(Sx=pPo in ti.x+dx;pStatio n k.e nd-y=pPo in ti.y+dy;pStatio nk.e nd-type

19、=2;14高斯正反算高斯正算：#in elude #i nclude #defi ne PI 3.14159265double DMSToRAD (double dDMS )int L1,L2 ;double T,L3;L仁(i nt) ( dDMS+0.3 );L2=(i nt)(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, d

20、ouble *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;B仁DMSToRAD(B);t=ta nB1;c=sqrt(e1*e1*cosB1*cos*B1);l=L-L0;N=a/(sqrt(1-e*e*si nB1*si nB1);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/

21、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=aO*B1-a2*(s in (2*B1)/2+a4*(si n( 4*B1)/4-a6*(si n( 6*B1)/6+a8*(s in (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*co