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

1、大地坐标与空间直角坐标的转换程序代码#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(

2、); 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_1

3、975椭球参数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

4、=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();b

5、reak; 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 s

6、p1; 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; do

7、uble 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

8、,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)

9、; 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=a

10、lfa*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(

11、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

12、*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; do

13、uble 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,&

14、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_

15、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

16、/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)*(alf

17、ap*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

18、*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(B

19、f); 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(nn 得到的大地坐标为:nn); printf(