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

#includestdio.h#includemath.h#includestdlib.h#includeiostream#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(请选择椭球参数(输入椭球序号):\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=(d0.0)?-1.0:1.0;if(d==0){sign=(f0.0)?-1.0:1.0;if(f==0){sign=(m0.0)?-1.0:1.0;}}if(d0)d=d*(-1.0);if(f0)f=f*(-1.0);if(m0)m=m*(-1.0);e=sign*(d*3600+f*60+m)*PI/(3600*180);returne;}voidRBD(doublehd){intt;intd,f;doublem;doublesign=(hd0.0)?-1.0:1.0;if(hd0)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=(m0)?m:m+2*PI;M=(M0)?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=(A20)?A2:A2+PI;A2=(A1PI)?A2:A2+PI;u2=atan(-cos(A2)*tan(M+sgm0));lmd1=atan(sin(u1)*tan(A1));lmd1=(lmd10)?lmd1:lmd1+PI;lmd1=(mPI)?lmd1:lmd1+PI;lmd2=atan(sin(u2)*tan(A2));lmd2=(lmd20)?lmd2:lmd2+PI;lmd2=(mPI)?(((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