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

#include stdio.h#include math.h#include stdlib.h#include iostream#define PI 3.1415926535897323double a,b,c,e2,ep2;int main()int m,n,t;double RAD(double d,double f,double m);void RBD(double hd);void BLH_XYZ();void XYZ_BLH();void B_ZS();void B_FS();void GUS_ZS();void GUS_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) case 1:a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.00673852541468;break;case 2:a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.00673950181947;break;case 3:a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.00673949677547;break;case 0: 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(nn输入错误!n请重新输入!nn);goto sp2 ;while(1)switch(m) case 1:BLH_XYZ();break;case 2:XYZ_BLH();break;case 3:B_ZS();break;case 4:B_FS();break;case 5:GUS_ZS();break;case 6:GUS_FS();break;default:printf(nn输入错误!n请重新输入!nn);goto sp1 ;printf(是否继续进行此功能计算? nn);printf(（ 若继续进行此功能计算，则输入1；n 若选择其他功能进行计算，则输入2；n 若退出， 则输入0. ）n);scanf(%d,&t);switch(t) case 1:break;case 2:goto sp1;case 0:exit(0);double RAD(double d,double f,double m)double e;double sign=(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);return e;void RBD(double hd)int t;int d,f;double m;double sign=(hd0.0)?-1.0:1.0;if(hd5*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=%lfnn,H);void B_ZS()double L1,B1,A1,s,d,f,mi;double u1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;double sgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2;printf(请输入已知点的大地坐标（输入格式为角度(例如：304050)，下同）：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=(m0)?lmd2:lmd2+PI;lmd2=(mPI)?(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(nn得到另一点的大地坐标和大地线在该点的大地方位角为：nn);printf(L2=);RBD(L2);printf(n);printf(B2=);RBD(B2);printf(n);printf(A2=);RBD(A2);printf(n);void B_FS()double L1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;double l,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2;printf(请输入第一个点大地坐标（输入格式为角度(例如：304050)，下同）：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=(A10)?A1:A1+PI;A1=(m0)?A1:A1+PI;M=atan(sin(u1)*tan(A1)/sin(m);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;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=(A10)?A1:A1+PI;A1=(m0)?A1:A1+PI;A2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(u1)*cos(u2);A2=(A20)?A2:A2+PI;A2=(m1*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+(