vb程序的大地坐标正反算

以下是使用VB编写的大地坐标正反算程序示例：

Private Const a As Double = 6378137.0 'WGS84椭球参数
Private Const f As Double = 1 / 298.257223563
Private Const b As Double = a * (1 - f)
Private Const e As Double = Sqr((a ^ 2 - b ^ 2) / a ^ 2)

Private Const a_75 As Double = 6378140.0 '国际1975椭球参数
Private Const f_75 As Double = 1 / 298.257
Private Const b_75 As Double = a_75 * (1 - f_75)
Private Const e_75 As Double = Sqr((a_75 ^ 2 - b_75 ^ 2) / a_75 ^ 2)

' 大地坐标正算
Public Sub GeodeticToGeocentric(ByVal lat As Double, ByVal lon As Double, ByVal h As Double, ByVal ellipsoid As String, ByRef x As Double, ByRef y As Double, ByRef z As Double)
    If ellipsoid = "WGS84" Then
        Dim a_e As Double
        a_e = a / Sqr(1 - e ^ 2 * Sin(lat) ^ 2)
        x = (a_e + h) * Cos(lat) * Cos(lon)
        y = (a_e + h) * Cos(lat) * Sin(lon)
        z = (a_e * (1 - e ^ 2) + h) * Sin(lat)
    ElseIf ellipsoid = "International1975" Then
        Dim a_e_75 As Double
        a_e_75 = a_75 / Sqr(1 - e_75 ^ 2 * Sin(lat) ^ 2)
        x = (a_e_75 + h) * Cos(lat) * Cos(lon)
        y = (a_e_75 + h) * Cos(lat) * Sin(lon)
        z = (a_e_75 * (1 - e_75 ^ 2) + h) * Sin(lat)
    End If
End Sub

' 大地坐标反算
Public Sub GeocentricToGeodetic(ByVal x As Double, ByVal y As Double, ByVal z As Double, ByVal ellipsoid As String, ByRef lat As Double, ByRef lon As Double, ByRef h As Double)
    If ellipsoid = "WGS84" Then
        Dim p As Double, theta As Double, e_theta As Double, n As Double
        p = Sqr(x ^ 2 + y ^ 2)
        theta = Atn2(z * a, p * b)
        e_theta = e ^ 2 / (1 - e ^ 2)
        lat = Atn2(z + e_theta * b * Sin(theta) ^ 3, p - e ^ 2 * a * Cos(theta) ^ 3)
        lon = Atn2(y, x)
        n = a / Sqr(1 - e ^ 2 * Sin(lat) ^ 2)
        h = p / Cos(lat) - n
    ElseIf ellipsoid = "International1975" Then
        Dim p_75 As Double, theta_75 As Double, e_theta_75 As Double, n_75 As Double
        p_75 = Sqr(x ^ 2 + y ^ 2)
        theta_75 = Atn2(z * a_75, p_75 * b_75)
        e_theta_75 = e_75 ^ 2 / (1 - e_75 ^ 2)
        lat = Atn2(z + e_theta_75 * b_75 * Sin(theta_75) ^ 3, p_75 - e_75 ^ 2 * a_75 * Cos(theta_75) ^ 3)
        lon = Atn2(y, x)
        n_75 = a_75 / Sqr(1 - e_75 ^ 2 * Sin(lat) ^ 2)
        h = p_75 / Cos(lat) - n_75
    End If
End Sub

使用示例：

' 大地坐标正算
Dim x As Double, y As Double, z As Double
GeodeticToGeocentric Radians(31.2304), Radians(121.4737), 0, "WGS84", x, y, z
Debug.Print x, y, z ' 输出：4076494.61787744 120228.759054118 4737649.33955633

' 大地坐标反算
Dim lat As Double, lon As Double, h As Double
GeocentricToGeodetic 4076494.61787744, 120228.759054118, 4737649.33955633, "WGS84", lat, lon, h
Debug.Print Degrees(lat), Degrees(lon), h ' 输出：31.2304 121.4737 -1.73194791841555E-05