http://www.gps-forum.ru/cgi-bin/forum/showpost.pl?Board=gpsgeneral&Number=9972

sk42-wgs84-gost-p-51794-2001.vb

' Функции рабочей книги для преобразования
' геодезических координат из координатной системы Пулково 1942
' в координатную систему WGS84

' Все угловые значения передаются и возвращаются в градусах,
' высоты передаются и возвращаются в метрах

' Сделано на основании ГОСТ Р 51794-2001 Российской Федерации

Const Pi As Double = 3.14159265358979 ' Число Пи
Const ro As Double = 206264.8062 ' Число угловых секунд в радиане

' Эллипсоид Красовского (Пулково 1942)
Const aP As Double = 6378245 ' Большая полуось
Const alP As Double = 1 / 298.3 ' Сжатие
Const e2P As Double = 2 * alP - alP ^ 2 ' Квадрат эксцентриситета

' Эллипсоид GRS80 (WGS84)
Const aW As Double = 6378137 ' Большая полуось
Const alW As Double = 1 / 298.257223563 ' Сжатие
Const e2W As Double = 2 * alW - alW ^ 2 ' Квадрат эксцентриситета

' Вспомогательные значения для преобразования эллипсоидов
Const a As Double = (aP + aW) / 2
Const e2 As Double = (e2P + e2W) / 2
Const da As Double = aW - aP
Const de2 As Double = e2W - e2P

' Линейные элементы трансформирования, в метрах
Const dx As Double = 23.92
Const dy As Double = -141.27
Const dz As Double = -80.9
' Угловые элементы трансформирования, в секундах
Const wx As Double = 0
Const wy As Double = -0.35
Const wz As Double = -0.82
' Дифференциальное различие масштабов
Const ms As Double = -0.12 * 10 ^ -6

Function WGS84Lat(Bd, Ld, H As Double) As Double
Dim B, L, M, N, dB As Double

B = Bd * Pi / 180
L = Ld * Pi / 180
M = a * (1 - e2) / (1 - e2 * Sin(B) ^ 2) ^ 1.5
N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5

dB = ro / (M + H) * (N / a * e2 * Sin(B) * Cos(B) * da _
+ (N ^ 2 / a ^ 2 + 1) * N * Sin(B) * Cos(B) * de2 / 2 _
- (dx * Cos(L) + dy * Sin(L)) * Sin(B) + dz * Cos(B)) _
- wx * Sin(L) * (1 + e2 * Cos(2 * B)) _
+ wy * Cos(L) * (1 + e2 * Cos(2 * B)) _
- ro * ms * e2 * Sin(B) * Cos(B)

WGS84Lat = Bd + dB / 3600
End Function

Function WGS84Long(Bd, Ld, H As Double) As Double
Dim B, L, N, dL As Double

B = Bd * Pi / 180
L = Ld * Pi / 180
N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5

dL = ro / ((N + H) * Cos(B)) * (-dx * Sin(L) + dy * Cos(L)) _
+ Tan(B) * (1 - e2) * (wx * Cos(L) + wy * Sin(L)) - wz

WGS84Long = Ld + dL / 3600
End Function

Function WGS84Alt(Bd, Ld, H As Double) As Double
Dim B, L, N, dH As Double

B = Bd * Pi / 180
L = Ld * Pi / 180
N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5

dH = -a / N * da + N * Sin(B) ^ 2 * de2 / 2 _
+ (dx * Cos(L) + dy * Sin(L)) * Cos(B) + dz * Sin(B) _
- N * e2 * Sin(B) * Cos(B) * (wx / ro * Sin(L) - wy / ro * Cos(L)) _
+ (a ^ 2 / N + H) * ms

WGS84Alt = H + dH
End Function