VBA lat long to UTM conversion

to_utm.vba

Function to_utm(x, y) As Collection
Dim coordinates As Collection
Set coordinates = New Collection

E = 0.00669438
R = 6378137

K0 = 0.9996

E2 = E * E
E3 = E2 * E
E_P2 = E / (1 - E)

SQRT_E = (1 - E) ^ (1 / 2)
E_1 = (1 - SQRT_E) / (1 + SQRT_E)
E_2 = E_1 * E_1
E_3 = E_2 * E_1
E_4 = E_3 * E_1
E_5 = E_4 * E_1

M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256)
M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024)
M3 = (15 * E2 / 256 + 45 * E3 / 1024)
M4 = (35 * E3 / 3072)

P2 = (3 / 2 * E_1 - 27 / 32 * E_3 + 269 / 512 * E_5)
P3 = (21 / 16 * E_2 - 55 / 32 * E_4)
P4 = (151 / 96 * E_3 - 417 / 128 * E_5)
P5 = (1097 / 512 * E_4)

Pi = Round((355 / 113), 6)
zone_number = Int(((x + 180) / 6) + 1)

lat_rad = (Pi / 180) * y
lat_sin = Sin(lat_rad)
lat_cos = Cos(lat_rad)
lat_tan = lat_sin / lat_cos
lat_tan2 = lat_tan * lat_tan
lat_tan4 = lat_tan2 * lat_tan2

lon_rad = (Pi / 180) * x
central_lon = ((zone_number - 1) * 6) - 180 + 3
central_lon_rad = (Pi / 180) * central_lon

N = R / (1 - E * lat_sin ^ 2) ^ (1 / 2)
c = E_P2 * lat_cos ^ 2

a = lat_cos * (lon_rad - central_lon_rad)
a2 = a * a
a3 = a2 * a
a4 = a3 * a
a5 = a4 * a
a6 = a5 * a

m = R * (M1 * lat_rad - M2 * Sin(2 * lat_rad) + M3 * Sin(4 * lat_rad) - M4 * Sin(6 * lat_rad))

easting = (K0 * N * (a + a3 / 6 * (1 - lat_tan2 + c) + a5 / 120 * (5 - 18 * lat_tan2 * lat_tan4 + 72 * c - 58 * E_P2)) + 500000)
northing = (K0 * (m + N * lat_tan * (a2 / 2 + a4 / 24 * (5 - lat_tan2 + 9 * c + 4 * c ^ 2) + a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2))))

coordinates.Add easting
coordinates.Add northing
coordinates.Add zone_number
Set to_utm = coordinates

End Function

Sub conversion()
Dim conversion As Collection

x = -122.0230604
y = 38.502032

Set conversion = to_utm(x, y)

Debug.Print conversion(1), conversion(2), conversion(3)

End Sub