// SofaJplDemo.cs
// 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 HirosePS. SofaJpl
// namespace.

// revision history
#if false
2018-09-07 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.
#endif


// Select the input time scale. Exactly one must be #defined, and the others commented out.
#define UTC
// #define UT1
// #define TT

// If #define is un-commented, convert a JPL ASCII ephemeris to binary.
// After the binary is created, this #define should be commented out.
// #define MAKE_EPHEMERIS

// If #define is un-commented, use the delta T you provide. Otherwise use the SofaJpl delta T model.
// Has no effect if time scale is UTC.
// #define MANUAL_DELTA_T

// If un-commented, 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.
// #define HORIZONS


using System;

// Everything in SofaJpl is in the HirosePS.SofaJpl namespace.
using SofaJpl = HirosePS.SofaJpl;


static class SofaJplDemo {

    // The leap second table. One is included in the SofaJpl distribution .zip file.
    static string _leapsecTable = @"C:\Users\Stan\Documents\astro\SofaJpl_2_1_2\leapSecTable.txt";

#if MAKE_EPHEMERIS
    // JPL ASCII ephemeris and header. Must be downloaded from JPL. Not needed after they are
    // converted to a binary ephemeris.
    static string _asciiEphemeris = @"C:\Users\Stan\Documents\astro\jpl\ascp2000.422";
    static string _ephemerisHeader = @"C:\Users\Stan\Documents\astro\jpl\header.422";
#endif

    // The binary JPL ephemeris to create and use.
    static string _binaryEphemeris = @"binp2000_2099.422";

    // mathematical constants
    const double _degPerHour = 15.0;

    // an enumeration of the three possible angle output formats
    enum angleFormat { D, DM, DMS };

    // output format desired by the user
    static angleFormat _format;

    // These control the output resolution of angle and time. The values are computed automatically
    // from the angle accuracy set by the user.
    static double _angleResolution;      // resolution units per degree
    static double _timeResolution;       // 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.
    static int _vectorResolution;


    static void Main(string[] args) {

        // Set the display format and angle accuracy.
        _format = angleFormat.DMS;
        double angleAccuracy = SofaJpl.Angle.DmsToRad(0, 0, 0.1);

        // Set the epoch.
        int year = 2018;
        int month = 6;
        int day = 6;
        int hour = 4;
        int minute = 0;
        double second = 0.0;
        bool julian = false;    // true if date is in Julian calendar


        SofaJpl.Duration deltaT;
#if UTC
        SofaJpl.Duration ut1MinusUtc = SofaJpl.Duration.FromSeconds(0.07522);
#else
#if MANUAL_DELTA_T
        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));
#endif
#endif

        // Observer position.
        string obsName = "Kitt Peak Observatory";    // name of topocenter (optional)
        double lon = SofaJpl.Angle.DmsToRad(-111, 36, 00.0);    // east longitude
        double lat = SofaJpl.Angle.DmsToRad(31, 57, 48.0);      // north latitude
        float height = 2100f;    // meters above ellipsoid

        // This instance of the Topocenter class will supply the observer position and velocity
        // with respect to the ICRS.
        SofaJpl.Topocenter obs = 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.
        float altimiterSetting = 29.9f * SofaJpl.Atmosphere.MillibarsPerInchHg;     // millibars
        float degC = 12f;       // Celsius
        float dewPointC = 0f;   // Celsius.
        SofaJpl.Atmosphere atm = new SofaJpl.Atmosphere(height, degC, altimiterSetting,
            dewPointC, false);

        // polar motion (radians)
        double poleX = SofaJpl.Angle.DmsToRad(0, 0, 0.0);
        double poleY = SofaJpl.Angle.DmsToRad(0, 0, 0.0);

        // deflection of the vertical at the observer (radians)
        double xi = SofaJpl.Angle.DmsToRad(0, 0, 0);
        double eta = SofaJpl.Angle.DmsToRad(0, 0, 0);


        // Create a binary ephemeris from JPL ASCII files.
#if MAKE_EPHEMERIS
        Console.WriteLine("\nCREATING BINARY JPL EPHEMERIS.\n");
        SofaJpl.JplEphemeris.AsciiToBinary(_ephemerisHeader, _asciiEphemeris, _binaryEphemeris);
#endif

        // 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.
        SofaJpl.JplEphemeris eph = new SofaJpl.JplEphemeris(_binaryEphemeris);

        // There are several ways to specify the target body.

        // Use the SofaJpl star catalog, a subset of the Hipparcos catalog complete to mag. 3.
#if false
        // Create a HipparcosCatalog object.
        SofaJpl.HipparcosCatalog catalog = 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.
         SofaJpl.Body body = catalog.GetStar("name Vega", eph);
#endif

        // Star, from data you supply. The first parameter is the epoch of the data, J2000.0 in
        // this case.
#if false
        SofaJpl.Body body = new SofaJpl.Star(SofaJpl.JulianDate.J2000Base, 293.0899579 / 15.0,
            +69.6611767, 173.77, 597.482, -1738.313, 26.78, eph, "sig Dra");
#endif

        // Solar system body in the JPL ephemeris. The correction for light time iterates until the
        // solution converges to the specified angle accuracy.
#if true
        SofaJpl.Body body = new SofaJpl.JplBody(SofaJpl.JplEphemeris.Body.Venus, eph,
             angleAccuracy);
#endif


        // 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;
        _vectorResolution = (int)(-Math.Log10(angleAccuracy) + 1.0);
        if (_vectorResolution < 1)
            _vectorResolution = 1;      // Always display at least 1 decimal place.

        // 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.
        SofaJpl.JulianDate ut1, tt;

        SofaJpl.Utc utc;
#if UTC
        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;
#endif
#if UT1
        ut1 = new SofaJpl.JulianDate(year, month, day, hour, minute, second, julian);
        tt = ut1 + deltaT;
#endif
#if TT
        tt = new SofaJpl.JulianDate(year, month, day, hour, minute, second, julian);
        ut1 = tt - deltaT;
#endif


        // Display date and time in UTC.

        // The TimeFields class breaks down a JulianDate into year, month, etc.
        SofaJpl.TimeFields tf1;

#if UTC
        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.
        SofaJpl.JulianDate tai = tt - SofaJpl.Duration.TTMinusTai;
        if (tai >= SofaJpl.Utc.DefaultTable.FirstTai && tai <= SofaJpl.Utc.DefaultTable.LastTai) {
            utc = new SofaJpl.Utc(tt);
            tf1 = utc.ToTimeFields(_timeResolution, julian);
            Console.WriteLine("{0} UTC", tf1);
        }
#endif

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

        // Display TT as date and time, and also as 2-part Julian date.
        tf1 = tt.ToTimeFields(_timeResolution, julian);
        Console.WriteLine("{0} TT", tf1);
        Console.WriteLine("Dates are {0} calendar.", julian ? "Julian" : "Gregorian");
        Console.WriteLine("JD {0} TT", tt);

        // Display delta T.
        SofaJpl.Sexagesimal sex1 = new SofaJpl.Sexagesimal(deltaT.ToHours(), _timeResolution);
        Console.WriteLine("{0:+fhms} delta T", sex1);

        // Display polar motion.
        Console.WriteLine("\npolar 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);
        SofaJpl.Vector obsVec = obs.ToVector();
        displayXyz(obsVec, "ITRS unit vector", "modulus (km)");

        // Display deflection of the vertical. The same format as polar motion is appropriate.
        Console.WriteLine("\ndeflection of the vertical");
        displayPolarMotion(xi, eta, "xi, eta");

        // Display atmosphere conditions.
        Console.WriteLine("\n{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)
            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);


        // Create rotation matrices for the coordinate transformations.

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

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

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

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

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

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

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

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

        // Compute the GCRS to terrestrial intermediate matrix
        double gast = SofaJpl.Angle.Gast94(ut1);
        SofaJpl.RMatrix 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.
        double eps = eps06;
        double dEps = dEps00;
        // GCRS to true equator/equinox matrix has already been computed.
        SofaJpl.RMatrix gcrsToTrue = gcrsToTrue06;

        // Get X and Y of the celestial intermediate pole.
        double cipX, cipY;
        gcrsToTrue.CipXY(out cipX, out cipY);
        // Form the GCRS to celestial intermediate matrix.
        SofaJpl.RMatrix gcrsToCirs = SofaJpl.RMatrix.GcrsToCirs(cipX, cipY,
            SofaJpl.RMatrix.S06(tt, cipX, cipY));

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

        // Compute the GCRS to terrestrial intermediate matrix
        SofaJpl.RMatrix gcrsToTirs = SofaJpl.RMatrix.GcrsToTirs(gcrsToCirs, era);
#endif
        // 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
        SofaJpl.RMatrix gcrsToEcliptic = SofaJpl.RMatrix.GcrsToEclip(gcrsToTrue, eps + dEps);

        // GCRS to ITRS, including polar motion
        SofaJpl.RMatrix gcrsToItrs = tirsToItrs * gcrsToTirs;

        // GCRS to horizontal, including deflection of the vertical.
        SofaJpl.RMatrix gcrsToHor = itrstoHor * gcrsToItrs;


        // Display barycentric coordinates of the body. 

        Console.WriteLine("\n{0} barycentric position & velocity", body.Name);
        SofaJpl.PVVector 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("\n{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("\n{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("\n{0} geocentric astrometric place", body.Name);
        SofaJpl.Vector vec1 = body.GeocentricAstrometric(tt);
        displayRaDec(vec1, "RA, dec (ICRS)");
        displayModulus(vec1, "astrometric distance (km)");

        Console.WriteLine("\n{0} geocentric apparent place", body.Name);
        vec1 = body.GeocentricApparent(tt);
        displayRaDec(vec1, "RA, dec (ICRS)");
        displayRaDec(gcrsToTrue * vec1, "equinox RA, dec");
#if ! HORIZONS
        displayRaDec(gcrsToCirs * vec1, "intermediate RA, dec");
#endif
        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.
        SofaJpl.Spherical sph1 = 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.
        double sd = SofaJpl.Angle.Semidiameter(vec1.Modulus(), body.Radius);
        displaySemidiameter(sd, "geocentric semidiameter");


        // Greenwich apparent sidereal time
        Console.WriteLine();
        displaySiderealTime(gast, "Greenwich apparent sidereal time");

#if ! HORIZONS
        // Earth rotation angle.
        displayHourAngle(era, "Earth rotation angle");
#endif

        // 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.
        SofaJpl.PVVector obsPV = obs.ToGcrsPV(gcrsToItrs.Transpose());

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

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

        Console.WriteLine("\n{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.
        SofaJpl.Spherical sphUnref = 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().
        SofaJpl.Spherical sphRefr = atm.Refract(sphUnref, angleAccuracy);

        Console.WriteLine("\n{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 && body.IsSun == false) {
            SofaJpl.JplBody earth = new SofaJpl.JplBody(
                SofaJpl.JplEphemeris.Body.Earth, eph, angleAccuracy);
            SofaJpl.JplBody sun = new SofaJpl.JplBody(
                SofaJpl.JplEphemeris.Body.Sun, eph, angleAccuracy);
            SofaJpl.Vector earthVec = body.TopocentricGeometric(tt, obsPV).Position;
            SofaJpl.Vector sunVec = body.Heliocentric(tt).Position;
            double pa = earthVec.SeparationAngle(sunVec);
            Console.WriteLine("{0:f0}° phase angle (0 = full, 180 = new)",
                SofaJpl.Angle.RadiansToDegrees(pa));
        }
    }   // end Main()



    // 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>
    static void displayModulus(SofaJpl.Vector vec, string label) {
        System.Text.StringBuilder sb1 = 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());
    }


    /// <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>
    static void displayHourAngle(double angle, string str) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:3c°'\"} {1}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1}";
            break;
        default:
            formatString = "{0:3a°'\"} {1}";
            break;
        }
        SofaJpl.Sexagesimal sex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(angle), _angleResolution, 360);
        Console.WriteLine(formatString, sex, str);
    }


    /// <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>Display degrees (not hours) unless format is DMS.</remarks>
    static void displaySiderealTime(double gast, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:3c°'\"} {1}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1}";
            break;
        default:
            formatString = "{0:2ahms} {1}";
            break;
        }
        SofaJpl.Sexagesimal sex;
        if (_format == angleFormat.DMS)
            sex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToHours(gast),
            _timeResolution, 24);
        else
            sex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToDegrees(gast),
            _angleResolution, 360);
        Console.WriteLine(formatString, sex, label);
    }


    /// <summary>
    /// Display semidiameter as a sexagesimal.
    /// </summary>
    /// <param name="sd">semidiameter (radians)</param>
    /// <param name="label">string to display after the angle</param>
    static void displaySemidiameter(double sd, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:c°'\"} {1}";
            break;
        case angleFormat.DM:
            formatString = "{0:e°'\"} {1}";
            break;
        default:
            formatString = "{0:f°'\"} {1}";
            break;
        }
        sd = SofaJpl.Angle.RadiansToDegrees(sd);
        SofaJpl.Sexagesimal sex = new SofaJpl.Sexagesimal(sd, _angleResolution);
        Console.WriteLine(formatString, sex, label);
    }


    /// <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>
    static void displayPolarMotion(double poleX, double poleY, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:+c°'\"} {1:+c°'\"}  {2}";
            break;
        case angleFormat.DM:
            formatString = "{0:+e°'\"} {1:+e°'\"}  {2}";
            break;
        default:
            formatString = "{0:+f°'\"} {1:+f°'\"}  {2}";
            break;
        }
        double xDeg = SofaJpl.Angle.RadiansToDegrees(poleX);
        double yDeg = SofaJpl.Angle.RadiansToDegrees(poleY);
        SofaJpl.Sexagesimal sexX = new SofaJpl.Sexagesimal(xDeg, _angleResolution);
        SofaJpl.Sexagesimal sexY = new SofaJpl.Sexagesimal(yDeg, _angleResolution);
        Console.WriteLine(formatString, sexX, sexY, label);
    }


    /// <summary>
    /// Display a vector as right ascension 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 format is DMS, display RA in degrees not hours.</remarks>
    static void displayRaDec(SofaJpl.Vector vec, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
            break;
        default:
            formatString = "{0:2ahms} {1:+2a°'\"} {2}";
            break;
        }

        SofaJpl.Spherical sph = new SofaJpl.Spherical(vec);

        SofaJpl.Sexagesimal raSex;
        if (_format == angleFormat.DMS)
            raSex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToHours(sph.RA),
                _timeResolution, 24);
        else
            raSex = new SofaJpl.Sexagesimal(SofaJpl.Angle.RadiansToDegrees(sph.RA),
                _angleResolution, 360);

        SofaJpl.Sexagesimal decSex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.Dec), _angleResolution);

        Console.WriteLine(formatString, raSex, decSex, label);
    }


    /// <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>
    static void displayAzEl(SofaJpl.Vector vec, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
            break;
        default:
            formatString = "{0:3a°'\"} {1:+2a°'\"} {2}";
            break;
        }

        SofaJpl.Spherical sph = new SofaJpl.Spherical(vec);

        SofaJpl.Sexagesimal azSex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.Az), _angleResolution, 360);

        SofaJpl.Sexagesimal elSex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.El), _angleResolution);

        Console.WriteLine(formatString, azSex, elSex, label);
    }


    /// <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>
    static void displayEcliptical(SofaJpl.Vector vec, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
            break;
        default:
            formatString = "{0:3a°'\"} {1:+2a°'\"} {2}";
            break;
        }

        SofaJpl.Spherical sph = new SofaJpl.Spherical(vec);

        SofaJpl.Sexagesimal lonSex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.LonEast), _angleResolution, 360);

        SofaJpl.Sexagesimal latSex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(sph.Lat), _angleResolution);

        Console.WriteLine(formatString, lonSex, latSex, label);
    }


    /// <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>
    static void displayGeodetic(double lon, double lat, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:3c°'\"} {1:+2c°'\"} {2}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1:+2b°'\"} {2}";
            break;
        default:
            formatString = "{0:3a°'\"} {1:+2a°'\"} {2}";
            break;
        }

        SofaJpl.Sexagesimal lonSex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(lon), _angleResolution, 360);

        SofaJpl.Sexagesimal latSex = new SofaJpl.Sexagesimal(
            SofaJpl.Angle.RadiansToDegrees(lat), _angleResolution);

        Console.WriteLine(formatString, lonSex, latSex, label);
    }


    /// <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>
    static void displayXyz(SofaJpl.Vector v, string xyzLabel, string modulusLabel) {
        System.Text.StringBuilder sb1 = 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("} ");

        SofaJpl.Vector uv = v.Unit();
        Console.WriteLine(sb1.ToString() + xyzLabel, uv.X, uv.Y, uv.Z);

        displayModulus(v, modulusLabel);
    }


    /// <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>
    static void displayAzEl(SofaJpl.Spherical sph, string label) {
        string formatString;
        switch (_format) {
        case angleFormat.D:
            formatString = "{0:3c°'\"} {1: 2c°'\"} {2}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1: 2b°'\"} {2}";
            break;
        default:
            formatString = "{0:3a°'\"} {1: 2a°'\"} {2}";
            break;
        }
        double az = SofaJpl.Angle.RadiansToDegrees(sph.Az);
        double el = SofaJpl.Angle.RadiansToDegrees(sph.El);
        SofaJpl.Sexagesimal azSex = new SofaJpl.Sexagesimal(az, _angleResolution, 360);
        SofaJpl.Sexagesimal elSex = new SofaJpl.Sexagesimal(el, _angleResolution);
        Console.WriteLine(formatString, azSex, elSex, label);
    }


    /// <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>
    static void displayAzZd(SofaJpl.Spherical sph, string label) {
        string formatString;
        switch (_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}";
            break;
        case angleFormat.DM:
            formatString = "{0:3b°'\"} {1:3b°'\"} {2}";
            break;
        default:
            formatString = "{0:3a°'\"} {1:3a°'\"} {2}";
            break;
        }
        double az = SofaJpl.Angle.RadiansToDegrees(sph.Az);
        double zd = SofaJpl.Angle.RadiansToDegrees(sph.ZenithDistance);
        SofaJpl.Sexagesimal azSex = new SofaJpl.Sexagesimal(az, _angleResolution, 360);
        SofaJpl.Sexagesimal zdSex = new SofaJpl.Sexagesimal(zd, _angleResolution);
        Console.WriteLine(formatString, azSex, zdSex, label);
    }
}   // end class

[back]

(last modified 2018-09-09)