' SofaJplDemoVb2_1.vb
' Demonstration program for SofaJpl DLL v. 2.1.2.
' Paul S. Hirose, 2018 August 26

' DOCUMENTATION
' The Visual Studio Object Browser displays the main SofaJpl documentation. After adding a
' reference to SofaJpl to your project, open Object Browser and expand the SofaJpl namespace.

' revision history
#If False Then
2018-09-09 Added a compatibility mode to duplicate the equatorial, ecliptical, and az/el of JPL
Horizons. The first two agree to 1 mas in my tests. Az/el agrees .0001 deg, which is the Horizons
output precision.
#End If


' The time scales you can use for input.
#Const UTC = 1
#Const UT1 = 2
#Const TT = 3

' Set this to one of the time scales above.
#Const TIMESCALE = UTC

' After the binary ephemeris is created, set this False to eliminate redundant processing of the
' multi megabyte ASCII ephemeris every time the program runs.
#Const MAKE_EPHEMERIS = False

' If #define is un-commented, use the delta T you provide. Otherwise use the SofaJpl delta T model.
' Not used if time scale is UTC.
#Const MANUAL_DELTA_T = False

' If True, use a precession nutation model compatible with JPL Horizons: IAU 1976 / 80
' precession / nutation, no frame bias, apply IERS pole offsets.
' NOTE: To obtain coordinates in the horizontal system or the ITRS, compatible with Horizons, it is
' also necessary to activate the MANUAL_DELTA_T option and supply the correct delta T of date.
' In addition, for compatible coordinates in the horizontal system set the deflection of the
' vertical parameters (xi and eta) to zero.
#Const HORIZONS = False


Imports System

' Everything in SofaJpl is in the HirosePS.SofaJpl namespace.
Imports HirosePS


Module Module1

    ' The leap second table. One is included in the SofaJpl distribution .zip file.
    Dim _leapsecTable = "C:\Users\Stan\Documents\astro\SofaJpl_2_1_2\leapSecTable.txt"

    ' JPL ASCII ephemeris and header. Must be downloaded from JPL. Not needed after they are
    ' converted to a binary ephemeris.
    Dim _asciiEphemeris As String = "C:\Users\Stan\Documents\astro\jpl\ascp2000.422"
    Dim _ephemerisHeader As String = "C:\Users\Stan\Documents\astro\jpl\header.422"

    ' The binary JPL ephemeris to create and use.
    Dim _binaryEphemeris As String = "binp2000_2099.422"

    ' mathematical constants
    Dim _degPerHour As Double = 15.0

    ' an enumeration of the three possible angle output formats
    Enum angleFormat
        D
        DM
        DMS
    End Enum

    ' output format desired by the user
    Dim _format As angleFormat

    ' These control the output resolution of angle and time. The values are computed automatically
    ' from the angle accuracy set by the user.
    Dim _angleResolution As Double      ' resolution units per degree
    Dim _timeResolution As Double       ' resolution units per hour

    ' Number of decimal places in the floating point 'f' format to achieve unit vector rectangular
    ' coordinate precision comparable to the angles.
    Dim _vectorResolution As Integer


    Sub Main()

        ' Set the display format and angle accuracy.
        _format = angleFormat.DMS
        Dim angleAccuracy As Double = SofaJpl.Angle.DmsToRad(0, 0, 0.1)

        ' Set the epoch.
        Dim year As Integer = 2018
        Dim month As Integer = 6
        Dim day As Integer = 6
        Dim hour As Integer = 4
        Dim minute As Integer = 0
        Dim second As Double = 0.0
        Dim julian As Boolean = False    ' true if date is in Julian calendar


        Dim deltaT As SofaJpl.Duration
#If TIMESCALE = utc Then
        Dim ut1MinusUtc As SofaJpl.Duration = SofaJpl.Duration.FromSeconds(0.07522)
#Else
#If MANUAL_DELTA_T Then
        deltaT = SofaJpl.Duration.FromSeconds(69.109)
#Else
        ' Get delta T from the SofaJpl model.
        deltaT = SofaJpl.DeltaT.GetDeltaT(
            New SofaJpl.JulianDate(year, month, day, hour, minute, second, julian))
#End If
#End If

        ' Observer position.
        Dim obsName = "Kitt Peak Observatory"    ' name of topocenter (optional)
        Dim lon As Double = SofaJpl.Angle.DmsToRad(-111, 36, 0.0)    ' east longitude
        Dim lat As Double = SofaJpl.Angle.DmsToRad(31, 57, 48.0)      ' north latitude
        Dim height As Single = 2100.0F    ' meters above ellipsoid

        ' This instance of the Topocenter class will supply the observer position and velocity
        ' with respect to the ICRS.
        Dim obs As SofaJpl.Topocenter = New SofaJpl.Topocenter(lat, lon, height)

        ' Construct an Atmosphere object to apply refraction. The parameters required by the
        ' constructor are type "float", thus the f suffix on the numbers.
        Dim altimiterSetting As Single = 29.9F *
            SofaJpl.Atmosphere.MillibarsPerInchHg     ' millibars
        Dim degC As Single = 12.0F       ' Celsius
        Dim dewPointC As Single = 0.0F   ' Celsius.
        Dim atm As New SofaJpl.Atmosphere(height, degC, altimiterSetting, dewPointC, False)

        ' polar motion (radians)
        Dim poleX As Double = SofaJpl.Angle.DmsToRad(0, 0, 0)
        Dim poleY As Double = SofaJpl.Angle.DmsToRad(0, 0, 0)

        ' deflection of the vertical at the observer (radians)
        Dim xi As Double = SofaJpl.Angle.DmsToRad(0, 0, 0)
        Dim eta As Double = SofaJpl.Angle.DmsToRad(0, 0, 0)


        ' Create a binary ephemeris from JPL ASCII files.
#If MAKE_EPHEMERIS Then
        Console.WriteLine(vbCrLf & "CREATING BINARY JPL EPHEMERIS." & vbCrLf)
        SofaJpl.JplEphemeris.AsciiToBinary(_ephemerisHeader, _asciiEphemeris, _binaryEphemeris)
#End If

        ' Open the binary JPL ephemeris. This loads the entire ephemeris into RAM. Even if the
        ' target body is a star, a solar system ephemeris is required for parallax, aberration,
        ' and light deflection due to solar gravitation.
        Dim eph As New SofaJpl.JplEphemeris(_binaryEphemeris)

        ' There are several mutually exclusive ways to specify the target body.

        ' Use the SofaJpl star catalog, a subset of the Hipparcos catalog complete to mag. 3.
#If False Then
        ' Create a HipparcosCatalog object.
        Dim catalog As New SofaJpl.HipparcosCatalog()
        ' Look up the star name in the catalog we just constructed, and construct a Star object.
        ' The current SofaJpl implementation requires a full match (not substring) to a designation
        ' in the name dictionary. It may be necessary to examine the dictionary (star_names.vot)
        ' with a text editor to find the correct name of a star.
        Dim body As SofaJpl.Body = catalog.GetStar("name Vega", eph)
#End If

        ' Star, from data you supply. The first parameter is the epoch of the data, J2000.0 in
        ' this case.
#If False Then
        Dim body As New SofaJpl.Star(SofaJpl.JulianDate.J2000Base, 293.0899579 / 15.0,
            +69.6611767, 173.77, 597.482, -1738.313, 26.78, eph, "sig dra")
#End If

        ' Solar system body in the JPL ephemeris. The correction for light time iterates until the
        ' solution converges to the specified angle accuracy.
#If True Then
        Dim body As New SofaJpl.JplBody(SofaJpl.JplEphemeris.Body.Venus, eph,
             angleAccuracy)
#End If


        ' Finished setting parameters. Calculate and display results.

        ' Compute the angle, time, and rectangular coordinate resolution consistent with the angle
        ' accuracy (radians) specified by user. These variables control output formatting.
        _angleResolution = 1.0 / SofaJpl.Angle.RadiansToDegrees(angleAccuracy)
        _timeResolution = _angleResolution * _degPerHour

        ' Since CInt rounds to the nearest integer, the expression yields, for example, 2 decimal
        ' places of vector resolution when the angle accuracy is 1/10 to 1/100 radian. 
        _vectorResolution = CInt(-Math.Log10(angleAccuracy) + 0.5)
        If (_vectorResolution < 1) Then
            _vectorResolution = 1      ' Always display at least 1 decimal place.
        End If

        ' Load the leap second table.
        SofaJpl.Utc.LoadTableFromFile(_leapsecTable)

        ' Calculate UT1 and TT. UTC is not calculated (yet) unless it's the input time scale.
        Dim ut1, tt As SofaJpl.JulianDate

        Dim utc As SofaJpl.Utc
#If TIMESCALE = utc Then
        utc = New SofaJpl.Utc(year, month, day, hour, minute, second, julian)
        tt = utc.TerrestrialTime
        ut1 = New SofaJpl.JulianDate(year, month, day, hour, minute, second, julian) + ut1MinusUtc
        deltaT = tt - ut1
#End If
#If TIMESCALE = ut1 Then
        ut1 = New SofaJpl.JulianDate(year, month, day, hour, minute, second, julian)
        tt = ut1 + deltaT
#End If
#If TIMESCALE = tt Then
        tt = New SofaJpl.JulianDate(year, month, day, hour, minute, second, julian)
        ut1 = tt - deltaT
#End If


        ' Display date and time in UTC.

        ' The TimeFields class breaks down a JulianDate into year, month, etc.
        Dim tf1 As SofaJpl.TimeFields

#If TIMESCALE = utc Then
        tf1 = utc.ToTimeFields(_timeResolution, julian)
        Console.WriteLine("{0} UTC", tf1)
#Else
        ' The time input from the user was UT1 or TT, so it's not safe to assume conversion to UTC
        ' is possible. It must be compared to the UTC table boundaries, which are in terms of TAI.
        Dim tai As SofaJpl.JulianDate = tt - SofaJpl.Duration.TTMinusTai
        If tai >= SofaJpl.Utc.DefaultTable.FirstTai And tai <= SofaJpl.Utc.DefaultTable.LastTai Then
            utc = New SofaJpl.Utc(tt)
            tf1 = utc.ToTimeFields(_timeResolution, julian)
            Console.WriteLine("{0} UTC", tf1)
        End If
#End If

        ' Display UT1.
        tf1 = ut1.ToTimeFields(_timeResolution, julian)
        Console.WriteLine("{0} UT1", tf1)

        ' Display TT as date and time.
        tf1 = tt.ToTimeFields(_timeResolution, julian)
        Console.WriteLine("{0} TT", tf1)
        If julian Then
            Console.WriteLine("Dates are Julian calendar.")
        Else
            Console.WriteLine("Dates are Gregorian calendar.")
        End If
        ' Display TT as 2-part Julian date.
        Console.WriteLine("JD {0} TT", tt)

        ' Display delta T.
        Dim sex1 As New SofaJpl.Sexagesimal(deltaT.ToHours(), _timeResolution)
        Console.WriteLine("{0:+fhms} delta T", sex1)

        ' Display polar motion.
        Console.WriteLine(vbCrLf & "polar motion")
        displayPolarMotion(poleX, poleY, "x, y")

        ' Display topocenter geodetic and rectangular coordinates.
        Console.WriteLine()
        Console.WriteLine(obsName)     ' name of observatory
        displayGeodetic(lon, lat, "E lon, N lat")
        Console.WriteLine("{0:f1} meters above ellipsoid", height)
        Dim obsVec As SofaJpl.Vector = obs.ToVector()
        displayXyz(obsVec, "ITRS unit vector", "modulus (km)")

        ' Display deflection of the vertical. The same format as polar motion is appropriate.
        Console.WriteLine(vbCrLf & "deflection of the vertical")
        displayPolarMotion(xi, eta, "xi, eta")

        ' Display atmosphere conditions.
        Console.WriteLine(vbCrLf & "{0:f0} C ({1:f0} F) at observer",
            atm.StationTemperature,
            atm.StationTemperature * 1.8 + 32.0)
        Console.WriteLine("{0,6:f1} mb ({1,5:f2}"" Hg) altimeter setting",
            atm.AltimeterSetting,
            atm.AltimeterSetting / SofaJpl.Atmosphere.MillibarsPerInchHg)
        Console.WriteLine("{0,6:f1} mb ({1,5:f2}"" Hg) station pressure",
            atm.StationPressure,
            atm.StationPressure / SofaJpl.Atmosphere.MillibarsPerInchHg)
        If (atm.HumidityIsRelative) Then
            Console.WriteLine("{0:f0}% relative humidity", atm.Humidity)
        Else
            Console.WriteLine("{0:f0} C ({1:f0} F) dew point",
                atm.Humidity,
                atm.Humidity * 1.8 + 32.0)
        End If

        ' Create rotation matrices for the coordinate transformations.

        ' Old code. To be deleted.
#If False Then
        ' GCRS to celestial intermediate system
        Dim gcrsToCirs = SofaJpl.RMatrix.GcrsToCirs06b(tt)

        ' GCRS to true equator and equinox of date
        Dim gcrsToTrue = SofaJpl.RMatrix.GcrsToTrue06b(tt)

        ' GCRS to ecliptic and true equinox
        Dim gcrsToEcliptic = SofaJpl.RMatrix.GcrsToEclipTrue06b(tt)

        ' GCRS to ITRS, including polar motion
        Dim gcrsToItrs = SofaJpl.RMatrix.GcrsToItrs06b(tt, ut1, poleX, poleY)

        ' ITRS to horizontal (east, north, zenith) system, including deflection of the vertical
        Dim itrstoHor = SofaJpl.RMatrix.ItrsToHor(lat, lon, xi, eta)

        ' GCRS to horizontal, is a product of two of the preceding matrices
        Dim gcrsToHor As SofaJpl.RMatrix = itrstoHor * gcrsToItrs
#End If

        ' ITRS to horizontal (east, north, zenith) system, including deflection of the vertical
        Dim itrstoHor = SofaJpl.RMatrix.ItrsToHor(lat, lon, xi, eta)

        ' polar motion matrix (terrestrial intermediate system to ITRS)
        Dim tirsToItrs = SofaJpl.RMatrix.TirsToItrs(tt, poleX, poleY)

        ' GCRS to true equator & equinox (IAU 2006 precession and 2000B nutation).
        Dim gcrsToMean06 = SofaJpl.RMatrix.Precess06(tt)
        Dim eps06 = SofaJpl.RMatrix.MeanObliq06(tt)     ' mean obliquity
        Dim dPsi00, dEps00 As Double      ' nutation in longitude and obliquity
        SofaJpl.RMatrix.NutationAngles00b(tt, dPsi00, dEps00)
        Dim gcrsToTrue06 = SofaJpl.RMatrix.Nutate(eps06, dPsi00, dEps00) * gcrsToMean06

#If HORIZONS Then
        ' Generate coordinates compatible with JPL Horizons (1976/80 precession/nutation).

        ' In SofaJpl, 1976 precession includes frame bias, which must be removed for
        ' Horizons compatibility.
        Dim gcrsToMean = SofaJpl.RMatrix.Precess76(tt) * SofaJpl.RMatrix.IcrsToJ2000.Transpose()

        ' Calculate the GCRS to ecliptic (mean equinox) rotation matrix.
        Dim eps = SofaJpl.RMatrix.MeanObliq80(tt)     ' mean obl.
        Dim gcrsToEclipticMean = SofaJpl.RMatrix.GcrsToEclip(gcrsToMean, eps)

        ' GCRS vector to the CIP, based on the IAU 2006 precession & 2000 nutation models.
        Dim cipGcrs06 = gcrsToTrue06.Row(3)
        ' Transform it to spherical coords in the 1976/80 ecliptic and mean equinox system.
        Dim cipEclip06 = gcrsToEclipticMean * cipGcrs06
        Dim cipSph = New SofaJpl.Spherical(cipEclip06)

        ' Derive and apply the nutation angles to obtain the GCRS to true equator/equinox matrix.
        Dim dPsi, dEps As Double
        dPsi = SofaJpl.Angle.HalfPi - cipSph.LonEast
        dEps = cipSph.NPD - eps
        Dim gcrsToTrue = SofaJpl.RMatrix.Nutate(eps, dPsi, dEps) * gcrsToMean

        ' Compute the GCRS to terrestrial intermediate matrix
        Dim gast = SofaJpl.Angle.Gast94(ut1)
        Dim gcrsToTirs = SofaJpl.RMatrix.GcrsToTirs(gcrsToTrue, gast)

#Else   ' Generate coordinates compatible with IAU 2006/00 precession/nutation.

        ' Mean obliquity and nutation in obliquity have already been computed.
        Dim eps = eps06
        Dim dEps = dEps00
        ' GCRS to true equator/equinox matrix has already been computed.
        Dim gcrsToTrue = gcrsToTrue06

        ' Get X and Y of the celestial intermediate pole.
        Dim cipX, cipY As Double
        gcrsToTrue.CipXY(cipX, cipY)
        ' Form the GCRS to celestial intermediate matrix.
        Dim gcrsToCirs = SofaJpl.RMatrix.GcrsToCirs(cipX, cipY,
            SofaJpl.RMatrix.S06(tt, cipX, cipY))

        ' Greenwich apparent sidereal time.
        Dim gast = SofaJpl.Angle.Gast06b(ut1, tt)
        ' Earth rotation angle
        Dim era = SofaJpl.Angle.Era00(ut1)

        ' Compute the GCRS to terrestrial intermediate matrix
        Dim gcrsToTirs = SofaJpl.RMatrix.GcrsToTirs(gcrsToCirs, era)
#End If
        ' End code compatible with IAU 2006/00 precession/nutation. The remaining rotation
        ' matrix computations are common to IAU 2006/00 and JPL Horizons modes.

        ' GCRS to ecliptic and true equinox
        Dim gcrsToEcliptic = SofaJpl.RMatrix.GcrsToEclip(gcrsToTrue, eps + dEps)

        ' GCRS to ITRS, including polar motion
        Dim gcrsToItrs = tirsToItrs * gcrsToTirs

        ' GCRS to horizontal, including deflection of the vertical.
        Dim gcrsToHor = itrstoHor * gcrsToItrs


        ' Display barycentric coordinates of the body. 

        Console.WriteLine(vbCrLf & "{0} barycentric position & velocity", body.Name)
        Dim pv1 = body.Barycentric(tt)
        displayRaDec(pv1.Position, "RA, dec (ICRS)")
        displayXyz(pv1.Position, "unit vector (ICRS)", "distance (km)")
        displayXyz(pv1.Velocity, "velocity unit vector", "km/day")

        ' heliocentric coordinates

        Console.WriteLine(vbCrLf & "{0} heliocentric position & velocity", body.Name)
        pv1 = body.Heliocentric(tt)
        displayRaDec(pv1.Position, "RA, dec (ICRS)")
        displayXyz(pv1.Position, "unit vector (ICRS)", "distance (km)")
        displayXyz(pv1.Velocity, "velocity unit vector", "km/day")

        ' geocentric coordinates

        Console.WriteLine(vbCrLf & "{0} geocentric geometric position & velocity", body.Name)
        pv1 = body.GeocentricGeometric(tt)
        displayRaDec(pv1.Position, "RA, dec (ICRS)")
        displayXyz(pv1.Position, "unit vector (ICRS)", "geometric distance (km)")
        displayXyz(pv1.Velocity, "velocity unit vector", "km/day")

        Console.WriteLine(vbCrLf & "{0} geocentric astrometric place", body.Name)
        Dim vec1 = body.GeocentricAstrometric(tt)
        displayRaDec(vec1, "RA, dec (ICRS)")
        displayModulus(vec1, "astrometric distance (km)")

        Console.WriteLine(vbCrLf & "{0} geocentric apparent place", body.Name)
        vec1 = body.GeocentricApparent(tt)
        displayRaDec(vec1, "RA, dec (ICRS)")
        displayRaDec(gcrsToTrue * vec1, "equinox RA, dec")
#If Not HORIZONS Then
        displayRaDec(gcrsToCirs * vec1, "intermediate RA, dec")
#End If
        displayEcliptical(gcrsToEcliptic * vec1, "ecliptic true lon, lat")

        ' Geographic position of the body in the ITRS. If the user provided the parameters, this
        ' includes polar motion.
        Dim sph1 As New SofaJpl.Spherical(gcrsToItrs * vec1)
        displayGeodetic(sph1.LonEast, sph1.Lat, "E lon, N lat (ITRS)")
        displayGeodetic(sph1.LonWest, sph1.Lat, "W lon, N lat (ITRS)")
        displayHourAngle(SofaJpl.Angle.NormPlus(sph1.LonWest + lon), "LHA")

        ' Semidiameter. All bodies have a property to give radius in km. Stars default to zero.
        ' Solar system bodies are automatically initialized to their adopted IAU radii. This
        ' property can be modified by the user.
        Dim sd = SofaJpl.Angle.Semidiameter(vec1.Modulus(), body.Radius)
        displaySemidiameter(sd, "geocentric semidiameter")


        ' Greenwich apparent sidereal time (IAU 2006 precession, 2000B nutation)
        Console.WriteLine()
        displaySiderealTime(SofaJpl.Angle.Gast06b(ut1, tt), "Greenwich apparent sidereal time")

        ' Earth rotation angle.
#If Not HORIZONS Then
        displayHourAngle(era, "Earth rotation angle")
#End If

        ' topocentric coordinates

        ' Get the GCRS position and velocity of the topocenter in a PVVector (position and velocity
        ' vector). The parameter passed to ToGcrsPV() must be the ITRS to GCRS rotation matrix, but
        ' what we calculated earlier does the opposite transformation. Thus it's transposed to
        ' reverse its sense.
        Dim obsPV = obs.ToGcrsPV(gcrsToItrs.Transpose())

        Console.WriteLine(vbCrLf & "{0} topocentric geometric place", body.Name)
        vec1 = body.TopocentricGeometric(tt, obsPV).Position
        displayRaDec(vec1, "RA, dec (ICRS)")
        displayModulus(vec1, "km")

        Console.WriteLine(vbCrLf & "{0} topocentric astrometric place", body.Name)
        vec1 = body.TopocentricAstrometric(tt, obsPV.Position)
        displayRaDec(vec1, "RA, dec (ICRS)")
        displayModulus(vec1, "km")

        Console.WriteLine(vbCrLf & "{0} topocentric apparent place", body.Name)
        vec1 = body.TopocentricApparent(tt, obsPV)
        displayRaDec(vec1, "RA, dec (ICRS)")

        ' Transform apparent place from the GCRS to the horizontal system. If the user supplied
        ' the parameters, this includes polar motion and deflection of the vertical.
        vec1 = gcrsToHor * vec1

        ' Convert to spherical coordinates.
        Dim sphUnref As New SofaJpl.Spherical(vec1)

        ' Use the Atmosphere object created earlier to generate refracted spherical coordinates.
        ' The application of refraction is iterative, so the desired accuracy must be passed
        ' to Refract().
        Dim sphRefr As SofaJpl.Spherical = atm.Refract(sphUnref, angleAccuracy)

        Console.WriteLine(vbCrLf & "{0} azimuth, zenith distance, elevation", body.Name)
        displayAzZd(sphUnref, "az, unrefracted ZD")
        displayAzEl(sphUnref, "az, unrefracted el")
        displayAzEl(sphRefr, "az, refracted el")

        ' Semidiameter. All bodies have a property to give radius in km. Stars default to zero, but
        ' solar system bodies are automatically initialized to their adopted IAU radii.

        sd = SofaJpl.Angle.Semidiameter(vec1.Modulus(), body.Radius)
        displaySemidiameter(sd, "topocentric semidiameter")


        ' If applicable, display phase angle: the separation angle, at the body, between
        ' vectors directed to the Sun and the observer. These are the negatives of vectors to the
        ' body's heliocentric and topocentric positions in the ICRS. If we omit the negations,
        ' both vectors are off by 180 degrees and thus the angle between them is still correct.

        If body.IsSolarSystemBody And body.IsSun = False Then
            Dim earth As New SofaJpl.JplBody(SofaJpl.JplEphemeris.Body.Earth, eph, angleAccuracy)
            Dim sun As New SofaJpl.JplBody(SofaJpl.JplEphemeris.Body.Sun, eph, angleAccuracy)
            Dim earthVec = body.TopocentricGeometric(tt, obsPV).Position
            Dim sunVec = body.Heliocentric(tt).Position
            Dim pa As Double = earthVec.SeparationAngle(sunVec)
            Console.WriteLine("{0:f0}° phase angle (0 = full, 180 = new)",
                SofaJpl.Angle.RadiansToDegrees(pa))
        End If
    End Sub



    ' helper methods to display data with the format and precision selected by user


    ''' <summary>
    ''' Display the modulus of a vector.
    ''' </summary>
    ''' <param name="vec">the vector</param>
    ''' <param name="label">string to display after the modulus</param>
    Sub displayModulus(ByVal vec As SofaJpl.Vector, ByVal label As String)
        Dim sb1 As New System.Text.StringBuilder("{0:e", 8)
        ' The _vectorResolution value is appropriate for the 'f' format. But the 'e' format gives
        ' one more significant digit, so subtract 1.
        sb1.Append(_vectorResolution - 1)
        sb1.Append("} ")
        Console.WriteLine(sb1.ToString() + label, vec.Modulus())
    End Sub


    ''' <summary>
    ''' Display hour angle as a sexagesimal.
    ''' </summary>
    ''' <param name="angle">hour angle (radians) </param>
    ''' <param name="str">string to display after the angle</param>
    Sub displayHourAngle(ByVal angle As Double, ByVal str As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:3c°'""} {1}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1}"
            Case Else
                formatString = "{0:3a°'""} {1}"
        End Select
        Dim sex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(angle), _angleResolution, 360)
        Console.WriteLine(formatString, sex, str)
    End Sub


    ''' <summary>
    ''' Display sidereal time as a sexagesimal.
    ''' </summary>
    ''' <param name="gast">sidereal time (radians)</param>
    ''' <param name="label">string to display after the angle</param>
    ''' <remarks>Unless the DMS format is selected, display degrees, not hours.</remarks>
    Sub displaySiderealTime(ByVal gast As Double, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:3c°'""} {1}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1}"
            Case Else
                formatString = "{0:2ahms} {1}"
        End Select

        Dim sex As SofaJpl.Sexagesimal
        If _format = angleFormat.DMS Then
            sex = New SofaJpl.Sexagesimal(
                SofaJpl.Angle.RadiansToHours(gast), _timeResolution, 24)
        Else
            sex = New SofaJpl.Sexagesimal(
                SofaJpl.Angle.RadiansToDegrees(gast), _angleResolution, 360)
        End If

        Console.WriteLine(formatString, sex, label)
    End Sub


    ''' <summary>
    ''' Display semidiameter as a sexagesimal.
    ''' </summary>
    ''' <param name="sd">semidiameter (radians)</param>
    ''' <param name="label">string to display after the angle</param>
    Sub displaySemidiameter(ByVal sd As Double, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:c°'""} {1}"
            Case angleFormat.DM
                formatString = "{0:e°'""} {1}"
            Case Else
                formatString = "{0:f°'""} {1}"
        End Select
        sd = SofaJpl.Angle.RadiansToDegrees(sd)
        Dim sex As New SofaJpl.Sexagesimal(sd, _angleResolution)
        Console.WriteLine(formatString, sex, label)
    End Sub


    ''' <summary>
    ''' Display polar motion angles as sexagesimals.
    ''' </summary>
    ''' <param name="poleX">pole X (radians)</param>
    ''' <param name="poleY">pole Y (radians)</param>
    ''' <param name="label">string to display after the angles</param>
    Sub displayPolarMotion(ByVal poleX As Double, ByVal poleY As Double, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:+c°'""} {1:+c°'""}  {2}"
            Case angleFormat.DM
                formatString = "{0:+e°'""} {1:+e°'""}  {2}"
            Case Else
                formatString = "{0:+f°'""} {1:+f°'""}  {2}"
        End Select
        Dim xDeg = SofaJpl.Angle.RadiansToDegrees(poleX)
        Dim yDeg = SofaJpl.Angle.RadiansToDegrees(poleY)
        Dim sexX As New SofaJpl.Sexagesimal(xDeg, _angleResolution)
        Dim sexY As New SofaJpl.Sexagesimal(yDeg, _angleResolution)
        Console.WriteLine(formatString, sexX, sexY, label)
    End Sub


    ''' <summary>
    ''' Display a vector as right ascension (hours) and declination.
    ''' </summary>
    ''' <param name="vec">vector to the body with respect to the equatorial system</param>
    ''' <param name="label">string to display after the angles</param>
    ''' <remarks>Unless the DMS format is selected, display RA in degrees, not hours.</remarks>
    Sub displayRaDec(ByVal vec As SofaJpl.Vector, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:3c°'""} {1:+2c°'""} {2}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1:+2b°'""} {2}"
            Case Else
                formatString = "{0:2ahms} {1:+2a°'""} {2}"
        End Select

        Dim sph As New SofaJpl.Spherical(vec)

        Dim raSex As SofaJpl.Sexagesimal
        If _format = angleFormat.DMS Then
            raSex = New SofaJpl.Sexagesimal(
                SofaJpl.Angle.RadiansToHours(sph.RA), _timeResolution, 24)
        Else
            raSex = New SofaJpl.Sexagesimal(
                SofaJpl.Angle.RadiansToDegrees(sph.RA), _angleResolution, 360)
        End If

        Dim decSex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.Dec), _angleResolution)

        Console.WriteLine(formatString, raSex, decSex, label)
    End Sub


    ''' <summary>
    ''' Display a vector as azimuth and elevation.
    ''' </summary>
    ''' <param name="vec">vector to the body with respect to the horizontal system</param>
    ''' <param name="label">string to display after the angles</param>
    Sub displayAzEl(ByVal vec As SofaJpl.Vector, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:3c°'""} {1:+2c°'""} {2}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1:+2b°'""} {2}"
            Case Else
                formatString = "{0:3a°'""} {1:+2a°'""} {2}"
        End Select

        Dim sph As New SofaJpl.Spherical(vec)

        Dim azSex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.Az), _angleResolution, 360)

        Dim elSex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.El), _angleResolution)

        Console.WriteLine(formatString, azSex, elSex, label)
    End Sub


    ''' <summary>
    ''' Display a vector as ecliptic longitude and latitude.
    ''' </summary>
    ''' <param name="vec">vector to the body with respect to the ecliptic system</param>
    ''' <param name="label">string to display after the angles</param>
    Sub displayEcliptical(ByVal vec As SofaJpl.Vector, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:3c°'""} {1:+2c°'""} {2}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1:+2b°'""} {2}"
            Case Else
                formatString = "{0:3a°'""} {1:+2a°'""} {2}"
        End Select

        Dim sph As New SofaJpl.Spherical(vec)

        Dim lonSex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.LonEast), _angleResolution, 360)

        Dim latSex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.Lat), _angleResolution)

        Console.WriteLine(formatString, lonSex, latSex, label)
    End Sub


    ''' <summary>
    ''' Display geodetic longitude and latitude.
    ''' </summary>
    ''' <param name="lon">east longitude (radians)</param>
    ''' <param name="lat">north latitude (radians)</param>
    ''' <param name="label">string to display after the angles</param>
    Sub displayGeodetic(ByVal lon As Double, ByVal lat As Double, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:3c°'""} {1:+2c°'""} {2}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1:+2b°'""} {2}"
            Case Else
                formatString = "{0:3a°'""} {1:+2a°'""} {2}"
        End Select

        Dim lonSex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(lon), _angleResolution, 360)

        Dim latSex As New SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(lat), _angleResolution)

        Console.WriteLine(formatString, lonSex, latSex, label)
    End Sub


    ''' <summary>
    ''' Display a vector as xyz components of a unit vector and modulus.
    ''' </summary>
    ''' <param name="v">vector</param>
    ''' <param name="xyzLabel">label for the unit vector</param>
    ''' <param name="modulusLabel">label for the modulus</param>
    Sub displayXyz(ByVal v As SofaJpl.Vector, ByVal xyzLabel As String,
    ByVal modulusLabel As String)
        Dim sb1 As New System.Text.StringBuilder("{0:f", 23)
        sb1.Append(_vectorResolution)
        sb1.Append("} {1:f")
        sb1.Append(_vectorResolution)
        sb1.Append("} {2:f")
        sb1.Append(_vectorResolution)
        sb1.Append("} ")

        Dim uv As SofaJpl.Vector = v.Unit()
        Console.WriteLine(sb1.ToString() & xyzLabel, uv.X, uv.Y, uv.Z)

        displayModulus(v, modulusLabel)
    End Sub


    ''' <summary>
    ''' Given a Spherical object, display azimuth and elevation as sexagesimals.
    ''' </summary>
    ''' <param name="sph">spherical coordinates in the horizontal system</param>
    ''' <param name="label">string to display after the angles</param>
    Sub displayAzEl(ByVal sph As SofaJpl.Spherical, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                formatString = "{0:3c°'""} {1: 2c°'""} {2}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1: 2b°'""} {2}"
            Case Else
                formatString = "{0:3a°'""} {1: 2a°'""} {2}"
        End Select
        Dim az As Double = SofaJpl.Angle.RadiansToDegrees(sph.Az)
        Dim el As Double = SofaJpl.Angle.RadiansToDegrees(sph.El)
        Dim azSex As New SofaJpl.Sexagesimal(az, _angleResolution, 360)
        Dim elSex As New SofaJpl.Sexagesimal(el, _angleResolution)
        Console.WriteLine(formatString, azSex, elSex, label)
    End Sub


    ''' <summary>
    ''' Given a Spherical object, display azimuth and zenith distance as sexagesimals.
    ''' </summary>
    ''' <param name="sph">spherical coordinates</param>
    ''' <param name="label">string to display after the angles</param>
    Sub displayAzZd(ByVal sph As SofaJpl.Spherical, ByVal label As String)
        Dim formatString As String
        Select Case _format
            Case angleFormat.D
                ' Zenith distance needs a different format than elevation. There's never a negative
                ' sign, and there can be up to 3 digits before the decimal point.
                formatString = "{0:3c°'""} {1:3c°'""} {2}"
            Case angleFormat.DM
                formatString = "{0:3b°'""} {1:3b°'""} {2}"
            Case Else
                formatString = "{0:3a°'""} {1:3a°'""} {2}"
        End Select
        Dim az As Double = SofaJpl.Angle.RadiansToDegrees(sph.Az)
        Dim zd As Double = SofaJpl.Angle.RadiansToDegrees(sph.ZenithDistance)
        Dim azSex As New SofaJpl.Sexagesimal(az, _angleResolution, 360)
        Dim zdSex As New SofaJpl.Sexagesimal(zd, _angleResolution)
        Console.WriteLine(formatString, azSex, zdSex, label)
    End Sub
End Module

[back]

(last modified 2019-05-27)