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

#include"stdio.h"#include"math.h"#include"stdlib.h"#include"iostream"#definePI3.1415926535897323doublea,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("请选择椭球 参数 转速和进给参数表a氧化沟运行参数高温蒸汽处理医疗废物pid参数自整定算法口腔医院集中消毒供应 (输入椭球序号):\n");printf("1.克拉索夫斯基椭球参数\n");printf("2.IUGG_1975椭球参数\n");printf("3.CGCS_2000椭球参数\n");printf("0.其他椭球参数(自行输入)\n");scanf("%d",&n);switch(n){case1:a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.00673852541468;break;case2:a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.00673950181947;break;case3:a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.00673949677547;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<0.0)?-1.0:1.0;if(d==0){sign=(f<0.0)?-1.0:1.0;if(f==0){sign=(m<0.0)?-1.0:1.0;}}if(d<0)d=d*(-1.0);if(f<0)f=f*(-1.0);if(m<0)m=m*(-1.0);e=sign*(d*3600+f*60+m)*PI/(3600*180);returne;}voidRBD(doublehd){intt;intd,f;doublem;doublesign=(hd<0.0)?-1.0:1.0;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("请输入大地坐标（输入 格式 pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载 为角度(例如：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)>2.8*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=(m0)?lmd2:lmd2+PI;lmd2=(mPI)?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=0.003351831*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-(1.5*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*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;B0=x/(a*(1-e2)*At);B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);while(fabs(B1-B0)>1*pow(10,-8)){B0=B1;B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);}Bf=B1;nf=sqrt(ep2)*cos(Bf);tf=tan(Bf);Vf=sqrt(1+ep2*cos(Bf)*cos(Bf));Nf=c/Vf;B=Bf-Vf*Vf*tf/2*((y/Nf)*(y/Nf)-(5+3*tf*tf+nf*nf-9*nf*nf*tf*tf)*pow((y/Nf),4)/12+(61+90*tf*tf+45*tf*tf)*pow((y/Nf),6)/360);L=((y/Nf)-(1+2*tf*tf+nf*nf)*(y/Nf)*(y/Nf)*(y/Nf)/6+(5+28*tf*tf+24*pow(tf,4)+6*nf*nf+8*nf*nf*tf*tf)*pow((y/Nf),5)/120)/cos(Bf);DH=Y/1000000;L3=3*PI/180*double(DH)+L;L6=6*PI/180*double(DH)-3*PI/180+L;printf("\n\n得到的大地坐标为:\n\n");printf("大地纬度B=");RBD(B);printf("\n");printf("若为6度带，大地经度L=");RBD(L6);printf("\n");printf("若为3度带，大地经度L=");RBD(L3);printf("\n");}