大地坐标与空间直角坐标的转换程序代码

大地坐标与空间直角坐标的转换程序代码大地坐标与空间直角坐标的转换程序代码大地坐标与空间直角坐标的转换程序代码大地坐标与空间直角坐标的转换程序代码编制仅供参考审核批准生效日期地址：电话：传真:邮编:#include""#include""#include""#include"iostream"#definePIdoublea,b,c,e2,ep2;intmain(){ intm,n,t; doubleRAD(doubled,doublef,doublem); voidRBD(doublehd); voidBLH_XYZ(); voidXYZ_BLH(); voidB_ZS(); voidB_FS(); voidGUS_ZS(); voidGUS_FS(); printf("大地测量学\n");sp1:printf("请选择功能:\n"); printf("1.大地坐标系到大地空间直角坐标的转换\n"); printf("2.大地空间直角坐标到大地坐标系的转换\n"); printf("3.贝塞尔大地问题正算\n"); printf("4.贝塞尔大地问题反算\n"); printf("5.高斯投影正算\n"); printf("6.高斯投影反算\n"); printf("0.退出程序\n"); scanf("%d",&m); if(m==0)exit(0);sp2:printf("请选择椭球参数(输入椭球序号):\n"); printf("1.克拉索夫斯基椭球参数\n"); printf("椭球参数\n"); printf("椭球参数\n"); printf("0.其他椭球参数(自行输入)\n"); scanf("%d",&n); switch(n) { case1:a=;b=;c=;e2=;ep2=;break; case2:a=;b=;c=;e2=;ep2=;break; case3:a=;b=;c=;e2=;ep2=;break; case0:{ printf("请输入椭球参数：\n"); printf("长半径a=");scanf("%lf",&a); printf("短半径b=");scanf("%lf",&b); c=a*a/b; ep2=(a*a-b*b)/(b*b); e2=(a*a-b*b)/(a*a); break; } default:printf("\n\n输入错误!\n请重新输入!\n\n");gotosp2; } while(1) { switch(m) { case1:BLH_XYZ();break; case2:XYZ_BLH();break; case3:B_ZS();break; case4:B_FS();break; case5:GUS_ZS();break; case6:GUS_FS();break; default:printf("\n\n输入错误!\n请重新输入!\n\n");gotosp1; } printf("是否继续进行此功能计算\n\n"); printf("（若继续进行此功能计算，则输入1；\n若选择其他功能进行计算，则输入2；\n若退出，则输入0.）\n"); scanf("%d",&t); switch(t) { case1:break; case2:gotosp1; case0:exit(0); } }}doubleRAD(doubled,doublef,doublem){ doublee; doublesign=(d<:; if(d==0) { sign=(f<:; if(f==0) { sign=(m<:; } } if(d<0) d=d*; if(f<0) f=f*; if(m<0) m=m*; e=sign*(d*3600+f*60+m)*PI/(3600*180); returne;}voidRBD(doublehd){ intt; intd,f; doublem; doublesign=(hd<:; if(hd<0) hd=fabs(hd); hd=hd*3600*180/PI;t=int(hd/3600); d=sign*t; hd=hd-t*3600; f=int(hd/60); m=hd-f*60; printf("%d'%d'%lf'\n",d,f,m);}voidBLH_XYZ(){ doubleB,L,H,N,W; doubled,f,m; doubleX,Y,Z; printf("请输入大地坐标（输入格式为角度(例如：30'40'50')）：\n"); printf("大地经度L="); scanf("%lf'%lf'%lf'",&d,&f,&m); L=RAD(d,f,m); printf("大地纬度B="); scanf("%lf'%lf'%lf'",&d,&f,&m); B=RAD(d,f,m); printf("大地高H="); scanf("%lf",&H); W=sqrt(1-e2*sin(B)*sin(B)); N=a/W; X=(N+H)*cos(B)*cos(L); Y=(N+H)*cos(B)*sin(L); Z=(N*(1-e2)+H)*sin(B); printf("\n\n转换后得到大地空间直角坐标为：\n\n"); printf("X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z);}voidXYZ_BLH(){ doubleB,L,H,N,W; doubleX,Y,Z; doubletgB0,tgB1; printf("请输入大地空间直角坐标：\n"); printf("X="); scanf("%lf",&X); printf("Y="); scanf("%lf",&Y); printf("Z="); scanf("%lf",&Z); printf("\n\n转换后得到大地坐标为：\n\n"); L=atan(Y/X); printf("大地经度为：L="); RBD(L); printf("\n"); tgB0=Z/sqrt(X*X+Y*Y); tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0)); while(fabs(tgB0-tgB1)>5*pow(10,-10)) { tgB0=tgB1; tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0)); } B=atan(tgB1); printf("大地纬度为：B=");RBD(B);printf("\n"); W=sqrt(1-e2*sin(B)*sin(B)); N=a/W; H=sqrt(X*X+Y*Y)/cos(B)-N; printf("大地高为：H=%lf\n\n",H);}voidB_ZS(){ doubleL1,B1,A1,s,d,f,mi; doubleu1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp; doublesgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2; printf("请输入已知点的大地坐标（输入格式为角度(例如：30'40'50')，下同）：\nL1="); scanf("%lf'%lf'%lf'",&d,&f,&mi); L1=RAD(d,f,mi); printf("\nB1="); scanf("%lf'%lf'%lf'",&d,&f,&mi); B1=RAD(d,f,mi); printf("请输入大地方位角：\nA1="); scanf("%lf'%lf'%lf'",&d,&f,&mi); A1=RAD(d,f,mi); printf("请输入该点至另一点的大地线长：\ns="); scanf("%lf",&s); u1=atan(sqrt(1-e2)*tan(B1)); m=asin(cos(u1)*sin(A1)); M=atan(tan(u1)/cos(A1)); m=(m>0)m:m+2*PI; M=(M>0)M:M+PI; k2=ep2*cos(m)*cos(m); alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b; bt=k2/4-k2*k2/8+37*k2*k2*k2/512; r=k2*k2/128-k2*k2*k2/128; sgm0=alfa*s; sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0); while(fabs(sgm0-sgm1)>*PI/180*pow(10,-7)) { sgm0=sgm1; sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0); } sgm0=sgm1; A2=atan(tan(m)/cos(M+sgm0)); A2=(A2>0)A2:A2+PI; A2=(A1>PI)A2:A2+PI; u2=atan(-cos(A2)*tan(M+sgm0)); lmd1=atan(sin(u1)*tan(A1)); lmd1=(lmd1>0)lmd1:lmd1+PI; lmd1=(m<PI)lmd1:lmd1+PI; lmd2=atan(sin(u2)*tan(A2)); lmd2=(lmd2>0)lmd2:lmd2+PI; lmd2=(m<PI)(((M+sgm0)<PI)lmd2:lmd2+PI):(((M+sgm0)>PI)lmd2:lmd2+PI); lmd=lmd2-lmd1; B2=atan(sqrt(1+ep2)*tan(u2)); kp2=e2*cos(m)*cos(m); alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128; btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32; rp=e2*kp2*kp2/256; l=lmd-sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0)+rp*sin(2*sgm0)*cos(4*M+2*sgm0)); L2=L1+l; printf("\n\n得到另一点的大地坐标和大地线在该点的大地方位角为：\n\n"); printf("L2="); RBD(L2);printf("\n"); printf("B2="); RBD(B2);printf("\n"); printf("A2="); RBD(A2); printf("\n");}voidB_FS(){ doubleL1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M; doublel,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2; printf("请输入第一个点大地坐标（输入格式为角度(例如：30'40'50')，下同）：\n大地经度L1="); scanf("%lf'%lf'%lf'",&du,&f,&mi); L1=RAD(du,f,mi); printf("大地纬度B1="); scanf("%lf'%lf'%lf'",&du,&f,&mi); B1=RAD(du,f,mi); printf("\n请输入第二个点大地坐标：\n大地经度：L2="); scanf("%lf'%lf'%lf'",&du,&f,&mi); L2=RAD(du,f,mi); printf("大地纬度：B2="); scanf("%lf'%lf'%lf'",&du,&f,&mi); B2=RAD(du,f,mi); l=L2-L1; u1=atan(sqrt(1-e2)*tan(B1)); u2=atan(sqrt(1-e2)*tan(B2)); sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l)); m0=asin(cos(u1)*cos(u2)*sin(l)/sin(sgm0)); dit_lmd=*sgm0*sin(m0); lmd0=l+dit_lmd; dit_sgm=sin(m0)*dit_lmd; sgm1=sgm0+dit_sgm; m=asin(cos(u1)*cos(u2)*sin(lmd0)/sin(sgm1)); A1=atan(sin(lmd0)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd0))); A1=(A1>0)A1:A1+PI; A1=(m>0)A1:A1+PI; M=atan(sin(u1)*tan(A1)/sin(m)); M=(M>0)M:M+PI; k2=ep2*cos(m)*cos(m); alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b; bt=k2/4-k2*k2/8+37*k2*k2*k2/512; r=k2*k2/128-k2*k2*k2/128; kp2=e2*cos(m)*cos(m); alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128; btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32; rp=e2*kp2*kp2/256; sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l)); sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0)))); while(fabs(sgm0-sgm1)>1*PI/180*pow(10,-8)) { sgm0=sgm1; sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0)))); } sgm=sgm1; lmd=l+sin(m)*(alfap*sgm+btp*sin(sgm)*cos(2*M+sgm)); s=(sgm-bt*sin(sgm)*cos(2*M+sgm)-r*sin(2*sgm)*cos(4*M+2*sgm))/alfa; A1=atan(sin(lmd)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd))); A1=(A1>0)A1:A1+PI; A1=(m>0)A1:A1+PI; A2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(u1)*cos(u2))); A2=(A2>0)A2:A2+PI; A2=(m<0)A2:A2+PI; printf("\n\n得到两点间大地线长S和大地正反方位角A1、A2如下：\n\n"); printf("s=%lf\n",s); printf("A1=");RBD(A1);printf("\n"); printf("A2=");RBD(A2);printf("\n");}voidGUS_ZS(){ doubleB,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t; doubleAt,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06; intDH3,DH6; printf("请输入大地坐标（输入格式为角度(例如：30'40'50')）：\n大地经度L="); scanf("%lf'%lf'%lf'",&du,&f,&mi); L=RAD(du,f,mi); printf("\n大地纬度B="); scanf("%lf'%lf'%lf'",&du,&f,&mi); B=RAD(du,f,mi); At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256; Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512; Ct=15*e2*e2/64+105*e2*e2*e2/256; Dt=35*e2*e2*e2/512; X=a*(1-e2)*(At*B-Bt*sin(2*B)/2+Ct*sin(4*B)/4-Dt*sin(6*B)/6); W=sqrt(1-e2*sin(B)*sin(B)); N=a/W; n=sqrt(ep2)*cos(B); t=tan(B); DH3=(L-*PI/180))/(3*PI/180)+1; DH6=L/(6*PI/180)+1; L03=DH3*(3*PI/180); L06=DH6*(6*PI/180)-(3*PI/180); l3=L-L03; l6=L-L06; m3=cos(B)*l3; m6=cos(B)*l6; x3=X+N*t*(m3*m3/2+(5-t*t+9*n*n+4*n*n*n*n)*m3*m3*m3*m3/24+(61-58*t*t+t*t*t*t)*m3*m3*m3*m3*m3*m3/720); x6=X+N*t*(m6*m6/2+(5-t*t+9*n*n+4*n*n*n*n)*m6*m6*m6*m6/24+(61-58*t*t+t*t*t*t)*m6*m6*m6*m6*m6*m6/720); y3=N*(m3+(1-t*t+n*n)*m3*m3*m3/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m3*m3*m3*m3*m3/120); y6=N*(m6+(1-t*t+n*n)*m6*m6*m6/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m6*m6*m6*m6*m6/120); Y3=DH3*1000000+500000+y3; Y6=DH6*1000000+500000+y6; printf("\n\n得到的高斯平面坐标为:\n\n"); printf("对于3度带:\n纵坐标x=%.3lf\n横坐标y=%.3lf(通用坐标Y=%.3lf)\n\n",x3,y3,Y3); printf("对于6度带:\n纵坐标x=%.3lf\n横坐标y=%.3lf(通用坐标Y=%.3lf)\n\n",x6,y6,Y6);}voidGUS_FS(){ doublex,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6; longDH; printf("请输入高斯平面坐标：\n\n"); printf("纵坐标X="); scanf("%lf",&x);printf("\n"); printf("自然坐标y="); scanf("%lf",&y);printf("\n"); printf("通用坐标Y="); scanf("%lf",&Y);printf("\n"); At=1+3*e