大地测量学基础实验报告

《大地测量学基础实验报告》由会员分享，可在线阅读，更多相关《大地测量学基础实验报告（40页珍藏版）》请在装配图网上搜索。

1、EAST CHINA INSTITUTE OF TECHNOLOGY大地测量学基础实习报告学号：姓名：班级： 专业：课程名称：大地测量学基础指导老师：2014年5月目录1.引言 31.1.实习的基本任务31.2 实验要求32.实习内容42.1 实习环境42.2 实习主要内容介绍42.2.1 空间大地坐标转化为空间直角坐标42.2.1.1 计算公式52.2.2 空间坐标转化为空间大地坐标52.2.2.1 计算公式62.2.3 高斯投影正算62.2.3.1 计算公式62.2.4 高斯投影反算72.2.4.1 计算公式73.1 空间大地坐标转化为空间直角坐标83.1.1 基本程序83.1.2 程序改

2、进83.1.3 程序完善93.1.4 程序扩展113.2 空间坐标转化为空间大地坐标113.2.1 基本程序113.2.2 程序改进113.2.3 程序完善123.2.4 程序扩展153.3 高斯投影正算153.3.1 基本程序153.3.2 程序改进153.3.3 程序完善163.3.4 程序扩展203.4 高斯投影反算203.4.1 基本程序203.4.2 程序改进213.4.3 程序完善233.4.4 程序扩展264.实习成果264.1 空间大地坐标转化为空间直角坐标264.1.1 初始调试及成功调试274.1.2 运行界面及运行成果284.2 空间直角坐标转化为空间大地坐标284.2.

3、1 初始调试及成功调试294.2.2 运行界面及运行成果294.3 高斯投影正算304.3.1 初始调试及成功调试304.3.2 运行界面及运行成果314.4 高斯投影反算314.4.1 初始调试及成功调试324.4.2 运行界面及运行成果325.编程扩展325.1 儒略日的计算公式335.2 儒略日正算335.2.1 程序代码335.2.2 初始调试及成功调试345.2.3 运行界面及运行成果345.3 儒略日反算355.3.1 程序代码355.3.2 初始调试及成功调试365.3.3 运行界面及运行成果366.总结 371.引言大地测量学基础的程序设计实验是教学中的重要环节，是测 绘工程专

4、业的学生必修的内容。编程计算在大地测量学基础的教 学中不可或缺，是大地测量技术理论联系实际的重要环节，可以培养 学生的动手能力，增强学生对大地测量数据处理在程序设计、算法阐 述，程序实现等方面的认识。它在培养严谨的治学态度、活跃的创新 意识、理论联系实际和适应科技发展的综合应用能力等方面具有其他 实践类课程不可替代的作用。1.1.实习的基本任务培养基本程序开发技能，提高自身的动手能力，加深对空间大地 坐标与空间坐标转化、高斯投影正算与反算等内容的理解。培养自身 的科学思维和创新意识，使自己掌握实验研究的基本方法，提高自身 的分析能力和创新能力。提高自身的科学素养，培养自己理论联系实际和实事求是

5、的科学 作风，认真严谨的科学态度，积极主动的探索精神，遵守纪律，团结 协作，爱护公共财产的优良品德。1.2 实验要求1）分析实习题目的要求。2）写出详细设计说明。（3）编写程序代码，调试程序使其能正确运行。（4）设计完成的软件要便于操作和使用。（5）设计完成后提交实习报告。2.实习内容2.1 实习环境本程序可以在 98/2000/XP 下运行,可以用 VC+6.0 执行。2.2 实习主要内容介绍2.2.1 空间大地坐标转化为空间直角坐标空间大地坐标以参考椭球的中心为坐标原点，椭球的短轴与参考 椭球旋转轴重合，以过地面点的椭球法线与椭球赤道面的夹角为大地 纬度B,以过地面点的椭球子午面与起始子午

6、面之间的夹角为大地经 度L，地面点沿椭球法线至椭球面的距离为大地高H,地面点的点位 用（B, L，H）表示。空间直角坐标以参心0为坐标原点，X轴与起 始子午面和赤道的交线重合，Y轴在赤道面上与X轴垂直，构成右手 直角坐标系0-XYZ，Z轴与参考椭球的短轴（旋转轴）相重合，地面 点的点位用（X，Y，Z）表示。在已知椭球的长半轴a、椭球的第一 偏心率e，空间大地坐标B、L、H，可以求出空间直角坐标X、Y、Z2.2.1.1 计算公式对应计算公式如下：X = (N + H) * cos B * cos LY = (N + H) * cos B * sin L Z 二N*(1 -e2) + H*sinB

7、_公式中，N为椭球面卯酉圈的曲率半径，e为椭球的第一偏心率, a、b椭球的长短半径，f椭球扁率，W为第一辅助系数W = (1 一 e 2 sin 2 B2.2.2 空间坐标转化为空间大地坐标空间直角坐标以参心0为坐标原点，X轴与起始子午面和赤道的 交线重合，Y轴在赤道面上与X轴垂直，构成右手直角坐标系0-XYZ, Z轴与参考椭球的短轴(旋转轴)相重合，地面点的点位用(X, Y， Z)表示。空间大地坐标以参考椭球的中心为坐标原点，椭球的短轴 与参考椭球旋转轴重合，以过地面点的椭球法线与椭球赤道面的夹角 为大地纬度B,以过地面点的椭球子午面与起始子午面之间的夹角为 大地经度L，地面点沿椭球法线至椭

8、球面的距离为大地高H,地面点 的点位用(B, L，H)表示。在已知椭球的长半轴a、椭球的第一偏 心率e,空间直角坐标X、Y、Z,可以求出空间大地坐标B、L、H。2.2.2.1 计算公式对应计算公式如下：YL = arctan()Z *( N + H)X*(1 - e 2) + HlB = arctan(；(X 2 + Y2)*X2 +Y2H = NcosB公式中，N为椭球面卯酉圈的曲率半径，e为椭球的第一偏心率，a、b椭球的长短半径，f椭球扁率，W为第一辅助系数xa2 一b2e =e = !ZH a 或fW =、(1 一 e 2 sin 2 BaN =W2.2.3 高斯投影正算高斯投影正算，已

9、知大地坐标（B,L）及中央子午线经度L0，计算高斯平面坐标（x,y）。2.2.3.1 计算公式x = X + sin B cos Bl2 +sin B cos3 B(5 -12 + 9n 2)14224NNy = N cos Bl +cos3 B(1 -12 +n 2)13 + cos5 B(5 一 18t2 +14 )156 120公式中，所用参数如下：aaaa一亠 sin 2 B H4 sin 4 B-6 sin 6 B + sin 8B2468t = tan B耳 2 = e2 C0S2 B2.2.4 高斯投影反算高斯投影反算，已知高斯平面直角坐标（X,y）及指定中央子午线经度L0,计算

10、大地坐标（B,L）。2.2.4.1 计算公式t3B = B -f2 MNf(5 + 3t 22 - 9q 212)y 4l =N cos B6 N 3 cos B(1 + 2t2 +n 2)y3 +f f120N5 cos B(5 + 28t i 2 + 24t 4 ) y 524M N 3f f f f公式中，所用参数如下：aaaa02X = a B - sin 2B + - sin 4B - 6 sin 6B + 8 sin 8 B4aN =Wt = tan B耳 N = a(1 - e2B )-2 = e2 cos2 B_3M = a(1 - e2)(1 - e2 sln2B )-23.

11、1 空间大地坐标转化为空间直角坐标3.1.1 基本程序依据空间大地坐标转化为空间直角坐标的公式，计算X、Y、Z,其中椭球选择克拉索夫斯基椭球体，程序如下：a=6378245;e2=0.006693421622966;N=a/sqrt(1-e2*sin(radianB)*sin(radianB);X=(N+H)*cos(radianB)*cos(radianL);Y=(N+H)*cos(radianB)*sin(radianL);Z=(N*(1-e2)+H)*sin(radianB);printf(X=%lfnY=%lfnZ=%lfn,X,Y,Z);3.1.2 程序改进空间大地坐标转化为空间直角

12、坐标的过程中椭球类型有4种：克 拉索夫斯基椭球体、1975年国际椭球体、WGS-84椭球体、2000中国 大地坐标系，通过switch语句进行选择，程序如下：printf(Your choose is n1(克拉索夫斯基椭球体)n2(1975年国际椭球 体)n3(WGS-84椭球体)n4(2000中国大地坐标系)：n)；scanf(%d,&tuoqiuti);switch (tuoqiuti)case 1：a=6378245；e2=0.006693421622966；break；case 2：a=6378140；e2=0.006694384999588；break；case 3：a=63781

13、37；e2=0.006694379990130；break；case 4：a=6378137；e2=0.006694380022900；break；N=a/sqrt(1-e2*sin(radianB)*sin(radianB);X=(N+H)*cos(radianB)*cos(radianL);Y=(N+H)*cos(radianB)*sin(radianL);Z=(N*(1-e2)+H)*sin(radianB); printf(X=%lfnY=%lfnZ=%lfn,X,Y,Z);3.1.3 程序完善选择椭球类型时可以使用通过switch语句进行选择，程序如下:# include# incl

14、ude# define PI 3.14159265358979323846void main ()int D1,M1,D2,M2,tuoqiuti,a;double X,Y,Z,e2;double MS1,S1,angleB,radianB,degree1;double MS2,S2,angleL,radianL,degree2;double N,H; printf(Please input angleB,angleL,Hn); scanf(%lf%lf%lf,&angleB,&angleL,&H);/* 纬度转化为弧度 * */D1=(int)(angleB+0.0001);MS1=angl

15、eB-D1;M1=(int)(MS1)*100+0.0001);S1=(MS1*100-M1)*100;degree1=D1+M1/60.0+S1/3600.0;radianB=degree1*PI/180.0;/* 经度转化为弧度 * */D2=(int)(angleL+0.0001);MS2=angleL-D2;M2=(int)(MS2+0.0001)*100;S2=(MS2*100-M2)*100;degree2=D2+M2/60.0+S2/3600.0;radianL=degree2*PI/180.0;printf(radianB=%lfnradianL=%lfn,radianB,ra

16、dianL);printf(Your choose is n1 克拉索夫斯基椭球体)n2(1975年国际椭 球体)n3(WGS-84椭球体)n4(2000中国大地坐标系)：n)；/* P100 表4-1 * */scanf(%d,&tuoqiuti)；switch (tuoqiuti)case 1：a=6378245；e2=0.006693421622966；break；case 2：a=6378140；e2=0.006694384999588；break；case 3：a=6378137；e2=0.006694379990130；break；case 4：a=6378137；e2=0.006

17、694380022900；break；N=a/sqrt(1-e2*sin(radianB)*sin(radianB)；/* P105 公式 4-30 * */ X=(N+H)*cos(radianB)*cos(radianL)； Y=(N+H)*cos(radianB)*sin(radianL)；Z=(N*(1-e2)+H)*sin(radianB); printf(X=%lfnY=%lfnZ=%lfn,X,Y,Z);3.1.4 程序扩展选择椭球类型时可以使用到判断、循环等语句，进行选择、重复运算。3.2 空间坐标转化为空间大地坐标3.2.1 基本程序依据空间坐标转化为空间大地坐标的公式，计算

18、B、L、H,程序如下：L=atan(Y/X);N=a/sqrt(1-e2*sin(B)*sin(B);H=(sqrt(X*X+Y*Y)/cos(B)-N;B=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(B)/sqrt(X*X+Y*Y);printf(B=%lfnL=%lfnH=%lfn,B, L,H);3.2.2 程序改进依据空间坐标转化为空间大地坐标的公式，计算B、L、H,此程序关键在于迭代计算(如黑色加出部分)程序如下：B=atan(Z/sqrt(X*X+Y*Y);printf(B=%lfn,B);while (fabs(B-B0)ERROR)B0=B;N=a/sqrt(1

19、-e2*sin(B0)*sin(B0);H=(sqrt(X*X+Y*Y)/cos(B0)-N;B=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(B0)/sqrt(X*X+Y*Y);Bf=B;printf(Bf=%lfn,Bf);H=(sqrt(X*X+Y*Y)/cos(Bf)-N;radianB=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(Bf)/sqrt(X*X+Y*Y);degreeB=radianB*180/PI;D1=(int)(degreeB+0.000001);MS1=degreeB-D;M1=(int)(MS1*60);s1=(MS1*60-M1)*

20、60;angleB=D1+M1/100.0+s1/10000.0;printf(angleB=%lfnangleL=%lfnH=%lfn,angleB,angleL,H);3.2.3 程序完善总述以上基本程序和迭代程序，可得如下程序：# # include# include# define PI 3.14159265358979323846# define ERROR le-10/ leTO 为 1* 1O(TO)void main ()int tuoqiuti,a;int D,M;int D1,M1;double X,Y,Z,e2;double B,B0,Bf;double H,N;doub

21、le radianL,degreeL,angleL,MS,s;double radianB,degreeB,angleB,MS1,s1;printf(Please input X,Y,Zn);scanf(%lf%lf%lf,&X,&Y,&Z);printf(Your choose is nl(克拉索夫斯基椭球体)n2(1975年国际椭球 体)n3(WGS-84椭球体)n4(2000中国大地坐标系)：n)；scanf(%d,&tuoqiuti);/* P100 表4-1 */switch (tuoqiuti)case 1：a=6378245；e2=0.006693421622966；break；

22、case 2：a=6378140；e2=0.006694384999588；break；case 3：a=6378137；e2=0.006694379990130；break；case 4：a=6378137；e2=0.006694380022900；break；/* P105 公式4-31 */radianL=atan(Y/X)；degreeL=radianL*180/PI；D=(int)(degreeL+0.000001)；MS=degreeL-D；M=(int)(MS*60)；s=(MS*60-M)*60；angleL=D+M/100.0+s/10000.0；printf(angleL=

23、%lfn,angleL);B=atan(Z/sqrt(X*X+Y*Y);printf(B=%lfn,B); /* 迭代计算-Bf 公式4-32 */while (fabs(B-B0)ERROR)B0=B;N=a/sqrt(1-e2*sin(B0)*sin(B0);H=(sqrt(X*X+Y*Y)/cos(B0)-N;B=atan(Z/sqrt(X*X+Y*Y)+N*e2*sin(B0)/sqrt(X*X+Y*Y);Bf=B;printf(Bf=%lfn,Bf);/* * P105 公式4-32/4-35 * */H=(sqrt(X*X+Y*Y)/cos(Bf)-N; radianB=atan(

24、Z/sqrt(X*X+Y*Y)+N*e2*sin(Bf)/sqrt(X*X+Y*Y); degreeB=radianB*180/PI;D1=(int)(degreeB+0.000001);MS1=degreeB-D;M1=(int)(MS1*60); s1=(MS1*60-M1)*60;angleB=D1+M1/100.0+s1/10000.0; printf(angleB=%lfnangleL=%lfnH=%lfn,angleB,angleL,H);3.2.4 程序扩展此程序的迭代运算中，常常要使用到判断、循环等语句，进行选择、重复运算，可加以扩展。3.3高斯投影正算3.3.1 基本程序根据

25、高斯投影正算公式，由大地坐标(B,L )及中央子午线经度L0,计算高斯平面坐标(x,y),选择椭球类型为WGS-84椭球体，程序如下：N=a/sqrt(1-e2*sin(radianB)*sin(radianB);t=tan(radianB);n=e22*cos(radianB);x=X+N*sin(radianB)*cos(radianB)*(radianj*radianj)/2+N*sin(radianB)*co s(radianB)*cos(radianB)*cos(radianB)*(5-t*t+9*n*n)*(radianj*radianj*r adianj*radianj)/24;

26、y=N*cos(radianB)*(radianj)+N*cos(radianB)*cos(radianB)*cos(radianB)*(1-t*t+n*n)*(radianj*radianj*radianj)/6+N*cos(radianB)*cos(radianB)*co s(radianB)*cos(radianB)*cos(radianB)*(5-18*t*t+n*n*n*n)*(radianj*radi anj*radianj*radianj*radianj)/120;printf(x=%lfny=%lfn,x,y);3.3.2 程序改进高斯投影正算中，椭球类型有4种：克拉索夫斯基椭

27、球体、1975年国际椭球体、WGS-84椭球体、2000中国大地坐标系，通过switch语句进行选择，程序如下：printf(Your choose is nl(克拉索夫斯基椭球体)n2(1975年国际椭球体)n3(WGS-84椭球体)n4(2000中国大地坐标系)：n)；scanf(%d,&tuoqiuti);switch (tuoqiuti)case1：a=6378245；e2=0.006693421622966；e22=0.006738525414683；break；case2：a=6378140；e2=0.006694384999588；e22=0.006739501819473；br

28、eak；case3：a=6378137；e2=0.006694379990130；e22=0.006739496742270；break；case4：a=6378137；e2=0.006694380022900；e22=0.006739496775480；break；N=a/sqrt(1-e2*sin(radianB)*sin(radianB)；t=tan(radianB)；n=e22*cos(radianB)；x=X+N*sin(radianB)*cos(radianB)*(radianj*radianj)/2+N*sin(radianB)*co s(radianB)*cos(radianB

29、)*cos(radianB)*(5-t*t+9*n*n)*(radianj*radianj*r adianj*radianj)/24；y=N*cos(radianB)*(radianj)+N*cos(radianB)*cos(radianB)*cos(radianB)*(1-t*t+n*n)*(radianj*radianj*radianj)/6+N*cos(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*(5-18*t*t+n*n*n*n)*(radianj*radi anj*radianj*radianj*radia

30、nj)/120；printf(x=%lfny=%lfn,x,y)；3.3.3 程序完善整合以上程序可得如下结果：# include# include# define PI 3.14159265358979323846void main ()int D1,M1,D2,M2,tuoqiuti,a,n1;double X,e2,e22;double MS1,S1,angleB,radianB,degree1;double MS2,S2,angleL,radianL,degree2,l1,l2;double radianj;double N;double m0,m2,m4,m6,m8;double a

31、0,a2,a4,a6,a8;double x,y;double t,n;printf(Please input angleB,angleLn);scanf(%lf%lf,&angleB,&angleL); /* * 纬度转化弧度 * */D1=(int)(angleB+0.0001);MS1=angleB-D1;M1=(int)(MS1)*100+0.0001);S1=(MS1*100-M1)*100;degree1=D1+M1/60.0+S1/3600.0;radianB=degree1*PI/180.0; /* * 经度转化弧度 * */D2=(int)(angleL+0.0001);MS

32、2=angleL-D2;M2=(int)(MS2+0.0001)*100;S2=(MS2*100-M2)*100;degree2=D2+M2/60.0+S2/3600.0; radianL=degree2*PI/180.0;/* * 计算经差 * */ n1=(int)(degree2+6)/6+0.0001);l1=6*n1-3;l2=(degree2-l1); radianj=l2*PI/180.0;printf(radianB=%lfnradianL=%lfnradianj=%lfn,radianB,radia nL,radianj);/输出弧度 B,L,Jprintf(Your cho

33、ose is n1 克拉索夫斯基椭球体)n2(1975年国际椭 球体)n3(WGS-84椭球体)n4(2000中国大地坐标系)：n)；scanf(%d,&tuoqiuti); /* * P100 表4-1 * */ switch (tuoqiuti)case1：a=6378245;e2=0.006693421622966;e22=0.006738525414683;break;case2：a=6378140;e2=0.006694384999588;e22=0.006739501819473;break;case3：a=6378137;e2=0.006694379990130;e22=0.00

34、6739496742270;break;case4：a=6378137;e2=0.006694380022900;e22=0.006739496775480;break; /* * P111 公式4-72 * */m0=a*(1-e2);m2=(3*e2*m0)/2.0;m4=(5*e2*m2)/4.0;m6=(7*e2*m4)/6.0;m8=(9*e2*m6)/8.0; /* P115 公式4-100 * */a0=m0+m2/2.0+3*m4/8.0+5*m6/16.0+35*m8/128.0;a2=m2/2.0+m4/2.0+15*m6/32.0+7*m8/16.0;a4=m4/8.0+

35、3*m6/16.0+7*m8/32.0;a6=m6/32.0+m8/16.0;a8=m8/128.0;X=a0*radianB-a2*sin(2*radianB)/2+a4*sin(4*radianB)/4-a6*sin(6*rad ianB)/6+a8*sin(8*radianB)/8;/* P100 公式 4-5 和 P104 公式 4-19*/N=a/sqrt(1-e2*sin(radianB)*sin(radianB);t=tan(radianB);n=e22*cos(radianB); /* P168 公式4-366 * */x=X+N*sin(radianB)*cos(radian

36、B)*(radianj*radianj)/2+N*sin(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*(5-t*t+9*n*n)*(radianj*radianj*r adianj*radianj)/24;/168 y=N*cos(radianB)*(radianj)+N*cos(radianB)*cos(radianB)*cos(radianB)*(1-t*t+n*n)*(radianj*radianj*radianj)/6+N*cos(radianB)*cos(radianB)*cos(radianB)*cos(radianB)*cos(

37、radianB)*(5-18*t*t+n*n*n*n)*(radianj*radi anj*radianj*radianj*radianj)/120;printf(x=%lfny=%lfn,x,y);3.3.4 程序扩展高斯投影正算要使用到判断、循环等语句，进行选择、重复运算，可以加以扩展。3.4高斯投影反算3.4.1 基本程序依据高斯投影反算的公式，由高斯平面直角坐标(x,y)及指定中央子午线经度L0，计算大地坐标(B,L)，基本程序如下：B0=x/a0;tp1=-a2 *sin(2.0 *B0)/2.0+a4 *sin(4.0 *B0)/4.0;tp2=-a6*sin(6.0*B0)/6.

38、0+a8*sin(8.0*B0)/8.0;B=(x-tp1-tp2)/a0;printf(a0=%lfnB0=%lfnB=%lfn,a0,B0,B);while(fabs(B-B0)ERROR)B0=B;tp1=-a2*sin(2.0*B0)/2.0+a4*sin(4.0*B0)/4.0;tp2=-a6*sin(6.0*B0)/6.0+a8*sin(8.0*B0)/8.0;B=(x-tp1-tp2)/a0;Bf=B;tf=tan(Bf);n=e22*cos(Bf)*cos(Bf);Mf=a*(1-e2)/sqrt(1-e2*sin(Bf)*sin(Bf)*(1-e2*sin(Bf)*sin(B

39、f)*(1 -e2*sin(Bf)*sin(Bf);Nf=a/sqrt(1-e2*sin(Bf)*sin(Bf);B=Bf-tf*y*y/(2*Mf*Nf)+tf*tf*tf*y*y*y*y*(5+3*tf*tf+n-9*n*tf*tf)/(2 4*Mf*Nf*Nf*Nf);degreeB=B*180/PI;D=(int)(degreeB+0.000001);MS=degreeB-D;M=(int)(MS*60);s=(MS*60-M)*60;angleB=D+M/100.0+s/10000.0;l=y/(Nf*cos(Bf)-(1+2*tf*tf+n)*y*y*y/(6*Nf*Nf*Nf*c

40、os(Bf)+(5+28*t f*tf+24*tf*tf*tf*tf)*y*y*y*y*y/(120*Nf*Nf*Nf*Nf*Nf*cos(Bf);degreeL=l*180/PI;prin tf(请输入带号：n);scanf(%d,&N);L0=6*N-3;angleu=L0+degreeL;printf(angleB=%lfnangleu=%lfn,angleB,angleu);3.4.2 程序改进依据高斯投影正算的公式，计算B、L此程序关键在于迭代计算(如黑色加出部分)程序如下：B0=x/a0; tp1=-a2*sin(2.0*B0)/2.0+a4*sin(4.0*B0)/4.0;tp2

41、=-a6*sin(6.0*B0)/6.0+a8*sin(8.0*B0)/8.0; B=(x-tp1-tp2)/a0;printf(a0=%lfnB0=%lfnB=%lfn,a0,B0,B); while(fabs(B-B0)ERROR)B0=B;tp1=-a2*sin(2.0*B0)/2.0+a4*sin(4.0*B0)/4.0;tp2=-a6*sin(6.0*B0)/6.0+a8*sin(8.0*B0)/8.0;B=(x-tp1-tp2)/a0;Bf=B;tf=tan(Bf);n=e22*cos(Bf)*cos(Bf);Mf=a*(1-e2)/sqrt(1-e2*sin(Bf)*sin(Bf

42、)*(1-e2*sin(Bf)*sin(Bf) *(1-e2*sin(Bf)*sin(Bf);Nf=a/sqrt(1-e2*sin(Bf)*sin(Bf);B=Bf-tf*y*y/(2*Mf*Nf)+tf*tf*tf*y*y*y*y*(5+3*tf*tf+n-9*n*tf*tf) /(24*Mf*Nf*Nf*Nf);degreeB=B*180/PI;D=(int)(degreeB+0.000001);MS=degreeB-D;M=(int)(MS*60);s=(MS*60-M)*60;angleB=D+M/100.0+s/10000.0;l=y/(Nf*cos(Bf)-(1+2*tf*tf+n

43、)*y*y*y/(6*Nf*Nf*Nf*cos(Bf)+(5+2 8*tf*tf+24*tf*tf*tf*tf)*y*y*y*y*y/(120*Nf*Nf*Nf*Nf*Nf*cos(Bf);degreeL=l*180/PI;prin tf(请输入带号：n);scanf(%d,&N);L0=6*N-3;angleu=L0+degreeL; printf(angleB=%lfnangleu=%lfn,angleB,angleu);3.4.3 程序完善总述以上基本程序和迭代程序，可得如下程序：# include# include# define PI 3.14159265358979323846#

44、define ERROR le-10/ leTO 为 1* 1O(TO)void main ()int tuoqiuti,a,N,LO;int D,M;double x,y,e2,e22;double mO,m2,m4,m6,m8;double aO,a2,a4,a6,a8;double BO,B,Bf,l,angleu;double tp1,tp2;double tf,Mf,Nf;double n;double degreeB,MS,s,angleB;double degreeL;printf(Please input x,yn);scanf(%lf%lf,&x,&y);printf(You

45、r choose is n1(克拉索夫斯基椭球体)n2(1975年国际椭 球体)n3(WGS-84椭球体)n4(2000中国大地坐标系)：n)；scanf(%d,&tuoqiuti);switch (tuoqiuti)case1:a=6378245;e2=0.006693421622966;e22=0.006738525414683;break;case2:a=6378140;e2=0.006694384999588;e22=0.006739501819473;break;case3:a=6378137;e2=0.006694379990130;e22=0.006739496742270;br

46、eak;case4:a=6378137;e2=0.006694380022900;e22=0.006739496775480;break;/* Pill 公式 472* */m0=a*(le2);m2=(3*e2*m0)/2.0;m4=(5*e2*m2)/4.0;m6=(7*e2*m4)/6.0;m8=(9*e2*m6)/8.0;/* Pii5 公式 4100*/a0=m0+m2/2.0+3*m4/8.0+5*m6/l6.0+35*m8/l28.0;a2=m2/2.0+m4/2.0+l5*m6/32.0+7*m8/l6.0;a4=m4/8.0+3*m6/l6.0+7*m8/32.0;a6=m6

47、/32.0+m8/l6.0;a8=m8/l28.0;/*P115 公式 4101*/B0=x/a0;tp1=-a2*sin(2.0*B0)/2.0+a4*sin(4.0*B0)/4.0; tp2=-a6*sin(6.0*B0)/6.0+a8*sin(8.0*B0)/8.0;B=(x-tp1-tp2)/a0; printf(a0=%lfnB0=%lfnB=%lfn,a0,B0,B); while(fabs(B-B0)ERROR)B0=B; tp1=-a2*sin(2.0*B0)/2.0+a4*sin(4.0*B0)/4.0; tp2=-a6*sin(6.0*B0)/6.0+a8*sin(8.0*

48、B0)/8.0; B=(x-tp1-tp2)/a0;Bf=B; /* * P100 公式4-4 * */ tf=tan(Bf);n=e22*cos(Bf)*cos(Bf); /* * P110 公式4-68/4-69 * */ Mf=a*(1-e2)/sqrt(1-e2*sin(Bf)*sin(Bf)*(1-e2*sin(Bf)*sin(Bf) *(1-e2*sin(Bf)*sin(Bf);Nf=a/sqrt(1-e2*sin(Bf)*sin(Bf); /* * P171 公式4-383 * */ B=Bf-tf*y*y/(2*Mf*Nf)+tf*tf*tf*y*y*y*y*(5+3*tf*t

49、f+n-9*n*tf*tf) /(24*Mf*Nf*Nf*Nf);degreeB=B*180/PI;D=(int)(degreeB+0.000001);MS=degreeB-D;M=(int)(MS*60);s=(MS*60-M)*60;angleB=D+M/100.0+s/10000.0;l=y/(Nf*cos(Bf)-(1+2*tf*tf+n)*y*y*y/(6*Nf*Nf*Nf*cos(Bf)+(5+28*tf*tf+24*tf*tf*tf*tf)*y*y*y*y*y/(120*Nf*Nf*Nf*Nf*Nf*cos(Bf); degreeL=l*180/PI;prin tf(请输入带号：

50、n);scanf(%d,&N);L0=6*N-3;angleu=L0+degreeL; printf(angleB=%lfnangleu=%lfn,angleB,angleu);3.4.4 程序扩展此程序的迭代运算中，常常要使用到判断、循环等语句，进行选择、重复运算，可加以扩展。4.实习成果4.1 空间大地坐标转化为空间直角坐标例子如下，其中椭球选择克拉索夫斯基椭球体，即：a=6378245;e2=0.006693421622966;BLH空间大地坐标44.5959999945999999.9987XYZ空间直角坐标3694472.4683694472.4685194534.424411初始调

51、试及成功调试g| M ! - Mera us h MiuJ t+ -hAMbi：LH.K*z | 血 | 丄i|0 由 E1叶* QkJ曰 rjtJd Xk idcw H4p当序RM4 CIt JB专輪21l如111电唱*r F l + 业盂鱼W:.m爾 i&y 口 I iEM 门 Bfjmi二|直锯已昌1 H占寸 r- - I 上BfRfMrTbaW场牲翻Z: C-3 .：=fmi.敗些砲J&yrrbdi in kVofkipnIfJ 列IgT IKIBiEW85 ： H lBWWtSWteV 發 apfi#曲Tjnritftl. 發 e4G&p|- anjLi?E -.DE41 J畸 un

52、ite TH 7ffl!| l K I.陋 1-1#4 .OHI J 2 51-i|H5.51fr Irfffl: dtrH-laHI/hl.MMfill.l： ri41cO4gr0 I P I/9III1 .;绘度揖化力 T 龍 Ifc-Kfc-S-irh- :m a LPI -. AMI ) ?lEJ-Mcgl-iL-U;-(!( -：.2 -*. 1 *；degiH? :F/idQn. iiF款Jl屮 W打呛fep|尹谒Q !priiifi FfailB diiE %! rnr dda dill .1 rrini訂 ids diiL j ;pr EM“F CMUS+丄,W诧时玉金判(理卄忖冉Elf湘制诛*询何小啾井休八皿怦0冲附ifei好总讥丁样 F h*ri W*b*Vh Wb*WTb* WkPhC dllKF ;a7-De7*TM9D fl;brrak.;Hua Lrl 斗1 rind jh fiL*n H 皿列上 I SQL D.bunn?nT库同丈城愛碎衍隹冋业闻iLio1矿济旳丁*$*打“F.曲“if吋412运行界面及运行成果- 、C:UsersAdininistratofXDesktop网上脣程D已berg空囘大地坐忘坟化为空间坐忘ELH_XZ (.Please input angleBj-angleL,H 44.59599999 4