全局代码
Const pi = 3.1415926535898
Private a, b, c, alpha, e, e2, w, V As Double
Private B1, L1, B2, L2 As Double
Private s As Double
Private A1, A2 As Double
Private Sub getellipseparameter()
a = 6378245
b = 6356752.3142
c = a ^ 2 / b
alpha = (a - b) / a
e = math.sqrt(a ^ 2 - b ^ 2) / a
e2 = math.sqrt(a ^ 2 - b ^ 2) / b
End Sub
Private Function computerw()
w = math.sqrt(1 - e ^ 2 * (math.sin(B1) ^ 2))
V = w * (a / b)
End Function
Public Function Computation(STARTLAT As Double, STARTLONG As Double, ANGLE1 As Double, DISTANCE As Double) As String '''''正算
Dim sinu1, cosu1, sinA0, cotq1, sin2q1, cos2q1, cos2A0 As Double
Dim k2, q0, sin2q1q0, cos2q1q0 As Double
Dim q As Double
Dim theta As Double
Dim aa, BB, cc, EE22, AAlpha, BBeta As Double
Dim sinu2, lamuda As Double
Dim e1 As Double
Dim W1 As Double
B1 = STARTLAT
L1 = STARTLONG
A1 = ANGLE1
s = DISTANCE
Call getellipseparameter
If B1 = 0 Then
If A1 = 90 Then
A2 = 270
B2 = 0
L2 = L1 + s / a * 180 / pi
End If
If A1 = 270 Then
A2 = 90
B2 = 0
L2 = L1 - s / a * 180 / pi
End If
Exit Function
End If
B1 = rad(B1)
L1 = rad(L1)
A1 = rad(A1)
Call computerw
e1 = e
W1 = w
sinu1 = math.sin(B1) * math.sqrt(1 - e1 * e1) / W1
cosu1 = math.cos(B1) / W1
sinA0 = cosu1 * math.sin(A1)
cotq1 = cosu1 * math.cos(A1)
sin2q1 = 2 * cotq1 / (cotq1 ^ 2 + 1)
cos2q1 = (cotq1 ^ 2 - 1) / (cotq1 ^ 2 + 1)
cos2A0 = 1 - sinA0 ^ 2
e2 = math.sqrt(a ^ 2 - b ^ 2) / b
k2 = e2 * e2 * cos2A0
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 * math.cos(2 * q0) + cos2q1 * math.sin(2 * q0)
cos2q1q0 = cos2q1 * math.cos(2 * q0) - sin2q1 * math.sin(2 * q0)
q = q0 + (BB + 5 * cc * cos2q1q0) * sin2q1q0 / aa
'theta = (AAlpha * q + BBeta * (sin2q1q0 - sin2q1)) * sinA0
theta = (AAlpha * q + BBeta * (sin2q1q0 - sin2q1)) * sinA0
sinu2 = sinu1 * math.cos(q) + cosu1 * math.cos(A1) * math.sin(q)
B2 = math.atan(sinu2 / (math.sqrt(1 - e1 * e1) * math.sqrt(1 - sinu2 * sinu2))) * 180 / pi
lamuda = math.atan(math.sin(A1) * math.sin(q) / (cosu1 * math.cos(q) - sinu1 * math.sin(q) * math.cos(A1))) * 180 / pi
If (math.sin(A1) > 0) Then
If (math.sin(A1) * math.sin(q) / (cosu1 * math.cos(q) - sinu1 * math.sin(q) * math.cos(A1)) > 0) Then
lamuda = math.abs(lamuda)
Else
lamuda = 180 - math.abs(lamuda)
End If
Else
If (math.sin(A1) * math.sin(q) / (cosu1 * math.cos(q) - sinu1 * math.sin(q) * math.cos(A1)) > 0) Then
lamuda = math.abs(lamuda) - 180
Else
lamuda = -math.abs(lamuda)
End If
End If
L2 = L1 * 180 / pi + lamuda - theta * 180 / pi
A2 = math.atan(cosu1 * math.sin(A1) / (cosu1 * math.cos(q) * math.cos(A1) - sinu1 * math.sin(q))) * 180 / pi
If (math.sin(A1) > 0) Then
If (cosu1 * math.sin(A1) / (cosu1 * math.cos(q) * math.cos(A1) - sinu1 * math.sin(q)) > 0) Then
A2 = 180 + math.abs(A2)
Else
A2 = 360 - math.abs(A2)
End If
Else
If (cosu1 * math.sin(A1) / (cosu1 * math.cos(q) * math.cos(A1) - sinu1 * math.sin(q)) > 0) Then
A2 = math.abs(A2)
Else
A2 = 180 - math.abs(A2)
End If
End If
Computation = format(L2, "0.00000000") & "," & format(B2, "0.00000000")
End Function
Private Function rad(ByVal angle_d As Double) As Double
rad = angle_d * pi / 180
End Function