坐标转换参数求取及坐标转换程序设计.doc

毕业设计
设计题目 坐标转换参数求取及坐标转换程序设计

学生姓名 张威
指导教师 杜继亮
专 业 测绘工程
班 级 测绘12-2班
填写日期 2023/4/6

矿业工程学院

摘 要
坐标系统是测量工作中定位的基础，坐标系统有多种形式和基准，由于各测量工作目的不同，所选用的坐标基准也会不同，根据不同的工作规定需要将不同的坐标系下的坐标进行互相转换。在这些坐标转换的过程中会用到很多坐标转换模型，但是坐标系转换模型过于复杂手算非常困难。本设计为了方便施工时碰到的坐标转换问题，设计运用Visual Basic 6.0编程语言编写程序，用来实现坐标系统之间的转换以及转换参数的求解，例如:大地坐标与空间直角坐标的互相转换、高斯投影正反算、二维坐标转换与四参数计算、三维坐标转换与七参数转换、同参考基准下的坐标换带计算，以及坐标数据的批量解决。
关键字：坐标系统，转换模型，坐标转换，程序设计

Abstract
The base of coordinate system in surveying work. there are many forms and benchmarks in the coordinate system. However, in general engineering, the control point and coordinate. System are the same. So It is necessary to transform the control point. coordinate during the construction process. Due to different purposes of each measurement and the selected. different coordinate references, there will be many different coordinate systems. Coordinate systems used in the measurement work are as follows: WGS-84 World Geodetic System, China Geodetic Coordinate System 2023, National Geodetic Coordinate System 1980, Beijing coordinate system 1954 and Local Coordinate System. There are space rectangular coordinate, geodetic coordinate and plane coordinate in the way of the reference in the same coordinate. According to the requirements of different tasks, we need to convert coordinates under the different coordinate systems. On condition that the coordinates of the reference standard can be obtained. the normal construction work can be done. A lot of coordinate transformation models are used in the process of the coordinate transformation. But the coordinate transformation model is very complex and difficult. Nowadays the conversion formula is suitable for the computerization whose language is easy to learn. So in the design I make use of Visual Basic 6 programming language to realize the transformation between the coordinate system and transformation parameters.
Key words : coordinate systems transformation model coordinate transform programming

目录
摘 要 I
Abstract II
第1章 绪论 1
1.1研究背景和意义 1
1.2国内外研究现状 1
1.3研究的重要内容 2
1.4程序设计思绪方法 2
第2章 坐标理论和知识 3
2.1地球椭球 3
2.2基准 3
2.3测量常用坐标系 4
2.3.1大地坐标系 4
2.3.2空间直角坐标系 4
2.3.3平面坐标系 5
2.3.4地方独立坐标系 5
2.4我国常用坐标系 6
2.4.1 1954年北京坐标系 6
2.4.2 1980国家大地坐标系 6
2.4.3 WGS-84世界大地坐标系 6
2.4.4 2023国家大地坐标系 7
第3章 坐标转换程序设计 8
3.1大地坐标系与空间直角坐标系的转换 9
3.1.1大地坐标转换成空间直角坐标 9
3.1.2空间直角坐标转成换大地坐标 10
3.2高斯平面坐标与大地坐标的转换 13
3.2.1 高斯正算 13
3.2.2 高斯反算 14
3.3 高斯投影换带计算 17
3.4 二维坐标转换与四参数计算 20
3.5 三维坐标转换与七参数计算 23
3.6 三、六度带带号与中央子午线计算 27
第4章 设计总结 31
致谢 32
参考文献 33

第1章 绪论
1.1研究背景和意义
随着大地测量学，摄影测量学的发展和电子计算机的普及，对各种坐标系的研究变得越来越重要了[1]。随着现在社会的快速发展，各种各样的大型工程的建设，凡是工程施工就必然需要坐标来定位才可以，而建造的地方又不同，又没有一个满足所有地方的坐标系统，所以产生多种不同的坐标系统，实际工作中测量人员必然会按照实际情况选择最为合适的坐标系统。而坐标系统之间的转换比较复杂，手算工作量巨大，因而各种坐标转换模型相继出现，运用计算机强大的数据计算能力可以轻松应对这些问题，提高工作效率。坐标转换的意义重大，不仅在我们熟知的工程领域中，在国防建设、航空航天科技、城市汇划等众多领域中都发挥着重要的作用，可以说对社会进步有着必不可少的作用。
1.2国内外研究现状
自60年代以来，各国大地测量学者，通过大量研究，提出了多种坐标转换模型及多种解算方法，北美1927基准面（基于克拉克1966椭球体与北美1983基准面（基于GRS1980椭球体）之间坐标转换是根据研究区内一系列己知点的大地坐标或网格坐标改正量进行插值进行的坐标系转换；英国采用北向与东向的双线性网格插值进行坐标转换；挪威在海岸带调查中，采用经纬度多项式用于坐标系转换这种方法进行新（ED87一欧洲1987基准面）、旧（ED50一欧洲1950基准面）坐标系之间的转换：欧洲石油勘探组织（EPSG）对新、旧坐标系采用“双线性插值” 进行坐标转换[2]。
国内空间三维直角坐标转换中，一般采用布尔莎七参数模型。一般有7个转换参数，即3个平移参数，3个旋转参数和1个尺度参数。需要三个及已经公共点时，才干运用平差的方法求出七参数。
1.3研究的重要内容
本坐标转换程序可实现功能有:1、大地坐标与空间直角坐标的互相转换， 2、高斯投影正反算，3、二维坐标转换与四参数计算，4、三维坐标转换与七参数转换，5、同参考基准下的坐标换带计算，以及坐标数据的批量解决。
1.4程序设计思绪方法
本程序名为万能坐标转换器。设计前期收集相关资料，参考一些成熟的坐标转换软件，拟定程序应有的功能以及界面设计。运用VB编写程序时，查阅相关书籍获取理论知识以及转换模型。完毕程序后将已知对的数据带入其中验证程序结果是否对的。若出现错误则检查每步代码，直到程序完美运营为止。

第2章 基础知识准备
2.1地球椭球
地球椭球体又称地球椭圆体或地球扁球体，代表地球大小和形状的数学曲面，以长半径和扁率表达，因它十分逼近于椭球体，故通常以参考椭球体表达地球椭球体的形状和大小。通常所说地球的形状和大小，事实上就是以参考椭球体的半长径、半短径和扁率来表达。1975年国际大地测量与地球物理联合会推荐的数据为：半长径6378140米，半短径6356755米，扁率1∶298.257。在众多椭球体中，WGS-84椭球体被认为符合上述条件最佳的椭球。
2.2基准
所谓基准是指为描述空间位置而定义的点、线和面。而大地测量基准是指用以描述地球形状的地球椭球参数，包含描述地球椭球几何特性的长短半轴和物理特性的有关参数、地球在空间的定位及定向以及描述这些位置所采用的单位长度的定义。不同的坐标系统会使用的基准也不同，根据参考椭球所选原点位置不同，可以分为地心坐标系和参心坐标系。
地心坐标系是以地球的质心为原点，有地心大地坐标系和地心空间直角坐标系两种表述方法。地心空间直角坐标系的定义为：以地球质心为原点，X轴指向格林尼治子午面与地球赤道的交点，Z轴指向北极，Y轴过原点垂直于平面 XOZ，构成右手空间直角坐标系。地心大地坐标系定义为：以地球的质心作为原点，以地球自转轴作为椭球的短轴，大地纬度B是过地面点的椭球法线与椭球赤道面之间的夹角，大地经度L为过地面点的椭球子午面与格林尼治子午面之间的夹角，大地高度H为地面点沿椭球法线到椭球面的最短距离。WGS-84坐标系，CGCS2023坐标系，GLONASS是采用PZ-90坐标，都是属于地心坐标系。
参心坐标系是选取一个参考椭球面作为基本的参考面，选一参考点作为大地测量的起算点，从而拟定参考椭球在地球面的位置和方向。这时参考椭球的原点不会和地球质心重合，所以称为参心北京54坐标系、西安80坐标系和新北京54坐标系，都是参心坐标系。它同样具有参心大地坐标系和参心空间直角坐标系两种表述方法，它们的定义与地心坐标系的定义相似。
2.3测量常用坐标系
2.3.1大地坐标系
空间一点的大地坐标用大地经度L、大地纬度B和大地高度H表达，地面上P地点的大地子午面NPS与起始大地子午面所构成的二面角L称P地点的大地经度， P地点对于椭球的法线与赤道面的夹角B称P地点的大地纬度。如图2-1所示

图2-1大地坐标系
P地点沿法线到椭球面的距离H称大地高，从椭球面起算，向外为正，向内为负[3]。
H = H正常 + ζ（高程异常）
H = H正 + N（大地水准面差距）
2.3.2空间直角坐标系
空间直角坐标系的坐标原点与参考椭球的中心重合，Z轴正向指向参考椭球的北极，X轴正向指向起始子午面与赤道的交点，Y轴按右手系与X轴呈90°夹角且位于赤道面上，某点在空间中的坐标可用该点在此空间坐标系的各个坐标轴上的投影来表达[4]。如图2-2所示：

图2-2空间直角坐标系
2.3.3平面坐标系
平面直角坐标系是运用投影，将空间坐标通过某种数学变换映射到平面上，这种变换称为投影变换[5]。在我国一般采用的是高斯一克吕格投影，是目前测量上广泛采用的正形投影，特点是没有角度变形，在不同点上的长度比随点位而异，但在同一点上各方向的长度比相同，也称为高斯投影[6]。
2.3.4地方独立坐标系
在我国平面坐标重要采用的是高斯投影，在该投影中，除中央子午线外，其它位置上的任何线段，投影后都会产生一定的长度变形，并且变形随离开中央子午线的距离增长而增长[7]。因此一般采用分带投影的办法，来限制长度变形，我国规定了采用3度带或6度带进行分带投影。在城市、工矿等工程测量中，假如直接在国家分带坐标系中建立控制网，会使地面长度投影的变形较大，当长度变形大于2．5 cm/km时，就难以满足工程上的需要[8]。另一些特殊的测量，比如大桥施工测量，水利水坝测量，滑坡变形监测等，采用国家坐标系精度不能满足工程规定，所以经常会建立适合本地区的地方独立坐标系[9]。
2.4我国常用坐标系
2.4.1 1954年北京坐标系
1954年北京坐标系，是采用先将我国一等锁与前苏联远东一等锁相联接，然后以连接处呼玛，吉拉林，东宁基线网扩大边端点的前苏联1942年普尔科沃坐标系的坐标为起算数据，平差我国东北及东部一等锁[10]。椭球参数：长半轴a=6378245m，短半轴b=6356863.0188m，扁率α=1/298.3，第一偏心率平方e2=0.54287，第二偏心率平方e’2=0.14652。
2.4.2 1980国家大地坐标系
1980国家大地坐标系采用地球椭球基本参数为1975年国际大地测量与地球物理联合会第十六届大会推荐的数据，大地原点设在我国中部的陕西省泾阳县永乐镇，位于西安市西北方向约60公里。椭球参数：长半轴a=6378140±5m，短半轴b=6356755.2882m，扁率α=1/298.257，第一偏心率平方e2 =0.59，第二偏心率平方e’2=0.47。
2.4.3 WGS-84世界大地坐标系
WGS-84坐标系是一种国际上采用的地心坐标系。坐标原点为地球质心，其地心空间直角坐标系的Z轴指向BIH （国际时间服务机构）1984.0定义的协议地球极（CTP）方向，X轴指向BIH 1984.0的零子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系称为1984年世界大地坐标系统。椭球参数：长半轴a=6378137，短半轴b=6356752.3142，扁率α=1/298.2572236，第一偏心率平方e2=0.13 ，第二偏心率平方e’2=0.227。
2.4.4 2023国家大地坐标系
CGCS2023是（中国）2023国家大地坐标系的缩写，该坐标系是通过中国GPS连续运营基准站、空间大地控制网以及天文大地网与空间地网联合平差建立的地心大地坐标系统[11]。Z轴指向BIH1984.0定义的协议极地方向（BIH国际时间局），X轴指向BIH1984.0定义的零子午面与协议赤道的交点，Y轴按右手坐标系拟定。椭球参数：长半轴a=6378137，短半轴b=6356752.31414，扁率α=1/298.2572236，第一偏心率平方e2=0.13 ，第二偏心率平方e’2=0.227。
2.5坐标转换模型
2.5.1大地坐标系与空间直角坐标系转换模型
将同一坐标系下的大地坐标（B，L，H）转换成空间直角坐标（X，Y，Z）的转换公式为：

式中：a为参考椭球长半轴，e为第一偏心率，N为卯酉圈的半径[12]。
将同一坐标系下的空间直角坐标（X，Y，Z）转换为大地坐标（B，L，H）的公式为:
L=arctan()
B=arctan()
H=
用公式进行空间直角坐标转换大地坐标时，需要采用迭代计算大地纬度B。具体计算时，可先根据下式求出大地纬度B的初值：

由于 , 带入B=arctan()
得:
令:
则式子可写成
然后运用该初值B代入公式右端tanB中，将等式左边的结果再次代入右端tanB，直到最后两次B值之差小与允许误差为止。当得到大地纬度B后，代入公式即可求出大地高H。
2.5.2高斯正反算转换模型
得到了点的大地坐标(L，B)，就可以将其转化为某投影带的高斯平面坐标，我们将椭球参数代入高斯投影正算公式得到更合用于电算的高斯坐标计算的实用公式：

其计算结果的精度可达0. 001m。
只要得到了高斯平面坐标(X, Y)后，便可通过高斯反算公式将其转换成大地坐标(B, L)，高斯投影反算公式为:

它们的计算精度,即平面坐标可达。
2.5.3坐标转换与参数计算转换模型
二维转换模型：
2个平移参数（原点不重合产生的）
1个旋转参数（坐标轴不平行产生的）
1个尺度参数（两坐标系间的尺度不一致产生的）

为某点在A空间直角坐标系中的坐标。为某点在B空间直角坐标系中的坐标。为某点从A空间直角坐标系转换到B空间直角坐标系的两个平移参数。β为从A空间直角坐标系转换到B空间直角坐标系中标系的两个平移参数。β为从A空间直角坐标系转换到B空间直角坐标系中一个旋转参数。m为从A空间直角坐标系转换到B空间直角坐标系中的一个尺度参数。
平面四参数求解环节如下:
运用公共点计算坐标参数，但至少有两个公共点，当有i个公共点时，可运用最小二乘原理求解参数。
将B直角坐标系中的坐标视为观测值，设A直角坐标系下的坐标视为无误差，列误差方程为：

写成矩阵形式既：
由于各点的坐标可视为同精度独立观测值，因此P=I。
把各点坐标已知值带入上述误差方程，然后按下列公式求解出四参数：

三维转换模型：
3个平移参数（原点不重合产生的）
3个旋转参数（坐标轴不平行产生的）
1个尺度参数（两坐标系间的尺度不一致产生的）

一般为微小转角，模型可简化为：

为某点在A空间直角坐标系中的坐标。为某点在B空间直角坐标系中的坐标。为某点从A空间直角坐标系转换到B空间直角坐标系的三个平移参数。为从A空间直角坐标系转换到B空间直角坐标系中三个旋转参数。m为从A空间直角坐标系转换到B空间直角坐标系中的一个尺度参数。
七参数求解环节如下:
运用公共点计算坐标参数，但至少有三个公共点，当有i个公共点时，可运用最小二乘原理求解参数。
将B空间直角坐标系中的坐标视为观测值，设A空间直角坐标系下的坐标视为无误差，列误差方程为：

写成矩阵形式既：
由于各点的坐标可视为同精度独立观测值，因此P=I。
把各点坐标已知值带入上述误差方程，然后按下列公式求解出七参数：

2.5.4三、六度带带号与中央子午线计算模型
求带号、中央子午线公式为：
3°带:带号 中央子午线 L为本地经度
6°带:带号 中央子午线 L为本地经度

第3章 坐标转换程序设计
计划程序有哪些功能

程序界面设计

搜集坐标转换公式

运用VB编写程序

检查能否正常运营
以及数据对的性
否
是
完毕程序

图3-1 编写程序流程图

图3-2程序主界面
3.1大地坐标系与空间直角坐标系的转换
重要代码如下：
椭球参数设立
If Combo1.Text = "WGS-84" Then
a1 = 6378137: b1 = 6356752.3142 'a1为长半轴，b1为短半轴
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "CGCS2023" Then
a1 = 6378137: b1 = 6356752.31414
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "北京54" Then
a1 = 6378245: b1 = 6356863.0188
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "西安80" Then
a1 = 6378140: b1 = 6356755.2882
Text7.Text = a1: Text8.Text = b1
ElseIf Combo1.Text = "自定义" Then
Text7.Enabled = True: Text8.Enabled = True
a1 = Text7.Text: b1 = Text8.Text
End If
大地坐标转空间直角坐标重要代码
B = Text1.Text: L = Text2.Text: H = Text3.Text
e12 = (a1 ^ 2 - b1 ^ 2) / a1 ^ 2
N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B)) ^ 2))
x = (N + H) * Cos(DuToHu(B)) * Cos(DuToHu(L))
y = (N + H) * Cos(DuToHu(B)) * Sin(DuToHu(L))
z = (N * (1 - e12) + H) * Sin(DuToHu(B))
Text4.Text = x: Text5.Text = y: Text6.Text = z
空间直角坐标转大地坐标重要代码
x = Text4.Text: y = Text5.Text: z = Text6.Text
L = HuToDu(Atn(y / x)) '计算经度
If y / x < 0 Then
L = HuToDu(Atn(y / x) + 3.979)
End If
e12 = (a1 ^ 2 - b1 ^ 2) / a1 ^ 2 '第一偏心率平方
TB = z / Sqr(x ^ 2 + y ^ 2) '初始纬度B正切值
c1 = (a1 ^ 2) / b1
t = z / Sqr(x ^ 2 + y ^ 2)
p = (c1 * e12) / Sqr(x ^ 2 + y ^ 2)
k = 1 + (a1 ^ 2 - b1 ^ 2) / b1 ^ 2
TB1 = t + (p * TB) / Sqr(k + TB ^ 2) '用TanB作为未知量
Do Until Abs(TB1 - TB) < 0. '迭代法
TB = TB1
TB1 = t + (p * TB) / Sqr(k + TB ^ 2)
Loop
B = Atn(TB)
N = a1 / Sqr(1 - e12 * (Sin(B) ^ 2)) '卯酉圈曲率半径
H = z / Sin(B) - N * (1 - e12) '大地高
Text1.Text = HuToDu(B): Text2.Text = L: Text3.Text = H

3.2高斯平面坐标与大地坐标的转换
高斯平面坐标(X, Y)与大地坐标(B ,L)的互相关系式分为两类：第一类称高斯投影正算公式，由(B,L)求(X,Y)；第二类称高斯投影反算公式，由(X,Y)求(B,L) [13]。由于电子计算机和各种可编程序电子计算器在测量上广泛应用，因此在这里给出适合电算的高斯正反算公式。
重要代码如下：
正算
Dim B, L, x, y As Double
Dim N, a1, b1, c, e12, e22, e1, e2, ll, l0, nn, t As Double
Dim XX, bb0, bb2, bb4, bb6, bb8 As Double
Dim yj As Double
Dim mf, nf, nnf, tf, bf, bf1 As Double
B = Text1.Text: L = Text2.Text
l0 = Text5.Text '输入中央子午线经度
yj = Text8.Text 'Y坐标加常数
a1 = Text6.Text: b1 = Text7.Text '椭球长短半轴
e12 = a1 ^ 2 - b1 ^ 2) / a1 ^ 2 '第一偏心率平方
e22 = (a1 ^ 2 - b1 ^ 2) / b1 ^ 2 '第二偏心率平方
e1 = Sqr((a1 ^ 2 - b1 ^ 2) / a1 ^ 2) '第一偏心率
e2 = Sqr((a1 ^ 2 - b1 ^ 2) / b1 ^ 2) '第二偏心率
c = a1 / Sqr(1 - e12) '极点子午线曲率半径
bb0 = 1 - (3 / 4) * e22 + (45 / 64) * e22 ^ 2 - (175 / 256) * e22 ^ 3 + (11025 / 16384) * e22 ^ 4
bb2 = bb0 - 1
bb4 = (15 / 32) * e22 ^ 2 - (175 / 384) * e22 ^ 3 + (3675 / 8192) * e22 ^ 4
bb6 = -(35 / 96) * e22 ^ 3 + (735 / 2048) * e22 ^ 4
bb8 = (315 / 1024) * e22 ^ 4
ll = (Du(L) - l0) / 57.
N = a1 / Sqr(1 - e12 * (Sin(DuToHu(B)) ^ 2)) '卯酉圈曲率半径
t = Sin(DuToHu(B)) / Cos(DuToHu(B))

nn = e2 * Cos(DuToHu(B))
XX = c * (bb0 * DuToHu(B) + (bb2 * Cos(DuToHu(B)) + bb4 * (Cos(DuToHu(B)) ^ 3) + bb6 * (Cos(DuToHu(B)) ^ 5) + bb8 * (Cos(DuToHu(B)) ^ 7)) * Sin(DuToHu(B))) '子午线弧长
x = XX + (N / 2) * Sin(DuToHu(B)) * Cos(DuToHu(B)) * (ll ^ 2) + (N / 24) * Sin(DuToHu(B)) * (Cos(DuToHuB)) ^ 3) * (ll ^ 4) * (5 - t ^ 2 + 9 * (nn ^ 2) + 4 * (nn ^ 4)) + (N / 720) * Sin(DuToHu(B)) *(Cos(DuToHu(B)) ^ 5) * (ll ^ 6) * (61 - 58 * (t ^ 2) + (t ^ 4))
y = N * Cos(DuToHu(B)) * ll + (N / 6) * (Cos(DuToHu(B)) ^ 3) * (ll ^ 3) * (1 - t ^ 2 + nn ^ 2) + (N / 120) * (Cos(DuToHu(B)) ^ 5) * (ll ^ 5) * (5 - 18 * (t ^ 2) + t ^ 4 + 14 * (nn ^ 2) - 58 * (t ^ 2) * (nn ^ 2))
反算
x = Text3.Text: y = Text4.Text
y = y - yj
bf = x / (c * bb0) '迭代
Do
bf1 = ((x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) ^ 3) + bb6 * (Cos(bf) ^ 5) + bb8 * (Cos(bf) ^ 7)) * Sin(bf)) / bb0
Loop Until Abs(bf - bf1) < 0.

mf = (a1 * (1 - e12)) / ((Sqr(1 - e12 * (Sin(bf) ^ 2))) ^ 3)
tf = Sin(bf) / Cos(bf)
nnf = e2 * Cos(bf)
nf = a1 / Sqr(1 - e12 * (Sin(bf) ^ 2))
B = bf - (tf * (y ^ 2)) / (2 * mf * nf) + (tf * (y ^ 4) * (5 + 3 * (tf ^ 2) + (nnf ^ 2) - 9 * (nnf ^ 2) * (tf ^ 2))) / (24 * mf * (nf ^ 3)) - (tf * (y ^ 6) * (61 + 90 * (tf ^ 2) + 45 * (tf ^ 4))) / (720 * mf * (nf ^ 5))
L = y / (nf * Cos(bf)) - (1 + 2 * tf ^ 2 + nnf ^ 2) * (y ^ 3) / (6 * nf ^ 3 * Cos(bf)) + (5 + 28 * (tf ^ 2) + 24 * (tf ^ 4) + 6 * (nnf ^ 2) + 8 * (nnf ^ 2) * (tf ^ 2)) * (y ^ 5) / (120 * (nf ^ 5) * Cos(bf))
3.3 高斯投影换带计算
高斯投影虽然保证了角度没有变形，但是长度变形比较严重，为了限制高斯投影的长度变形，必须以中央子午线进行分带，把投影范围限制在中央子午线东、西侧一定的范围内；但是这样又使得统一的坐标系分割成各带的独立坐标系[14]。这样就会产生新的问题，需要邻带换算来解决。
重要代码如下：
l1 = Text3.Text : yy1 = Text4.Text '加常数和中央子午线
x = Text9.Text : y = Text10.Text
反算
y = y - yy1
e12 = (a1 ^ 2 - b1 ^ 2) / a1 ^ 2 '第一偏心率平方
e22 = (a1 ^ 2 - b1 ^ 2) / b1 ^ 2 '第二偏心率平方
e1 = Sqr((a1 ^ 2 - b1 ^ 2) / a1 ^ 2) '第一偏心率
e2 = Sqr((a1 ^ 2 - b1 ^ 2) / b1 ^ 2) '第二偏心率
c = a1 / Sqr(1 - e12) '极点子午线曲率半径
bb0 = 1 - (3 / 4) * e22 + (45 / 64) * e22 ^ 2 - (175 / 256) * e22 ^ 3 + (11025 / 16384) * e22 ^ 4
bb2 = bb0 - 1
bb4 = (15 / 32) * e22 ^ 2 - (175 / 384) * e22 ^ 3 + (3675 / 8192) * e22 ^ 4
bb6 = -(35 / 96) * e22 ^ 3 + (735 / 2048) * e22 ^ 4
bb8 = (315 / 1024) * e22 ^ 4
bf1 = x / (c * bb0) '迭代 单位弧度
Do
bf = bf1
bf1 = ((x / c) - (bb2 * Cos(bf) + bb4 * (Cos(bf) ^ 3) + bb6 * (Cos(bf) ^ 5) + bb8 * (Cos(bf) ^ 7)) * Sin(bf)) / bb0
Loop Until Abs(bf - bf1) < 0.
mf = (a1 * (1 - e12)) / ((Sqr(1 - e12 * (Sin(bf) ^ 2))) ^ 3)
tf = Sin(bf) / Cos(bf)
nnf = e2 * Cos(bf)
nf = a1 / Sqr(1 - e12 * (Sin(bf) ^ 2))
'B L单位弧度
B = bf - (tf * (y ^ 2)) / (2 * mf * nf) + (tf * (y ^ 4) * (5 + 3 * (tf ^ 2) + (nnf ^ 2) - 9 * (nnf ^ 2) * (tf ^ 2))) / (24 * mf * (nf ^ 3)) - (tf * (y ^ 6) * (61 + 90 * (tf ^ 2) + 45 * (tf ^ 4))) / (720 * mf * (nf ^ 5))
L = y / (nf * Cos(bf)) - (1 + 2 * tf ^ 2 + nnf ^ 2) * (y ^ 3) / (6 * nf ^ 3 * Cos(bf)) + (5 + 28 * (tf ^ 2) + 24 * (tf ^ 4) + 6 * (nnf ^ 2) + 8 * (nnf ^ 2) * (tf ^ 2)) * (y ^ 5) / (120 * (nf ^ 5) * Cos(bf))
L = L DuToHu(l1)
正算
l2 = Text7.Text: yy2 = Text8.Text '加常数和中央子午线
ee12 = (a1 ^ 2 - b1 ^ 2) / a1 ^ 2 '第一偏心率平方
ee22 = (a1 ^ 2 - b1 ^ 2) / b1 ^ 2 '第二偏心率平方
ee1 = Sqr((a1 ^ 2 - b1 ^ 2) / a1 ^ 2) '第一偏心率
ee2 = Sqr((a1 ^ 2 - b1 ^ 2) / b1 ^ 2) '第二偏心率
nn = a1 / Sqr(1 - ee12 * (Sin(B) ^ 2)) '卯酉圈曲率半径
tt = Sin(B) / Cos(B)
nnn = ee2 * Cos(B)
cc = a1 / Sqr(1 - ee12) '极点子午线曲率半径
bbb0 = 1 - (3 / 4) * ee22 + (45 / 64) * ee22 ^ 2 - (175 / 256) * ee22 ^ 3 + (11025 / 16384) * ee22 ^ 4
bbb2 = bbb0 - 1
bbb4 = (15 / 32) * ee22 ^ 2 - (175 / 384) * ee22 ^ 3 + (3675 / 8192) * ee22 ^ 4
bbb6 = -(35 / 96) * ee22 ^ 3 + (735 / 2048) * ee22 ^ 4
bbb8 = (315 / 1024) * ee22 ^ 4
L = HuToDu(L)
ll = (Du(L) - l2) / 57.
XX = cc * (bbb0 * B + (bbb2 * Cos(B) + bbb4 * (Cos(B) ^ 3) + bbb6 * (Cos(B) ^ 5) + bbb8 * (Cos(B) ^ 7)) * Sin(B)) '子午线弧长
x = XX + (nn / 2) * Sin(B) * Cos(B) * (ll ^ 2) + (nn / 24) * Sin(B) *