武汉大学测绘学院大地测量计算与实习报告

1、课程编号： 课程性质：必修大地测量计算与实习实 习 报 告学 院： 测绘学院 专业方向： 测绘工程A方向 实习地点： 武汉大学 班级组号： X 班第X小组 学生姓名： XXX 学号： XXX 指导教师： 刘宗泉、苏新洲、丁士俊向 东、史俊波 等2013年 6月24日 至 2013年 7月12日目 录目 录.2一 技术报告.41.1.实习的目的和意义.41.2. 实习的任务.4二 二等精密水准测量外业观测与计算 .4 2.1. 测区概况.4 2.2 已知高程.4 2.3 作业依据.5 2.4 水准测量详述.5 2.4.1踏勘、选点.6 2.4.2使用的仪器和仪器检验.6 (1) i角检验.6 (

2、2) 一对水准尺零点差检验. 62.4.3 外业数据采集与概算 7 2.5 外业成果表格. 7三 大地测量计算课程设计.73.1 编程语言简介.73.2 编程任务、程序框图及基本数学模型81）高斯投影坐标正反算.82）实测斜距化算到高斯平面距离.10 3）大地主题正反算（高斯平均引数+白塞尔）.11 3.3 大地测量计算成果14四 实习总结15附录1 数字水准仪i角检验结果表21附录2 数字水准尺的零点差检验结果表22附录 3 水准点点之记.23附录4 自测二等水准外业观测高差与高程概算表.24附录5 大地测量计算源程序.25附录6 国测二等水准外业观测高差与高程概算表.40一、技术报告1.1

3、 实习的目的和意义大地测量计算与实习课程，是学习大地测量学基础之后的后续课程，是一次重要的集中教学实习。 通过这次集中地教学实习，旨在巩固我们在课堂上学过的理论知识，锻炼和提高我们的仪器操作能力和动手能力，增强综合、灵活运用所学知识的能力；在我们自主编程计算的过程中，提高发现问题，分析问题和解决相关问题的能力；开拓我们的视野，把基本理论与实践紧密的结合在一起，对理论知识进一步消化和吸收，激化和培养我们的创新思维的能力，为今后走上工作岗位打下较为扎实的专业基础。1.2 实习的任务本次实习共有两项任务1、二等精密水准测量外业观测与概算。(1) 踏勘和选点：了解水准路线的布设和水准点位置的选择方法，

4、根据安排的水准路线，各组独立完成本组的选点工作。(2) 填写水准点之记：每人完成一个水准点“点之记”的绘制，且每个组的成员间不能同点。(3) 仪器检验：每人完成数字水准仪i角测定成果和一对水准标尺零点差的测定结果各一份。(4) 数字水准观测和记录：每人使用数字水准仪完成至少80站的观测和记簿，每小组完成珞珈山环或樱园环的往返闭合水准环线。(5) 光学水准仪观测记录：每人使用光学精密水准仪完成至少10站的观测和记簿，每组构成闭合环。(6) 外业观测数据检查与概算。2、大地测量编程实习(1) 高斯投影正反算的计算程序。(2) 实测斜距化算至高斯投影平面边长的边长改正程序。(3) 大地主题正反算的计

5、算程序。(4) 水准测量外业高差与概略高程表的编制。二、二等精密水准测量外业观测与计算2.1 测区概况本次我们5班第2小组的测区是樱园环，这段水准线路线，地质为混凝土和柏油马路，部分测段穿过马路，来往车辆较多，教务部至行政楼高差起伏较大，其余路段较为平坦。2.2 已知高程环 名点 名已知高程(m)珞珈山环未名湖119.8882.3 作业依据国家测绘局，国家一、二等水准测量规范2006-05-24 测绘出版社，2010.2010仪器的技术指标 (1) 水准仪的i角限差为15(2) 标尺的零点不等差为0.10mm(3) 测站观测限差仪器类型基辅分划读数的差（mm）基辅分划所测高差的差(mm)检测间

6、歇点高差的差(mm)数字水准仪0.40.61.0光学水准仪0.40.61.0仪器类型视距限值前后视距差(m)累计视距差(m)数字水准仪=3且=501.56光学水准仪=5013(4) 观测程序光学水准仪数字水准仪备注测站程序奇数站后前前后奇数站后前前后左边表格均为往测，光学仪器返测时与往测相反，数字水准仪返测与往测相同偶数站前后后前偶数站前后后前(5) 闭合差闭合差往返高差不符值 4mm( k为测段长度，单位km,0.1km按0.1km算)环闭合差 4mm ( F为环线长度，单位为km )2.4 水准测量详述2.4.1 踏勘、选点樱园环点号位置1未名湖2行政楼3新闻与传播学院4城市设计学院5田园

7、小观园6大自强超市7生科院8图书馆2.4.2 使用的仪器和仪器检验使用的仪器为DL200。进行了下面两项检验。（1）i角检验1.在较平坦的地方选定适当距离（约取40m）的两点A、B，并用尺垫确定不动。2.置水准仪于A 点附近，如图：此时测量A，B尺的读数a1,b1，并记录；3.然后将水准仪置于B 点的右边，在线段AB 的延长线上，如图读取此时A,B尺的读数a2,b2,并记录;4.利用相关公式解算i角大小；（2）一对水准尺零点差检验在距离水准仪约20m30m 的等距离处选择3 个点，使得各个点之间存在着约20cm 的高差。此项检验进行3 个测回，每一测回中，分别在3 个点上一次安置一对标尺，分别

8、用仪器进行读数，测回间应变换仪器高。由于使用的是数字水准仪，因而设置重复测量次数为5 次，每测回每点连续观测4 次。水准尺零点差即两次读数中数之差。检验的结果见附录。2.4.3 外业数据采集与概算水准测量路线图：我观测的是行政楼到新闻与传播学院的往返测，与6月27日下午完成。往测采用手簿记录，返测为电子记录。大约有40站（加上后来补测的一些站数总数达到要求的80站）。往返测闭合差+2.0mm，具体成果见附录4。2.5 外业成果表格成果表格有：1） 数字水准仪i角检验结果表(附录1)。2） 数字水准尺的零点差检验结果表(附录2)。3） 水准点点之记(附录3)。4） 自测测量外业观测高差与高程概算

9、成果表(附录4)。三、大地测量计算课程设计3.1 编程语言简介本次大地测量编程实习，我采用的是C+语言，主要利用了MFC。C+语言是一种优秀的面向对象程序设计语言，它在C语言的基础上发展而来，但它比C语言更容易为人们学习和掌握。C+以其独特的语言机制在计算机科学的各个领域中得到了广泛的应用。面向对象的设计思想是在原来结构化程序设计方法基础上的一个质的飞跃，C+完美地体现了面向对象的各种特性。MFC(Microsoft Foundation Classes)，是一个微软公司提供的类库（class libraries），以C+类的形式封装了Windows的API，并且包含一个应用程序框架，以减少应

10、用程序开发人员的工作量。其中包含的类包含大量Windows句柄封装类和很多Windows的内建控件和组件的封装类。3.2 编程任务、程序框图及基本数学模型本次大地测量实习主要编写三个程序，分别是高斯投影坐标正反算和实测斜距归算至高斯平面和大地主题正反算。1） 高斯投影坐标正反算1. 高斯投影坐标正算输入B，L计算带号N,中央纬度L0，纬度差 lN=L/6(进一取整)，L0=6N，l=L-L0.按照正算公式计算x, y转换为国家统一坐标X , YX=x, Y=N+y输出X，Y高斯投影坐标正算公式 2. 高斯投影坐标反算输入国家统一坐标X，Y确定带号N，算出x, yN=Y/(向下取整) x=X,

11、y=Y-N-求中央子午线L0L0=6N按照反算公式计算B，l求出大地经度LL=L0+l输出B，L高斯投影坐标反算公式 主要代码见附录5.2） 实测斜距归算至高斯平面边长先利用高斯坐标计算大地方位角A，之后利用下面的公式求出实测斜距在椭球面上的大地线长度，S=D*1-(h2-h1)/D*(h2-h1)/D/(1+h1/Ra)*(1+h2/Ra)然后利用距离改化公式，d=S*(1+ym*ym/2/Rm/Rm+(y2-y1)*(y2-y1)/24/Rm/Rm+ym*ym*ym*ym/24/Rm/Rm/Rm/Rm);主要代码见附录5.3） 大地测量主题正反算1. 高斯平均引数正算首先把勒让德级数在P1

12、点展开改在大底线长度中低昂M处展开，使收敛级数减少，收敛快，精度高;其次，考虑到求中点M的复杂性，将M点用大地线两单点平均纬度及平均方位角相对应的m点来代替，并借助迭代计算，便可顺利的实现大地主题正算。2. 高斯平均引数反算由于已知L1、B1及L2、B2，这时经度差、纬度差及平均纬度均为已知，故可以正算公式很容易导出反算公式。3. 白塞尔大地主题正算（1） 计算起点的归化纬度（2） 计算辅助函数值，解球面三角形（3） 按公式计算相关系数A,B,C以及,（4） 计算球面长度（5） 计算纬度差改正数（6） 计算终点大地坐标及大地方位角4.白塞尔大地主题反算（1） 辅助计算（2） 用逐次趋近法同时计

13、算起点大地方位角、球面长度及经差，第一次趋近时，取。计算下式,重复上述计算过程2.（3） 计算大地线长度S（4） 计算反方位角 主要代码见附录5.3.3 大地测量计算成果3.3.1 高斯投影正反算计算结果（采用第2组第2个数据）： 已知数据x = .5697y = .5108反算结果B = 412342.L = 1250624.正算结果x = .5697y = .51083.3.2 实测斜距化算至高斯投影平面边长计算结果（采用克拉索夫斯基椭球）：D12 = 578.8686m。3.3.3 大地主题解算计算结果（采用5班第2组第2个数据）：正算反算已知数据计算结果检核结果B1 = 43.B2 =

14、 44.S=L1 = 132.L2 = 132.A12= 1.A12= 1.A21= 181.A21= 181.S = 3.3.4 国测数据外业高差与概略高程表计算结果见附录6四实习总结附录1 数字水准仪的i角检验结果表i角的检校仪器：DL-201/2007 方法：I1ABI2 观测者：王蔷日期：2013-06-26 标尺： 记录者：黄式春时间：09:00 呈像：清晰稳定 检查者：尹建鹏仪器距近标尺距离6.4 m 仪器距远标尺距离45.6m仪器站观测次序A尺读数B尺读数A尺读数B尺读数1311.8266.8272.8227.32311.6266.8273.0227.33311.6266.727

15、2.9227.44311.6266.7272.9227.4中数311.65266.75272.9227.35高差（a-b）,mm44.945.55方法：， 0.325 mm方法： mm1.710-0.8372=0.87”校正： 附录二 数字水准仪的零点差检验结果表一对水准尺零点差测定标尺： 观测者：陆直仪器： DL-201/2007 记录者：揭元翔日期： 2013-06-29 检查者：尹建鹏桩号No. 标尺读数No. 标尺读数 I测回（mm）II测回（mm）III测回（mm）I测回（mm）II测回（mm）III测回（mm）11365.51406.81253.71365.41406.71253.

16、41365.51406.71253.81365.51406.61253.51365.61406.71253.71365.51406.71253.51365.51406.71253.61365.51406.71253.621052.91094.0940.71052.61093.8940.71052.81093.9940.71052.71093.8940.61052.81093.8940.81052.61093.8940.61052.91093.9940.81052.71093.8940.73789.8830.9677.7789.8830.9677.6789.8830.9677.7789.8830

17、.8677.6789.8830.9677.8789.7830.8677.6789.9830.9677.7789.7830.8677.7平均1069.41110.51957.391069.291110.43957.27总平均1045.771045.67 一对标尺零点不等差=-0.10mm附录3 水准点点之记二 等 水 准 点 之 记 武汉大学樱园 线 点名：B027.5米10.8米40.6米招生就业处政管院中国传统文化研究中心行政楼详细位置图所在图幅H-50-49标石类型混凝土基本水准标石经纬度L：1142150 B：303232标石质料混凝土、铜标志所在地湖北省武汉大学文理学部土地使用者武汉大

18、学地别土质水泥路地下水深度10米交通路线武汉大学文理学部学府路行800米至奥场旁边点 位详细说明1.东40.6米至政治与共管理学院。2.南7.5米至行政楼。3.西南10.8米至招生就业处。接管单位武汉大学测绘学院保管人丁士俊（老师）选点单位11级测绘X班第X组埋石单位11级测绘X班第X组维修单位11级测绘X班第X组选点者XXX埋石者XXX维修者XXX选点日期2013年6月25日埋石日期2013年6月25日维修日期2013年7月3日备注未名湖至行政楼，高差起伏较大。行政楼至新闻与传播学院，地形较为平坦。附录4 自测二等水准测量外业观测高差与高程概算成果表点名测段距离R/km往测高差/mm返测高差

19、/mm往返测不符值/mm不符值累积/mm高差中数闭合改正mm概略高程/mm备注01119 8880.45+19 335.5-19 335.3+0.2+0.2+19 335.4 -0.202139 2230.38- 1 216.9+1 218.7+1.8+2.0-1 217.8 -0.203138 0050.28-10 047.0+10 047.1+0.1+2.1-10 047.0 -0.104127 9580.71-13 821.0+13 823.0+2.0+4.1-13 822.0 -0.305114 1360.22-4 772.4+4 771.2-1.2+2.9-4 771.8-0.106

20、109 3641.06+4 350.0-4 348.0+2.0+4.9+4 349.0 -0.507113 7120.17+ 707.8- 706.7+1.1+6.0+ 707.2 -0.108114 4200.47+5 468.4-5 469.0-0.6+5.4+5 468.7 -0.201119 888R=3.47km n=8 M=R4n=1.0mm h=0.17mm注：01：未名湖 02：行政楼 03：新闻与传播学院 04：城市设计学院 05：田园小观园 06：大自强超市 07：生科院 08：图书馆附录5 大地测量计算源程序1.自定义函数/四舍五入保留n位小数double round(d

21、ouble a,int n=0)for (int i=0,times=1;i=0.5)return (b+1)/times;elsereturn b/times;/度分秒表示法转换为度表示法double dmsTod(double a)double degree = 0, minute = 0, b = 0;degree = floor(a);a -= degree;a *= 100;minute = floor(a);a -= minute;a *= 100;b = degree + minute / 60 + a / 3600;return b;/度表示法转换为度分秒表示法double d

22、Todms(double a)double degree = 0, minute = 0, b = 0;degree = floor(a);a -= degree;a *= 60;minute = floor(a);a -= minute;a *= 60;b = degree + minute / 100 + a / 10000;return b;/椭球参数的获取及相关量计算int getpara(int sign,double &a,double &b,double &e_2,double &ee_2,double &a0,double &a2,double &a4,double &a6,d

23、ouble &a8)a=0;e_2=0;ee_2=0; /椭球参数长半轴、第一偏心率平方、第二偏心率平方switch(sign)Case 0:a = ;b=.;e_2 = 0.2966;ee_2=0.4683;break;case 1:a = ;b=.;e_2 = 0.9588;ee_2=0.9473;break;case 2:a = ;b=.3142 ;e_2 = 0.013;ee_2 = 0.227;break;case 3:a = ;b=.3141 ;e_2 = 0.290;ee_2 = 0.548;break;default:return 1;double m0=0,m2=0,m4=0

24、,m6=0,m8=0;m0 = a * (1 - e_2);m2 = 3 * e_2 * m0 / 2;m4 = 5 * e_2 * m2 / 4;m6 = 7 * e_2 * m4 / 6;m8 = 9 * e_2 * m6 / 8;a0 = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128;a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;a6 = m6 / 32 + m8 / 16; a8 = m8

25、/ 128;return 0;2.高斯投影坐标正反算/高斯投影坐标正算 sign表示椭球，m63表示分带方式(0代表6度带，其他代表3度带)int GaussPro1(int ellipsoid,int m63,double B,double L,double &X,double &Y) const double pi=4*atan(1);double a=0,b=0,e_2=0,ee_2=0; /椭球参数X=0;Y=0; /高斯坐标参数double a0=0,a2=0,a4=0,a6=0,a8=0;getpara(ellipsoid,a,b,e_2,ee_2,a0,a2,a4,a6,a8);

26、B = dmsTod(B); B = B * pi / 180; /B转换为弧度表示 L = dmsTod(L); /L转换为度表示/计算中央子午线经度L0和带号Noint L0=0,No=0; if (!m63)No = static_cast(ceil(L/6);L0 = 6 * No - 3;elseif (L 3)L = L + 360;No = static_cast(round(L/3);L0 = 3 * No; double l=0;l = L - L0; l = l * pi / 180; /l转换为弧度表示X = a0 * B - a2 * sin(2*B) / 2 + a4

27、 * sin(4*B) / 4 -a6 * sin(6*B) / 6 + a8 * sin(8*B) / 8;double sinB=0,cosB=0,sinB_2=0,cosB_2=0,t=0,eta_2=0,N=0,m=0;sinB=sin(B);cosB=cos(B);cosB_2=cosB*cosB;sinB_2=sinB*sinB;t = tan(B);eta_2 = ee_2 * cosB_2;N = a * sqrt(1.0/(1 - e_2 * sinB_2);m = cosB * l;X=X+N*t*(1.0/2+(5-t*t+9*eta_2+4*eta_2*eta_2)/2

28、4+(61-58*t*t+t*t*t*t)*m*m/720)*m*m)*m*m);Y=N*(1+(1-t*t+eta_2)/6+(5-18*t*t+t*t*t*t+14*eta_2-58*eta_2*t*t)*m*m/120)*m*m)*m);Y=No * + + Y;return 0;/高斯投影坐标反算int GaussPro2(int ellipsoid,int m63,double X,double Y,double &B,double &L) const double pi=4*atan(1);double a=0,b=0,e_2=0,ee_2=0; /椭球参数B=0;L=0; /大地

29、坐标参数double a0=0,a2=0,a4=0,a6=0,a8=0;getpara(ellipsoid,a,b,e_2,ee_2,a0,a2,a4,a6,a8);int No=0,L0=0;double x=0,y=0;No = static_cast(floor(Y / );y = Y - No * - ;x = X;if (!m63)L0 = 6 * No - 3;else L0 = (3 * No)%360;/通过迭代求底点纬度Bfdouble eps=1.0E-15,Bf=0,Bf1=0;Bf = X / a0;do Bf1 = Bf;Bf= (X+a2 * sin(2 * Bf1

30、) / 2 - a4 * sin(4 * Bf1) / 4 +a6 * sin(6 * Bf1) / 6 - a8 * sin(8 * Bf1)/ 8)/a0; while (fabs(Bf-Bf1)eps);double sinBf=0,cosBf=0,sinBf_2=0,cosBf_2=0;double tf=0,Nf=0,Mf=0,etaf_2=0,Vf=0,yN=0;sinBf = sin(Bf);cosBf = cos(Bf);cosBf_2 = cosBf * cosBf;sinBf_2 = sinBf * sinBf;tf = tan(Bf);Nf = a / sqrt(1 -

31、e_2 * sinBf_2);Mf=a*(1-e_2)/pow(1-e_2*sinBf_2),1.5);etaf_2 = ee_2 * cosBf_2;Vf=sqrt(etaf_2);yN = y / Nf; double tf_2=0,Nf_3=0,l=0;tf_2=tf*tf;Nf_3 = Nf * Nf * Nf;B = Bf - tf * y * y / 2 / Mf / Nf + tf * (5 + 3 * tf_2 + etaf_2 - 9 * etaf_2 * tf_2) * y * y * y * y/24/Mf/Nf_3 - tf * (61 + 90 * tf_2 + 45

32、 * tf_2 * tf_2) * y * y * y / Nf_3* y * y * y / 720 / Mf /Nf/Nf;l = y / Nf / cosBf - (1 + 2 * tf_2 + etaf_2) * y * y * y / 6 / Nf_3 / cosBf + (5 + 28 * tf_2 + 24 * tf_2 * tf_2 + 6 * etaf_2 + 8 * etaf_2 * tf_2) * y * y * y * y * y / 120 / Nf_3/Nf/Nf / cosBf;B = B * 180 / pi; l = l * 180 / pi;L = L0 +

33、 l;B = dTodms(B);L = dTodms(L);return 0;void CGaussProDlg:OnButton1() / TODO: Add your control notification handler code hereUpdateData(true);if (m_63=-1)MessageBox(请选择投影带！);return;if (m_ellipsoid=-1)MessageBox(请选择椭球！);return;/m_B=37.;m_L=112.;GaussPro1(m_ellipsoid,m_63,m_B,m_L,m_X,m_Y); VarR8Round(

34、m_X,4,&m_X); VarR8Round(m_Y,4,&m_Y);UpdateData(false);void CGaussProDlg:OnButton2() / TODO: Add your control notification handler code hereUpdateData(true);if (m_63=-1)MessageBox(请选择投影带！);return;if (m_ellipsoid=-1)MessageBox(请选择椭球！);return;GaussPro2(m_ellipsoid,m_63,m_X,m_Y,m_B,m_L);VarR8Round(m_B,1

35、0,&m_B);VarR8Round(m_L,10,&m_L);UpdateData(false);3.实测斜距化算至高斯平面边长void CDistanceDlg:OnButton1() / TODO: Add your control notification handler code hereUpdateData(true);const double pi=4*atan(1);double x1=m_x1,y1=m_y1,B1=m_B1,H1=m_H1+50,x2=m_x2,y2=m_y2,B2=m_B2,H2=m_H2+50,D=m_D;double a=0,b=0,e_2=0,ee_

36、2=0,a0=0,a2=0,a4=0,a6=0,a8=0,A=0,N=0,RA=0;getpara(0,a,b,e_2,ee_2,a0,a2,a4,a6,a8);/根据高斯坐标反算大地坐标double B11=0,B22=0,L11=0,L22=0;GaussPro2(0,0,x1,y1+,B11,L11);GaussPro2(0,0,x2,y2+,B22,L22);/利用高斯平均引数法反算点1到点2的大地方位角AGeoSolution2(0,0,B11,L11,B22,L22,A,a0,a2);A=dmsTod(A)*pi/180;B1=dmsTod(B1)*pi/180;B2=dmsTod

37、(B2)*pi/180;double Bm=(B1+B2)/2,Nm=0,Vm=0,Rm=0;Vm=sqrt(1+ee_2*cos(Bm)*cos(Bm); Nm=a/sqrt(1-e_2*sin(Bm)*sin(Bm);Rm=Nm/Vm;N=a/sqrt(1-e_2*sin(B1)*sin(B1);RA=N/(1+ee_2*cos(B1)*cos(B1)*cos(A)*cos(A);double ym=(y1+y2)/2-,S=0,d=0;S=D*sqrt(1-(H2-H1)*(H2-H1)/D/D)/(1+H1/RA)/(1+H2/RA)+D*D*D/24/RA/RA;d=S*(1+ym*

38、ym/2/Rm/Rm+(y2-y1)*(y2-y1)/24/Rm/Rm+ym*ym*ym*ym/24/Rm/Rm/Rm/Rm);VarR8Round(d,4,&d);m_d=d;UpdateData(false);4.大地测量主题正反算/大地主题正算 ellopsoid代表椭球体,solution是解算方法(0是高斯平均引数,其他是白塞尔)int GeoSolution1(int ellipsoid,int solution,double B1,double L1,double A12,double S,double &B2,double &L2,double &A21) const doub

39、le pi=4*atan(1);double a=0,b=0,e_2=0,ee_2=0,a0=0,a2=0,a4=0,a6=0,a8=0; getpara(ellipsoid,a,b,e_2,ee_2,a0,a2,a4,a6,a8);B1 = dmsTod(B1)*pi/180;L1 = dmsTod(L1)*pi/180;A12 = dmsTod(A12)*pi/180;if (!solution)double Bm=B1,Am=A12,Vm=0,Nm=0,tm=0,etam_2=0,deltaB=0,deltaL=0,deltaA=0,deltaB0=0,deltaL0=0,deltaA0

40、=0,sinAm=0,cosAm=0;B2=B1;A21=A12;do deltaB0=deltaB;deltaL0=deltaL;deltaA0=deltaA;Bm=(B1+B2)/2;Am=(A12+A21)/2;sinAm=sin(Am);cosAm=cos(Am);Vm=sqrt(1+ee_2*cos(Bm)*cos(Bm);Nm=a/sqrt(1-e_2*sin(Bm)*sin(Bm);tm=tan(Bm);etam_2=ee_2*cos(Bm)*cos(Bm);deltaB=Vm*Vm/Nm*S*cosAm*(1+S*S/24/Nm/Nm*(sinAm*sinAm*(2+3*tm*

41、tm*(1+etam_2)+3*etam_2*cosAm*cosAm*(-1+tm*tm-etam_2-4*tm*tm*etam_2); deltaL=S/Nm/cos(Bm)*sinAm*(1+S*S/24/Nm/Nm*(sinAm*sinAm*tm*tm-cosAm*cosAm*(1+etam_2-9*tm*tm*etam_2+etam_2*etam_2);deltaA=S/Nm*sinAm*tm*(1+S*S/24/Nm/Nm*(cosAm*cosAm*(2+7*etam_2+9*tm*tm*etam_2+5*etam_2*etam_2)+sinAm*sinAm*(2+tm*tm+2*e

42、tam_2);B2=B1+deltaB;L2=L1+deltaL;A21=A12+deltaA; while (fabs(deltaB-deltaB0)1.0E-15|fabs(deltaL-deltaL0)1.0E-15|fabs(deltaA-deltaA0)1.0E-15);if(A12pi)A21-=pi;elseA21+=pi;elsedouble sinB1=sin(B1),cosB1=cos(B1),sinA12=sin(A12),cosA12=cos(A12),W1=0,sinu1=0,cosu1=0;W1 = sqrt(1 - (e_2 * sinB1 * sinB1);si

43、nu1 = sinB1 * sqrt(1 -e_2) / W1;cosu1 = cosB1 / W1;double sinA0=0,cosA0_2=0,k_2=0,A=0,B=0,C=0,alpha=0,beta=0;sinA0 = cosu1 * sinA12;cosA0_2 = 1 - sinA0 * sinA0;k_2 = ee_2 * cosA0_2;A = (1 + k_2 / 4 - 3 * k_2 * k_2 / 64 + 5 * k_2 * k_2 * k_2 / 256) * b;B = (k_2 / 8 - k_2 * k_2 / 32 + 15 * k_2 * k_2 *

44、 k_2 / 1024) * b;C = (k_2 * k_2 / 128 - 3 * k_2 * k_2 * k_2 / 512) * b;alpha = (e_2 / 2 + e_2 * e_2 / 8 + e_2 * e_2 * e_2 / 16) - (e_2 * e_2 + e_2 * e_2 * e_2) / 16 * cosA0_2 + 3 * e_2 * e_2 * e_2 / 128 * cosA0_2 * cosA0_2;beta = (e_2 * e_2 + e_2 * e_2 * e_2) / 32 * cosA0_2 - e_2 * e_2 * e_2 * cosA0

45、_2 * cosA0_2 / 64;double cotsigma1=0,sin2sigma1=0,cos2sigma1=0;cotsigma1 = cosu1 * cosA12 / sinu1;sin2sigma1 = 2 * cotsigma1 / (cotsigma1 * cotsigma1 + 1);cos2sigma1 = (cotsigma1 * cotsigma1 - 1) / (cotsigma1 * cotsigma1 + 1);double sigma0=0,sin2sigma0=0,cos2sigma0=0;sigma0 = (S - (B + C * cos2sigma

46、1) * sin2sigma1) / A;sin2sigma0=sin(2*sigma0);cos2sigma0=cos(2*sigma0);double sin2=0,cos2=0,sigma=0,cossigma=0,sinsigma=0,delta=0;sin2 = sin2sigma1 * cos2sigma0 + cos2sigma1 * sin2sigma0;cos2 = cos2sigma1 * cos2sigma0 - sin2sigma1 * sin2sigma0;sigma = sigma0 + (B + 5 * C * cos2) * sin2 / A;cossigma = cos(sigma);sinsigma = sin(sigma);delta = (alpha * sigma + beta * (sin2 - sin2sigma1) * sinA0;double sinu2=0;sinu2 = sinu1 * cossigma + cosu1 * cosA12