大地主题解算(正算)代码

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

新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：
起点经度：116.235(度)
终点纬度：37.435(度)
方向角：50（度）
长度：500（米）
终点经纬度（"经度,纬度"）=Computation(37.435,116.235,50,500)

Const Pi = 3.1415926535898
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))) * 1

80 / 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
lamuda = 180 - Abs(lamuda)
End If
Else
If (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Then
lamuda = Abs(lamuda) - 180
Else
lamuda = -Abs(lamuda)
End If
End If
L2 = L1 * 180 / Pi + lamuda - theta * 180 / Pi

A2 = Atn(cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q))) * 180 / Pi
If (Sin(A1) > 0) Then
If (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) Then
A2 = 180 + Abs(A2)
Else
A2 = 360 - Abs(A2)
End If
Else
If (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) Then
A2 = Abs(A2)
Else
A2 = 180 - Abs(A2)
End If
End If

Computation = L2 & "," & B2
End Function

Private Function rad(ByVal angle_d As Double) As Double
rad = angle_d * Pi / 180
End Function