Public Class ConvertGPS
    Inherits CollectionBase

    'Constants
    Private Const k0 = 0.9996

    'Private Variables
    Private a, b As Double
    Private _myDatum As MyDatums

    'Public Variables
    Public dtm(10) As Datum

    '**************'
    '  STRUCTURES  '
    '**************'

    Public Structure UTM
        Dim Northing As Double
        Dim Easting As Double
        Dim zone As Integer
    End Structure

    Public Structure XY
        Dim X As Double
        Dim Y As Double
    End Structure

    Public Structure LatLon
        Dim Latitude As String
        Dim Longitude As String
    End Structure

    Public Structure Datum
        Dim Name As String
        Dim a As Double
        Dim b As Double
    End Structure

    'Datum Type
    Public Enum MyDatums
        NAD83 = 0
        GRS80 = 1
        WGS72 = 2
        Austrailain1965 = 3
        Krasovsky1940 = 4
        International1924 = 5
        Clake1880 = 6
        Clarke1866 = 7
        Airy1830 = 8
        Bessel1841 = 9
        Everest1830 = 10
    End Enum

    Public Sub New()
        'NAD83/WGS83
        dtm(0).Name = "NAD83/WGS84 (Global)"
        dtm(0).a = 6378137
        dtm(0).b = 6353752.3142
        'GRS 80
        dtm(1).Name = "GRS 80 (US)"
        dtm(1).a = 6378137
        dtm(1).b = 6356752.3141
        'WGS72
        dtm(2).Name = "WGS72 (NASA, DOD)"
        dtm(2).a = 6378135
        dtm(2).b = 6356750.5
        'Australian 1965
        dtm(3).Name = "Australian 1965 (Australia)"
        dtm(3).a = 6378160
        dtm(3).b = 6356774.7
        'Krasovsky 1940
        dtm(4).Name = "Krasovsky 1940 (Soviet Union)"
        dtm(4).a = 6378245
        dtm(4).b = 6356863.0
        'International (1924) - Hayford (1909)
        dtm(5).Name = "International (1924) - Hayford (1909) (Global Except as listed)"
        dtm(5).a = 6378388
        dtm(5).b = 6356911.9
        'Clake 1880
        dtm(6).Name = "Clake 1880 (France, Africa)"
        dtm(6).a = 6378249.1
        dtm(6).b = 6356514.9
        'Clarke 1866
        dtm(7).Name = "Clarke 1866 (North America)"
        dtm(7).a = 6378206.4
        dtm(7).b = 6356583.8
        'Airy 1830 
        dtm(8).Name = "Airy 1830 (Great Britian)"
        dtm(8).a = 6377563.4
        dtm(8).b = 66356256.9
        'Bessel 1841
        dtm(9).Name = "Bessel 1841 (Centreal Europe, Chile, Indonesia)"
        dtm(9).a = 6377397.2
        dtm(9).b = 6356079.0
        'Everest 1830
        dtm(10).Name = "Everest 1830 (South Asia)"
        dtm(10).a = 6377276.3
        dtm(10).b = 6353075.4
    End Sub

    'Property
    Public Property datumType() As MyDatums
        Get
            Return _myDatum
        End Get
        Set(ByVal value As MyDatums)
            _myDatum = value
            Me.a = dtm(_myDatum).a
            Me.b = dtm(_myDatum).b
        End Set
    End Property

    Public Function DMSToDecimal(ByVal dd As Double, ByVal mm As Double, Optional ByVal ss As Double = 0) As Double
        Return (dd + mm / 60) + (ss / 3600)
    End Function

    Public Function DecimalToDMS(ByVal decimalValue As Double) As Integer()
        Dim intDMS(2) As Integer
        Dim dblRemainder, dblMM As Double
        Dim intDD, intMM, intSS As Integer
        intDD = Math.Floor(decimalValue)
        dblRemainder = decimalValue - intDD
        dblMM = dblRemainder * 60
        intMM = Math.Floor(dblMM)
        dblRemainder = dblMM - intMM
        intSS = Math.Round(dblRemainder * 60)
        intDMS(0) = intDD
        intDMS(1) = intMM
        intDMS(2) = intSS
        Return intDMS
    End Function

    Private Function DegreesToRadians(ByVal degrees As Double) As Double
        Dim radians As Double
        radians = (degrees / 180 * Math.PI)
        Return radians
    End Function

    Private Function RadiansToDegrees(ByVal Radians As Double) As Double
        Dim degrees As Double
        degrees = (Radians / Math.PI * 180)
        Return degrees
    End Function

    Private Function MeridianArc(ByVal lat As Double) As Double
        Dim alpha, beta, gamma, delta, epsilon, n, result As Double

        'Calculate n
        n = (a - b) / (a + b)
        'Calculate alpha
        alpha = ((a + b) / 2) * (1 + ((n ^ 2) / 4) + ((n ^ 4) / 64))
        'Calculate beta
        beta = (-3 * n / 2) + (9 * (n ^ 3) / 16) + (-3 * (n ^ 5) / 32)
        'calculate gamma
        gamma = (15 * (n ^ 2) / 16) + (-15 * (n ^ 4) / 32)
        'calculate Delta
        delta = (-35.0 * (n ^ 3.0) / 48.0) + (105.0 * (n ^ 5.0) / 256.0)
        'calculate epsilon
        epsilon = (315.0 * (n ^ 4.0) / 512.0)
        result = alpha * (lat + (beta * Math.Sin(2 * lat)) + (gamma * Math.Sin(4 * lat)) + (delta * Math.Sin(6 * lat)) + (epsilon * Math.Sin(8 * lat)))
        Return result

    End Function

    Private Function UTMCentralMeridian(ByVal Zone As Integer) As Double
        Dim cMeridian As Double
        Dim cm As Double = (Zone * 6) - 183
        cMeridian = DegreesToRadians(cm)
        Return cMeridian
    End Function

    Private Function LatitudeFootPoint(ByVal lat As Double) As Double
        Dim y, alpha, beta, gamma, delta, epsilon, n, result As Double
        'Calculate n
        n = (a - b) / (a + b)
        'calculate alpha
        alpha = ((a + b) / 2) * (1 + ((n ^ 2) / 4) + ((n ^ 4) / 64))
        'Calculate y
        y = lat / alpha
        'Calculate beta
        beta = (3.0 * n / 2.0) + (-27.0 * (n ^ 3.0) / 32.0) + (269.0 * (n ^ 5.0) / 512.0)
        'Calculate gamma
        gamma = (21.0 * (n ^ 2.0) / 16.0) + (-55.0 * (n ^ 4.0) / 32.0)
        'Calculate delta
        delta = (151.0 * (n ^ 3.0) / 96.0) + (-417.0 * (n ^ 5.0) / 128.0)
        'Calculate epsilon
        epsilon = (1097.0 * (n ^ 4.0) / 512.0)
        'Calculate the Result
        result = y + (beta * Math.Sin(2.0 * y)) + (gamma * Math.Sin(4.0 * y)) + (delta * Math.Sin(6.0 * y)) + (epsilon * Math.Sin(8.0 * y))

        Return result
    End Function

    Private Function LatitudeLongitudeToXY(ByVal lat As Double, ByVal longitude As Double, ByVal CentralMeridianLong As Double) As XY
        Dim xy As XY
        Dim tmp, n, nu2, ep2, t, t2, l, coef13, coef14, coef15, coef16, coef17, coef18 As Double

        'Calculate ep2
        ep2 = ((a ^ 2.0) - (b ^ 2.0)) / (b ^ 2.0)
        'Calculate nu2
        nu2 = ep2 * ((Math.Cos(lat)) ^ 2)
        'Calculate n
        n = Math.Pow(a, 2) / (b * Math.Sqrt(1 + nu2))
        'Calculate t
        t = Math.Tan(lat)
        t2 = t * t
        tmp = (t2 * t2 * t2) - (t ^ 6)
        'Calculate l
        l = longitude - CentralMeridianLong
        'Calculate the Coefficioents
        coef13 = 1 - t2 + nu2
        coef14 = 5 - t2 + 9 * nu2 + 4 * (nu2 * nu2)
        coef15 = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - 58.0 * t2 * nu2
        coef16 = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - 330.0 * t2 * nu2
        coef17 = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2)
        coef18 = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2)
        'Calculate Easting
        xy.X = n * Math.Cos(lat) * l + (n / 6.0 * ((Math.Cos(lat)) ^ 3) * coef13 * (l ^ 3)) _
            + (n / 120.0 * ((Math.Cos(lat)) ^ 5.0) * coef15 * (l ^ 5.0)) _
            + (n / 5040.0 * ((Math.Cos(lat)) ^ 7.0) * coef17 * (l ^ 7.0))
        'Calculate Northing
        xy.Y = MeridianArc(lat) + (t / 2.0 * n * ((Math.Cos(lat)) ^ 2.0) * (l ^ 2.0)) _
            + (t / 24.0 * n * ((Math.Cos(lat)) ^ 4.0) * coef14 * (l ^ 4.0)) _
            + (t / 720.0 * n * ((Math.Cos(lat)) ^ 6.0) * coef16 * (l ^ 6.0)) _
            + (t / 40320.0 * n * ((Math.Cos(lat)) ^ 8.0) * coef18 * (l ^ 8.0))

        Return xy
    End Function

    Private Function XYToLatitudeLongitude(ByVal xy As XY, ByVal CentralMeridianLong As Double) As LatLon
        Dim latlon As LatLon
        Dim phif, Nf, Nfpow, nuf2, ep2, tf, tf2, tf4, cf, x1frac, x2frac, x3frac, x4frac As Double
        Dim x5frac, x6frac, x7frac, x8frac, x2poly, x3poly, x4poly, x5poly, x6poly, x7poly, x8poly As Double

        'get phif
        phif = LatitudeFootPoint(xy.Y)
        'Calculate ep2
        ep2 = ((a ^ 2.0) - (b ^ 2.0)) / (b ^ 2.0)
        'Precalculate cos (phif)
        cf = Math.Cos(phif)
        'Precalculate nuf2
        nuf2 = ep2 * (cf ^ 2.0)
        'Precalculate Nf and initialize Nfpow 
        Nf = Math.Pow(a, 2.0) / (b * Math.Sqrt(1 + nuf2))
        Nfpow = Nf
        'Precalculate tf 
        tf = Math.Tan(phif)
        tf2 = tf * tf
        tf4 = tf2 * tf2

        x1frac = 1.0 / (Nfpow * cf)

        Nfpow *= Nf   'now equals Nf**2)
        x2frac = tf / (2.0 * Nfpow)

        Nfpow *= Nf   'now equals Nf**3)
        x3frac = 1.0 / (6.0 * Nfpow * cf)

        Nfpow *= Nf   'now equals Nf**4)
        x4frac = tf / (24.0 * Nfpow)

        Nfpow *= Nf   'now equals Nf**5)
        x5frac = 1.0 / (120.0 * Nfpow * cf)

        Nfpow *= Nf   'now equals Nf**6)
        x6frac = tf / (720.0 * Nfpow)

        Nfpow *= Nf   'now equals Nf**7)
        x7frac = 1.0 / (5040.0 * Nfpow * cf)

        Nfpow *= Nf   'now equals Nf**8)
        x8frac = tf / (40320.0 * Nfpow)

        'Precalculate polynomial coefficients for x**n.
        'x**1 does not have a polynomial coefficient.
        x2poly = -1.0 - nuf2

        x3poly = -1.0 - 2 * tf2 - nuf2

        x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 * nuf2) - 9.0 * tf2 * (nuf2 * nuf2)

        x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2

        x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2

        x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2)

        x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2)

        'Calculate latitude
        latlon.Latitude = phif + x2frac * x2poly * (xy.X * xy.X) _
                        + x4frac * x4poly * (xy.X ^ 4.0) _
                        + x6frac * x6poly * (xy.X ^ 6.0) _
                        + x8frac * x8poly * (xy.X ^ 8.0)

        'Calculate longitude
        latlon.Longitude = CentralMeridianLong + x1frac * xy.X _
                        + x3frac * x3poly * (xy.X ^ 3.0) _
                        + x5frac * x5poly * (xy.X ^ 5.0) _
                        + x7frac * x7poly * (xy.X ^ 7.0)
        Return latlon
    End Function

    Public Function LatitudeLongitudeToUTM(ByVal latitude As Double, ByVal longitude As Double, ByVal zone As Integer) As UTM

        latitude = DegreesToRadians(latitude)
        longitude = DegreesToRadians(longitude)
        Dim xy As XY = LatitudeLongitudeToXY(latitude, longitude, UTMCentralMeridian(zone))
        Dim utm As UTM
        xy.X = xy.X * k0 + 500000.0
        xy.Y = xy.Y * k0
        If xy.Y < 0.0 Then
            xy.Y = xy.Y + 10000000.0
        End If
        utm.Northing = xy.Y
        utm.Easting = xy.X
        utm.zone = zone
        Return utm

    End Function

End Class