Function VN2000_To_WGS84_DMS(X As Double, Y As Double) As String
Dim a As Double: a = 6378137#
Dim f As Double: f = 1 / 298.257223563
Dim k0 As Double: k0 = 0.9999
Dim FE As Double: FE = 500000#
Dim FN As Double: FN = 0#
Dim lon0 As Double: lon0 = 108.5
Dim b As Double: b = a * (1 - f)
Dim e2 As Double: e2 = (a ^ 2 - b ^ 2) / a ^ 2
Dim ep2 As Double: ep2 = e2 / (1 - e2)
Dim M As Double: M = (X - FN) / k0
Dim mu As Double
mu = M / (a * (1 - e2 / 4 - 3 * e2 ^ 2 / 64 - 5 * e2 ^ 3 / 256))
Dim e1 As Double
e1 = (1 - Sqr(1 - e2)) / (1 + Sqr(1 - e2))
Dim J1 As Double, J2 As Double, J3 As Double, J4 As Double
J1 = (3 * e1 / 2 - 27 * e1 ^ 3 / 32)
J2 = (21 * e1 ^ 2 / 16 - 55 * e1 ^ 4 / 32)
J3 = (151 * e1 ^ 3 / 96)
J4 = (1097 * e1 ^ 4 / 512)
Dim fp As Double
fp = mu + J1 * Sin(2 * mu) + J2 * Sin(4 * mu) + J3 * Sin(6 * mu) + J4 * Sin(8 *
mu)
Dim C1 As Double, T1 As Double, N1 As Double, R1 As Double, D As Double
C1 = ep2 * Cos(fp) ^ 2
T1 = Tan(fp) ^ 2
N1 = a / Sqr(1 - e2 * Sin(fp) ^ 2)
R1 = N1 * (1 - e2) / (1 - e2 * Sin(fp) ^ 2)
D = (Y - FE) / (N1 * k0)
Dim lat_rad As Double, lon_rad As Double
lat_rad = fp - (N1 * Tan(fp) / R1) * (D ^ 2 / 2 - (5 + 3 * T1 + 10 * C1 - 4 *
C1 ^ 2 - 9 * ep2) * D ^ 4 / 24)
lon_rad = (D - (1 + 2 * T1 + C1) * D ^ 3 / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 ^
2 + 8 * ep2 + 24 * T1 ^ 2) * D ^ 5 / 120) / Cos(fp)
lon_rad = lon_rad + lon0 * [Link]() / 180
Dim lat_deg As Double, lon_deg As Double
lat_deg = lat_rad * 180 / [Link]()
lon_deg = lon_rad * 180 / [Link]()
Dim N As Double
N = a / Sqr(1 - e2 * Sin(lat_rad) ^ 2)
Dim X0 As Double, Y0 As Double, Z0 As Double
X0 = N * Cos(lat_rad) * Cos(lon_rad)
Y0 = N * Cos(lat_rad) * Sin(lon_rad)
Z0 = (N * (1 - e2)) * Sin(lat_rad)
Dim dx As Double: dx = -191.904414
Dim dy As Double: dy = -39.303182
Dim dz As Double: dz = -111.450328
Dim wx As Double: wx = 0.00928836 * [Link]() / (180 * 3600)
Dim wy As Double: wy = 0.01975479 * [Link]() / (180 * 3600)
Dim wz As Double: wz = 0.00427372 * [Link]() / (180 * 3600)
Dim m As Double: m = 1 + 1.0000631E-06
Dim X1 As Double, Y1 As Double, Z1 As Double
X1 = dx + m * (X0 + wz * Y0 - wy * Z0)
Y1 = dy + m * (-wz * X0 + Y0 + wx * Z0)
Z1 = dz + m * (wy * X0 - wx * Y0 + Z0)
Dim p As Double, theta As Double, e_dash2 As Double
p = Sqr(X1 ^ 2 + Y1 ^ 2)
theta = Atn(Z1 * a / (p * b))
e_dash2 = (a ^ 2 - b ^ 2) / b ^ 2
lat_deg = Atn((Z1 + e_dash2 * b * Sin(theta) ^ 3) / (p - e2 * a * Cos(theta) ^
3)) * 180 / [Link]()
lon_deg = Atn2(Y1, X1) * 180 / [Link]()
VN2000_To_WGS84_DMS = ToDMS(lat_deg) & "," & ToDMS(lon_deg)
End Function
Private Function ToDMS(deg As Double) As String
Dim d As Integer, m As Integer, s As Double
d = Int(Abs(deg))
m = Int((Abs(deg) - d) * 60)
s = ((Abs(deg) - d) * 60 - m) * 60
ToDMS = Format(d, "000") & "d" & Format(m, "00") & "m" & Format(s,
"00.0000000000") & "s"
End Function
Private Function Atn2(y As Double, x As Double) As Double
If x > 0 Then
Atn2 = Atn(y / x)
ElseIf x < 0 And y >= 0 Then
Atn2 = Atn(y / x) + [Link]()
ElseIf x < 0 And y < 0 Then
Atn2 = Atn(y / x) - [Link]()
ElseIf x = 0 And y <> 0 Then
Atn2 = [Link]() / 2 * Sgn(y)
Else
Atn2 = 0
End If
End Function
