Imports ASCOM.Utilities

Module JDadjust

    Dim J2000 As Double = 2451545.0
    Dim JD As Double
    Dim Now As Date
    Public Yr As Integer
    Dim Mo As Integer
    Dim Da As Integer
    Public NewRaH As String
    Public NewRaM As String
    Public NewRaS As String
    Public NewDecD As String
    Public NewDecM As String
    Public NewDecS As String
    Public NewRA As String
    Public NewDEC As String
    Public RaIn As String
    Public DecIn As String
    Dim utl As New ASCOM.Utilities.Util

    Public Function JDate(ByVal Y As Integer, ByVal M As Integer, ByVal D As Single)
        Dim JD As Double
        Dim A As Double
        Dim B As Double

        If M <= 2 Then
            Y = Y - 1
            M = M + 12
        End If

        A = Math.Truncate(Y / 100)
        B = 2 - A + Math.Truncate(A / 4)
        Dim Y1 As Integer = Y + 4716
        Dim M1 As Integer = 30.6001 * (M + 1)

        JD = Math.Truncate(365.25 * (Y + 4716)) + Math.Truncate(30.6001 * (M + 1)) + D + B - 1524.5

        Return JD

    End Function

    Public Sub RaDecJ2000()
        ' Convert RA and DEC from J2000 to Current Date JD

        Now = Date.Now
        Yr = Now.Year
        Mo = Now.Month
        Da = Now.Day
        JD = JDate(Yr, Mo, Da)

        Dim Ra As String = RaIn
        Dim Dec As String = DecIn

        REM Input starting epoch
        Dim JD0 As Double = J2000

        REM Input RA and Dec
        Dim ra0 As Double = 11.0 + (38.0 / 60.0) + (24.96 / 3600.0)
        Dim dec0 As Double = 89.0 + (22.0 / 60.0) + (15.1 / 3600.0)

        ra0 = RAtoHours(Mid(Ra, 1, 2), Mid(Ra, 4, 2), Mid(Ra, 7, 3))
        dec0 = DECtoDegs(Mid(Dec, 1, 1), Mid(Dec, 2, 2), Mid(Dec, 5, 2), Mid(Dec, 8, 2))
        REM override test parameters using inputs to function

        REM Find the transformation matrix
        Dim T0, T, Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz As Double
        T0 = (JD0 - 2433282.423) / 36524.219878
        T = (JD - JD0) / 36524.219878
        Xx = 1.0 - ((29696.0 + 26.0 * T0) * T * T - 13.0 * T * T * T) * 0.00000001
        Xy = ((2234941.0 + 1355.0 * T0) * T - 676.0 * T * T + 221.0 * T * T * T) * 0.00000001
        Xz = ((971690.0 - 414.0 * T0) * T + 207.0 * T * T + 96.0 * T * T * T) * 0.00000001
        Yx = -Xy
        Yy = 1.0 - ((24975.0 + 30.0 * T0) * T * T - 15.0 * T * T * T) * 0.00000001
        Yz = -((10858.0 + 2.0 * T0) * T * T) * 0.00000001
        Zx = -Xz
        Zy = Yz
        Zz = 1.0 - ((4721.0 - 4.0 * T0) * T * T) * 0.00000001

        REM Transform coordinates
        Dim pi, rar0, decr0, x0, y0, z0, x, y, z, rarad, decrad As Double
        pi = 4.0 * Math.Atan(1.0)
        rar0 = 15.0 * pi / 180.0 * ra0
        decr0 = pi / 180.0 * dec0
        x0 = Math.Cos(decr0) * Math.Cos(rar0)
        y0 = Math.Cos(decr0) * Math.Sin(rar0)
        z0 = Math.Sin(decr0)
        x = Xx * x0 + Yx * y0 + Zx * z0
        y = Xy * x0 + Yy * y0 + Zy * z0
        z = Xz * x0 + Yz * y0 + Zz * z0
        rarad = Math.Atan(y / x)
        If x < 0.0 Then rarad = rarad + pi
        If (y < 0.0) And (x > 0.0) Then rarad = rarad + 2.0 * pi
        decrad = Math.Atan(z / Math.Sqrt(1.0 - z * z))

        REM Print transformed coordinates in DMS format
        Dim ra1, d1, dec1 As Double
        ra1 = 180.0 / pi * rarad
        ra1 = ra1 / 15.0
        d1 = Math.Abs(ra1)

        NewRA = utl.DegreesToDMS(d1, ":", ":", "", 1)

        dec1 = 180.0 / pi * decrad
        d1 = Math.Abs(dec1)
        Dim Dsign As String = "+"
        If dec1 < 0 Then Dsign = "-"

        NewDEC = Dsign & utl.DegreesToDMS(d1, "*", ":", "", 0)

    End Sub

    Public Sub RaDecToJ2000()
        ' Convert Current mount RA and DEC to J2000

        Now = Date.Now
        Yr = Now.Year
        Mo = Now.Month
        Da = Now.Day
        JD = J2000

        Dim Ra As String = RaIn
        Dim Dec As String = DecIn

        REM Input starting epoch
        Dim JD0 As Double = JDate(Yr, Mo, Da)

        REM Input RA and Dec
        Dim ra0 As Double = 11.0 + (38.0 / 60.0) + (24.96 / 3600.0)
        Dim dec0 As Double = 89.0 + (22.0 / 60.0) + (15.1 / 3600.0)

        ra0 = RAtoHours(Mid(Ra, 1, 2), Mid(Ra, 4, 2), Mid(Ra, 7, 3))
        dec0 = DECtoDegs(Mid(Dec, 1, 1), Mid(Dec, 2, 2), Mid(Dec, 5, 2), Mid(Dec, 8, 2))
        REM override test parameters using inputs to function

        REM Find the transformation matrix
        Dim T0, T, Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz As Double
        T0 = (JD0 - 2433282.423) / 36524.219878
        T = (JD - JD0) / 36524.219878
        Xx = 1.0 - ((29696.0 + 26.0 * T0) * T * T - 13.0 * T * T * T) * 0.00000001
        Xy = ((2234941.0 + 1355.0 * T0) * T - 676.0 * T * T + 221.0 * T * T * T) * 0.00000001
        Xz = ((971690.0 - 414.0 * T0) * T + 207.0 * T * T + 96.0 * T * T * T) * 0.00000001
        Yx = -Xy
        Yy = 1.0 - ((24975.0 + 30.0 * T0) * T * T - 15.0 * T * T * T) * 0.00000001
        Yz = -((10858.0 + 2.0 * T0) * T * T) * 0.00000001
        Zx = -Xz
        Zy = Yz
        Zz = 1.0 - ((4721.0 - 4.0 * T0) * T * T) * 0.00000001

        REM Transform coordinates
        Dim pi, rar0, decr0, x0, y0, z0, x, y, z, rarad, decrad As Double
        pi = 4.0 * Math.Atan(1.0)
        rar0 = 15.0 * pi / 180.0 * ra0
        decr0 = pi / 180.0 * dec0
        x0 = Math.Cos(decr0) * Math.Cos(rar0)
        y0 = Math.Cos(decr0) * Math.Sin(rar0)
        z0 = Math.Sin(decr0)
        x = Xx * x0 + Yx * y0 + Zx * z0
        y = Xy * x0 + Yy * y0 + Zy * z0
        z = Xz * x0 + Yz * y0 + Zz * z0
        rarad = Math.Atan(y / x)
        If x < 0.0 Then rarad = rarad + pi
        If (y < 0.0) And (x > 0.0) Then rarad = rarad + 2.0 * pi
        decrad = Math.Atan(z / Math.Sqrt(1.0 - z * z))

        REM Print transformed coordinates in DMS format
        Dim ra1, d1, dec1 As Double
        ra1 = 180.0 / pi * rarad
        ra1 = ra1 / 15.0
        d1 = Math.Abs(ra1)

        NewRA = utl.DegreesToDMS(d1, ":", ":", "", 1)

        dec1 = 180.0 / pi * decrad
        d1 = Math.Abs(dec1)
        Dim Dsign As String = "+"
        If dec1 < 0 Then Dsign = "-"

        NewDEC = Dsign & utl.DegreesToDMS(d1, "*", ":", "", 0)

    End Sub


    Function RAtoHours(ByVal RaH As Integer, ByVal RaM As Integer, ByVal RaS As Single)
        Dim Dec As Double
        Dec = RaH + (RaM / 60) + (RaS / 3600)
        Return Dec

    End Function

    Function DECtoDegs(ByVal DegSign As String, ByVal DegD As Integer, ByVal DecM As Integer, ByVal DecS As Single)
        Dim Dec As Double
        Dim SignDec As Integer = 1
        If DegSign = "-" Then
            SignDec = -1
        End If
        Dec = (Math.Abs(DegD) + (DecM / 60) + (DecS / 3600)) * SignDec
        Return Dec

    End Function


End Module
