Option Explicit On

Public Class DomeControl

    REM Verson 5.0.7 25-Mar-2013


    '---------------------------------------------------------------------
    ' Copyright © 2012 Chris Rowland
    '
    ' Permission is hereby granted to use this Software for any purpose
    ' including combining with commercial products, creating derivative
    ' works, and redistribution of source or binary code, without
    ' limitation or consideration. Any redistributed copies of this
    ' Software must include the above Copyright Notice.
    '
    ' THIS SOFTWARE IS PROVIDED "AS IS". I MAKE NO
    ' WARRANTIES REGARDING THIS SOFTWARE, EXPRESS OR IMPLIED, AS TO ITS
    ' SUITABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
    '---------------------------------------------------------------------
    '
    '   ==============
    '   DomeControl.cls
    '   ==============
    '
    ' Written:  2012/05/06   Chris Rowland <chris.rowland@dsl.pipex.com>
    '
    ' This class gets the dome pointing direction using Paul Bourke's
    ' algorithim for determining the intersections of a sphere and a line.
    '
    ' The line starts at the scope OTA position and goes in the direction
    ' the scope is looking.
    ' The sphere represents the dome, there should be one intersection which
    ' is where the scope pointing direction intercepts the sphere.
    '
    ' Edits:
    '
    ' When       Who     What
    ' ---------  ---     --------------------------------------------------
    ' 2012/05/06 cdr     Initial edit
    ' -----------------------------------------------------------------------------

    ' initialisation variables
    Private m_mountPosn As XYZ
    Private m_domeRadius As Double
    Private m_otaOffset As Double
    Private m_otaOffsetEW As Double
    Private m_latitude As Double
    Private m_isGem As Boolean

    Private Const PI As Double = 3.1415926535

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

    Private Structure XYZ
        Public X As Double
        Public Y As Double
        Public z As Double
    End Structure

    Private m_domeAltitude As Double

    ''' <summary>
    ''' Set the dome and mount properties, these only need to be done once.
    ''' The units for the distances do not matter but must be the same for all items.
    ''' </summary>
    ''' <param name="DomeRadius">The dome radius.</param>
    ''' <param name="OtaOffset">The ota offset.</param>
    ''' <param name="mountOffsetN">The mount offset N.</param>
    ''' <param name="mountoffsetE">The mountoffset E.</param>
    ''' <param name="mountOffsetUp">The mount offset up.</param>
    ''' <param name="latitude">The latitude in degrees</param>
    ''' <param name="isGem">if set to <c>true</c> [is gem].</param>
    Public Sub InitDome(ByVal DomeRadius As Double, _
    ByVal OtaOffset As Double, _
    ByVal OtaOffsetEW As Double, _
    ByVal mountOffsetN As Double, _
    ByVal mountoffsetE As Double, _
    ByVal mountOffsetUp As Double, _
    ByVal Latitude As Double, _
    ByVal isGem As Boolean)
        m_domeRadius = DomeRadius
        m_otaOffset = OtaOffset
        m_otaOffsetEW = OtaOffsetEW
        m_mountPosn = SetXYZ(mountOffsetN, mountoffsetE, mountOffsetUp)
        m_latitude = Latitude
        m_isGem = isGem
    End Sub

    ''' <summary>
    ''' Gets the dome azimuth that is required so the scope can see out using the
    ''' scope pointing direction, the mount parameters and the dome radius
    ''' </summary>
    ''' <param name="scopeAzimuth">The scope azimuth in degrees</param>
    ''' <param name="scopeAltitude">The scope altitude in degrees</param>
    ''' <param name="hourAngle">The hour angle in degrees</param>
    ''' <param name="isPierWest">if set to <c>true</c> [is pier west].</param>
    ''' <returns>the dome azimuth in degrees</returns>
    Public Function DomeAzimuth(ByVal scopeAzimuth As Double, ByVal scopeAltitude As Double, ByVal hourAngle As Double, ByVal isPierWest As Boolean) As Double
        Dim scopePos As XY
        Dim domePos As XY
        Dim otaPosn As XYZ

        ' get the ota position WRT the centre of the dome, allowing for the offsets
        ' the scope direction and the latitude
        otaPosn = OtaPosition(hourAngle, isPierWest)

        ' set the scope alt azm in degrees
        scopePos = SetXY(scopeAzimuth, scopeAltitude)

        domePos = DomeAltAzm(scopePos, otaPosn)
        ' set the altitude to the range -180 to +180
        Do While domePos.Y >= 180.0#
            domePos.Y = domePos.Y - 360.0#
        Loop
        Do While domePos.Y < -180.0#
            domePos.Y = domePos.Y + 360.0#
        Loop
        ' save this so it can be read by the DomeAltitude Function
        m_domeAltitude = domePos.Y

        ' set the azimuth to the range 0 to 360 degrees
        Do While (domePos.X >= 360.0#)
            domePos.X = domePos.X - 360.0#
        Loop
        Do While (domePos.X < 0)
            domePos.X = domePos.X + 360.0#
        Loop
        DomeAzimuth = domePos.X
        Debug.Print("DomeAzimuth Azm:" & scopeAzimuth & ", Alt:" & scopeAltitude & ", Ha:" & hourAngle & ", isPierWest:" & isPierWest & " = " & DomeAzimuth)
    End Function

    ''' <summary>
    ''' Get the dome altitude in degrees.
    ''' The altitude is set by the DomeAzimuth function
    ''' so this must only be called immediately after this.
    ''' </summary>
    ''' <returns></returns>
    Public Function DomeAltitude() As Double
        DomeAltitude = m_domeAltitude
    End Function

    ''' <summary>
    ''' Return the OTA position wrt the dome centre
    ''' X is to the N,
    ''' Y is to the W,
    ''' Z is up
    ''' Currently assumes the Nothern hemisphere
    ''' </summary>
    ''' <param name="ha"></param>
    ''' <param name="isPierWest"></param>
    ''' <returns></returns>
    Private Function OtaPosition(ByVal ha As Double, ByVal isPierWest As Boolean) As XYZ
        ' get the position of the mount
        Dim p1 As XYZ
        p1 = m_mountPosn
        p1.Y = -p1.Y        ' change the sign of the y direction to match the WPF coordinate system

        Dim os As XYZ
        os = SetXYZ(0, m_otaOffset, 0)      ' ota offset is along the dec axis
        os = SetXYZ(m_otaOffsetEW, m_otaOffset, 0)      ' ota offset is along the dec axis and east / west offset

        If (m_otaOffset > 0.00001) Then
            If (m_isGem And isPierWest) Then
                ' reverse the dec ota offset
                os.Y = -os.Y
            End If
            'os.Y = -os.Y
            ' the order of the rotations matters
            ' rotate by dec
            '        os = Rotate(os, SetXYZ(0, 1, 0), haDec.Y)
            ' rotate by Ha
            os = Rotate(os, SetXYZ(-1, 0, 0), ha)
            ' rotate by latitude
            os = Rotate(os, SetXYZ(0, -1, 0), Math.Abs(m_latitude))
        End If
        OtaPosition = SetXYZ(p1.X + os.X, p1.Y + os.Y, p1.z + os.z)
        Debug.Print("OtaPosition X:" & OtaPosition.X & ", Y:" & OtaPosition.Y & ", Z:" & OtaPosition.z)
    End Function

    ''' <summary>
    ''' Determine the position of the dome opening using Paul Bourke's intersection of
    ''' a ray and a sphere code.  The ray represents the scope pointing direction and the
    ''' sphere the dome.  The OTA offset is determined from the mount offset and the
    ''' scope pointing direction specifies the direction of the ray.
    ''' </summary>
    ''' <param name="scopeAltAz">scope pointing direction in Alt(y) Azm(x), as radians</param>
    ''' <param name="otaPosn">OTA position wrt the centre of the dome, allowing for the Ha and Dec and the pier side</param>
    ''' <returns>
    ''' Dome Altitude(y), Azimuth(x), in degrees
    ''' </returns>
    Private Function DomeAltAzm(ByVal scopeAltAz As XY, ByVal otaPosn As XYZ) As XY
        ' get the position of the OTA
        Dim p1 As XYZ
        p1 = otaPosn

        ' get the offset from the OTA position towards the scope pointing direction,
        Dim dp As XYZ
        dp = PolarToCartesian(scopeAltAz, m_domeRadius * 4)

        ' the conversion is simplified because we are using the centre of the dome as the origin
        Dim a As Double
        a = dp.X * dp.X + dp.Y * dp.Y + dp.z * dp.z
        Dim b As Double
        b = 2 * (dp.X * p1.X + dp.Y * p1.Y + dp.z * p1.z)
        Dim c As Double
        c = p1.X * p1.X + p1.Y * p1.Y + p1.z * p1.z
        c = c - m_domeRadius * m_domeRadius
        Dim bb4ac As Double
        bb4ac = b * b - 4 * a * c
        If Math.Abs(a) < 0.000001 Or bb4ac < 0 Then
            Err.Raise(123, "DomeControl", "scope is not pointing out of the dome")
        End If
        Dim mu As Double
        mu = (-b + Math.Sqrt(bb4ac)) / (2 * a)
        If (mu <= 0 Or mu > 1.0#) Then
            mu = (-b - Math.Sqrt(bb4ac)) / (2 * a)
        End If
        ' There are potentially two points of intersection given by
        ' p = p1 + mu1 (p2 - p1)
        ' p = p1 + mu2 (p2 - p1)
        Dim pos As XYZ
        pos = SetXYZ(p1.X + mu * dp.X, p1.Y + mu * dp.Y, p1.z + mu * dp.z)
        DomeAltAzm = CartesianToPolar(pos)
        Debug.Print("DomeAltAzm Azm:" & DomeAltAzm.X & ", Alt:" & DomeAltAzm.Y)
    End Function

    ''' <summary>
    ''' Rotates the specified val.
    ''' </summary>
    ''' <param name="vector">The vector that's rotated</param>
    ''' <param name="axis">A unit vector defining the axis direction, must be in the X or Y direction</param>
    ''' <param name="angle">The rotation angle in degrees</param>
    ''' <returns></returns>
    Private Function Rotate(ByVal vector As XYZ, ByVal axis As XYZ, ByVal angle As Double) As XYZ
        Dim ar As Double
        ar = D2R(angle)     ' convert the angle to radians
        ' choose the axis to rotate about
        If Math.Abs(axis.X) > Math.Abs(axis.Y) Then
            If axis.X > 0 Then ar = -ar
            Rotate.X = vector.X
            Rotate.z = Math.Cos(ar) * vector.z - Math.Sin(ar) * vector.Y
            Rotate.Y = Math.Sin(ar) * vector.z + Math.Cos(ar) * vector.Y
        Else
            If axis.Y < 0 Then ar = -ar
            Rotate.Y = vector.Y
            Rotate.X = Math.Cos(ar) * vector.X - Math.Sin(ar) * vector.z
            Rotate.z = Math.Sin(ar) * vector.X + Math.Cos(ar) * vector.z
        End If
    End Function

    ''' <summary>
    ''' converts angle in degres to radians.
    ''' </summary>
    ''' <param name="angle">The angle in degrees</param>
    ''' <returns>the angle in radians</returns>
    Private Function D2R(ByVal angle As Double) As Double
        D2R = angle / 180.0# * PI
    End Function

    ''' <summary>
    ''' converts angle in radians to degrees
    ''' </summary>
    ''' <param name="rad">The angle in radians</param
    ''' ><returns>the angle in degrees</returns>
    Private Function R2D(ByVal rad As Double) As Double
        R2D = rad * 180.0# / PI
    End Function

    Private Function SetXY(ByVal X As Double, ByVal Y As Double) As XY
        Dim ret As XY
        ret.X = X
        ret.Y = Y
        SetXY = ret
    End Function

    Private Function SetXYZ(ByVal X As Double, ByVal Y As Double, ByVal z As Double) As XYZ
        Dim ret As XYZ
        ret.X = X
        ret.Y = Y
        ret.z = z
        SetXYZ = ret
    End Function

    ''' <summary>
    ''' Converts an alt azm vector in radians and a distance in radius to an xyz coordinate
    ''' </summary>
    ''' <param name="azmAlt">The azm(X) and alt(Y) angles in degrees</param>
    ''' <param name="radius">The radius.</param>
    ''' <returns>XYZ vector, length is radius</returns>
    Private Function PolarToCartesian(ByVal azmAlt As XY, ByVal radius As Double) As XYZ
        Dim azm As Double
        Dim alt As Double
        azm = D2R(azmAlt.X)
        alt = D2R(azmAlt.Y)
        PolarToCartesian = SetXYZ(radius * Math.Cos(alt) * Math.Cos(azm), _
                                  radius * Math.Cos(alt) * Math.Sin(azm), _
                                  radius * Math.Sin(alt))
    End Function

    ''' <summary>
    ''' Converts a vector to azimuth and altitude in degrees
    ''' </summary>
    ''' <param name="vector">The vector.</param>
    ''' <returns>azimuth and altitude in degrees</returns>
    Private Function CartesianToPolar(ByVal vector As XYZ) As XY
        CartesianToPolar = SetXY(R2D(Atn2(vector.Y, vector.X)), _
                                 90 - R2D(Atn2(Math.Sqrt(vector.X * vector.X + vector.Y * vector.Y), vector.z)))
    End Function

    ''' <summary>
    ''' Atn2 function, returns angle in radians
    ''' </summary>
    ''' <param name="Y">The Y.</param>
    ''' <param name="X">The X.</param>
    ''' <returns>angle in radians</returns>
    Private Function Atn2(ByVal Y As Double, ByVal X As Double) As Double
        Select Case X
            Case Is > 0.0#
                Atn2 = Math.Atan(Y / X)
            Case Is < 0.0#
                Atn2 = Math.Atan(Y / X) + PI * IIf(Y > 0, 1, -1)
            Case Else
                Atn2 = (PI / 2.0#) * Math.Sign(Y)
        End Select
    End Function





    REM OLD Version below

    ''---------------------------------------------------------------------
    '' Copyright © 2012 Chris Rowland
    ''
    '' Permission is hereby granted to use this Software for any purpose
    '' including combining with commercial products, creating derivative
    '' works, and redistribution of source or binary code, without
    '' limitation or consideration. Any redistributed copies of this
    '' Software must include the above Copyright Notice.
    ''
    '' THIS SOFTWARE IS PROVIDED "AS IS". I MAKE NO
    '' WARRANTIES REGARDING THIS SOFTWARE, EXPRESS OR IMPLIED, AS TO ITS
    '' SUITABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
    ''---------------------------------------------------------------------
    ''
    ''   ==============
    ''   DomeControl.cls
    ''   ==============
    ''
    '' Written:  2012/05/06   Chris Rowland <chris.rowland@dsl.pipex.com>
    ''
    '' This class gets the dome pointing direction using Paul Bourke's
    '' algorithim for determining the intersections of a sphere and a line.
    ''
    '' The line starts at the scope OTA position and goes in the direction
    '' the scope is looking.
    '' The sphere represents the dome, there should be one intersection which
    '' is where the scope pointing direction intercepts the sphere.
    ''
    '' Edits:
    ''
    '' When       Who     What
    '' ---------  ---     --------------------------------------------------
    '' 2012/05/06 cdr     Initial edit
    '' -----------------------------------------------------------------------------


    'Private m_mountPosn As XYZ
    'Private m_domeRadius As Double
    'Private m_otaOffset As Double

    'Private Const PI As Double = 3.1415926535

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

    'Public Structure XYZ
    '    Public X As Double
    '    Public Y As Double
    '    Public z As Double
    'End Structure

    'Private m_domeAltitude As Double

    '' set the dome and mount properties, these only need to be
    '' done once.
    '' The units for the distances do not matter but must be the same for
    '' all items.
    'Public Sub InitDome(ByVal DomeRadius As Double, _
    'ByVal OtaOffset As Double, _
    'ByVal mountOffsetN As Double, _
    'ByVal mountoffsetE As Double, _
    'ByVal mountOffsetUp As Double)
    '    m_domeRadius = DomeRadius
    '    m_otaOffset = OtaOffset
    '    m_mountPosn = SetXYZ(mountOffsetN, mountoffsetE, mountOffsetUp)
    'End Sub

    '' Gets the dome azimuth that is required so the scope can see out using the
    '' scope pointing direction, the mount parameters and the dome radius
    ''
    '' scopeAzimuth  azimuth of scope in degrees
    '' scopeAltitude altitude of scope in degrees
    '' hourAngle     hour angle in degrees
    '' isPierWest    true is the pier is west of the meridian, don't care for fork mounts
    ''
    '' returns the dome azimuth in degrees
    'Public Function DomeAzimuth(ByVal scopeAzimuth As Double, ByVal scopeAltitude As Double, ByVal hourAngle As Double, ByVal isPierWest As Boolean) As Double
    '    Dim scopePos As XY
    '    Dim domePos As XY
    '    Dim angle As Double

    '    ' convert scope alt azm to radians
    '    scopePos = SetXY(scopeAzimuth * PI / 180.0#, scopeAltitude * PI / 180.0#)

    '    domePos = DomeAltAzm(scopePos, hourAngle * PI / 180.0#, isPierWest)
    '    ' get the altitude and convert to degrees
    '    angle = domePos.Y * 180.0# / PI
    '    Do While angle >= 180.0#
    '        angle = angle - 360.0#
    '    Loop
    '    Do While angle < -180.0#
    '        angle = angle + 360.0#
    '    Loop
    '    ' save this so it can be read by the DomeAltitude Function
    '    m_domeAltitude = angle
    '    ' get the azimuth and convert to degrees
    '    angle = domePos.X * 180.0# / PI
    '    ' force it to the range 0 to 360 degrees
    '    Do While (angle >= 360.0#)
    '        angle = angle - 360.0#
    '    Loop
    '    Do While (angle < 0)
    '        angle = angle + 360.0#
    '    Loop
    '    DomeAzimuth = angle
    'End Function

    '' Get the dome altitude in degrees.
    '' The altitude is set by the DomeAzimuth function
    '' so this must only be called immediately after this.
    'Public Function DomeAltitude() As Double
    '    DomeAltitude = m_domeAltitude
    'End Function

    '''' <summary>
    '''' Determine the position of the dome opening using Paul Bourke's intersection of
    '''' a ray and a sphere code.  The ray represents the scope pointing direction and the
    '''' sphere the dome.  The OTA offset is determined from the mount offset and the
    '''' scope pointing direction specifies the direction of the ray.
    '''' </summary>
    '''' <param name="scopeAltAz">scope pointing direction in Alt(y) Azm(x), as radians</param>
    '''' <param name="scopeHaAxis">The scope hour angle axis in radians</param>
    '''' <param name="isPierWest">if set to <c>true</c> if the scope is on the west side of the pier.</param>
    '''' <returns>
    '''' Dome Altitude(y), Azimuth(x), in radians
    '''' </returns>
    'Private Function DomeAltAzm(ByVal scopeAltAz As XY, ByVal scopeHaAxis As Double, ByVal isPierWest As Boolean) As XY
    '    ' get the position of the OTA
    '    Dim p1 As XYZ
    '    p1 = SetXYZ(m_mountPosn.X, m_mountPosn.Y, m_mountPosn.z)

    '    ' apply the ota offset in the direction of the hour angle
    '    If Math.Abs(m_otaOffset) > 0.000001 Then
    '        Dim os As XYZ
    '        os = PolarToCartesian(SetXY(scopeHaAxis, 0), IIf(isPierWest, m_otaOffset, -m_otaOffset))
    '        p1.X = p1.X + os.Y
    '        p1.Y = p1.Y - os.X
    '        p1.z = p1.z + os.z
    '    End If
    '    ' get the offset from the OTA position towards the scope pointing direction,
    '    Dim dp As XYZ
    '    dp = PolarToCartesian(scopeAltAz, m_domeRadius * 4)

    '    ' the conversion is simplified because we are using the centre of the dome as the origin
    '    Dim a As Double
    '    a = dp.X * dp.X + dp.Y * dp.Y + dp.z * dp.z
    '    Dim b As Double
    '    b = 2 * (dp.X * p1.X + dp.Y * p1.Y + dp.z * p1.z)
    '    Dim c As Double
    '    c = p1.X * p1.X + p1.Y * p1.Y + p1.z * p1.z
    '    c = c - m_domeRadius * m_domeRadius
    '    Dim bb4ac As Double
    '    bb4ac = b * b - 4 * a * c
    '    If Math.Abs(a) < 0.000001 Or bb4ac < 0 Then
    '        Err.Raise(123, "DomeControl", "scope is not pointing out of the dome")
    '    End If
    '    Dim mu As Double
    '    mu = (-b + Math.Sqrt(bb4ac)) / (2 * a)
    '    If (mu <= 0 Or mu > 1.0#) Then
    '        mu = (-b - Math.Sqrt(bb4ac)) / (2 * a)
    '    End If
    '    ' There are potentially two points of intersection given by
    '    ' p = p1 + mu1 (p2 - p1)
    '    ' p = p1 + mu2 (p2 - p1)
    '    Dim pos As XYZ
    '    pos = SetXYZ(p1.X + mu * dp.X, p1.Y + mu * dp.Y, p1.z + mu * dp.z)
    '    DomeAltAzm = CartesianToPolar(pos)
    'End Function

    'Private Function SetXY(ByVal X As Double, ByVal Y As Double) As XY
    '    Dim ret As XY
    '    ret.X = X
    '    ret.Y = Y
    '    SetXY = ret
    'End Function

    'Private Function SetXYZ(ByVal X As Double, ByVal Y As Double, ByVal z As Double) As XYZ
    '    Dim ret As XYZ
    '    ret.X = X
    '    ret.Y = Y
    '    ret.z = z
    '    SetXYZ = ret
    'End Function

    'Private Function PolarToCartesian(ByVal polar As XY, ByVal radius As Double) As XYZ
    '    PolarToCartesian = SetXYZ(radius * Math.Cos(polar.Y) * Math.Cos(polar.X), radius * Math.Cos(polar.Y) * Math.Sin(polar.X), radius * Math.Sin(polar.Y))
    'End Function

    'Private Function CartesianToPolar(ByVal cartes As XYZ) As XY
    '    CartesianToPolar = SetXY(Atn2(cartes.Y, cartes.X), Atn2(Math.Sqrt(cartes.X * cartes.X + cartes.Y * cartes.Y), cartes.z))
    'End Function

    'Private Function Atn2(ByVal Y As Double, ByVal X As Double) As Double
    '    Select Case X
    '        Case 0
    '            Atn2 = PI / 2.0# + Math.Sign(Y)
    '        Case Is > 0.0#
    '            Atn2 = Math.Atan(Y / X)
    '        Case Is < 0
    '            Atn2 = Math.Atan(Y / X) + PI
    '    End Select
    'End Function


End Class
