大地主题解算深度干货超精

大地主题解算(正算)代码与白塞尔大地主题解算

大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标

新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：

起点经度：116.235(度)

终点纬度：37.435(度)

方向角：50（度）

长度：500（米）

终点经纬度（"经度,纬度"）

=Computation(37.435,116.235,50,500)

Const Pi = 3.14

Private a, b, c, alpha, e, e2, W, V As Double

'a 长轴半径

'b 短轴

'c 极曲率半径

'alpha 扁率

'e 第一偏心率

'e2 第二偏心率

'W 第一基本纬度函数

'V 第二基本纬度函数

Private B1, L1, B2, L2 As Double

'B1 点1的纬度

'L1 点1的经度

'B2 点1的纬度

'L2 点2的经度

Private S As Double '''''大地线长度

Private A1, A2 As Double

'A1 点1到点2的方位角

'A2 点2到点1的方位角

Function Computation(STARTLAT, STARTLONG, ANGLE1, DISTANCE As Double) As String

B1 = STARTLAT

L1 = STARTLONG

A1 = ANGLE1

S = DISTANCE

a = 6378245

b = 6356752.3142

c = a ^ 2 / b

alpha = (a - b) / a

e = Sqr(a ^ 2 - b ^ 2) / a

e2 = Sqr(a ^ 2 - b ^ 2) / b

B1 = rad(B1)

L1 = rad(L1)

A1 = rad(A1)

W = Sqr(1 - e ^ 2 * (Sin(B1) ^ 2))

V = W * (a / b)

Dim W1 As Double

E1 = e ''''第一偏心率

'// 计算起点的归化纬度

W1 = W ''Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 )) sinu1 = Sin(B1) * Sqr(1 - E1 * E1) / W1

cosu1 = Cos(B1) / W1

'// 计算辅助函数值

sinA0 = cosu1 * Sin(A1)

cotq1 = cosu1 * Cos(A1)

sin2q1 = 2 * cotq1 / (cotq1 ^ 2 + 1)

cos2q1 = (cotq1 ^ 2 - 1) / (cotq1 ^ 2 + 1)

'// 计算系数AA,BB,CC及AAlpha, BBeta的值。

cos2A0 = 1 - sinA0 ^ 2

e2 = Sqr(a ^ 2 - b ^ 2) / b

k2 = e2 * e2 * cos2A0

Dim aa, BB, CC, EE22, AAlpha, BBeta As Double

aa = b * (1 + k2 / 4 - 3 * k2 * k2 / 64 + 5 * k2 * k2 * k2 / 256)

BB = b * (k2 / 8 - k2 * k2 / 32 + 15 * k2 * k2 * k2 / 1024)

CC = b * (k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512)

e2 = E1 * E1

AAlpha = (e2 / 2 + e2 * e2 / 8 + e2 * e2 * e2 / 16) - (e2 * e2 / 16 + e2 * e2 * e2 / 16) * cos2A0 + (3 * e2 * e2 * e2 / 128) * cos2A0 * cos2A0

BBeta = (e2 * e2 / 32 + e2 * e2 * e2 / 32) * cos2A0 - (e2 * e2 * e2 / 64) * cos2A0 * cos2A0

'// 计算球面长度

q0 = (S - (BB + CC * cos2q1) * sin2q1) / aa

sin2q1q0 = sin2q1 * Cos(2 * q0) + cos2q1 * Sin(2 * q0)

cos2q1q0 = cos2q1 * Cos(2 * q0) - sin2q1 * Sin(2 * q0)

q = q0 + (BB + 5 * CC * cos2q1q0) * sin2q1q0 / aa

'// 计算经度差改正数

theta = (AAlpha * q + beta * (sin2q1q0 - sin2q1)) * sinA0

'// 计算终点大地坐标及大地方位角

sinu2 = sinu1 * Cos(q) + cosu1 * Cos(A1) * Sin(q)

B2 = Atn(sinu2 / (Sqr(1 - E1 * E1) * Sqr(1 - sinu2 * sinu2))) * 180 / Pi

lamuda = Atn(Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1))) * 180 / Pi

If (Sin(A1) > 0) Then

If (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Then

lamuda = Abs(lamuda)

Else