//-------------------------------------------------------------------------------------
// SID Monitor - Utility for Sudden Ionospheric Disturbances Monitoring Stations
// Copyright (C) 2006 - Lionel Loudet
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
//-------------------------------------------------------------------------------------


using System;
using System.Collections.Generic;
using System.Text;

namespace SID_monitor
{
    class SunriseSunset
    {
        #region Constants
        private const Double degs = 180.0 / Math.PI;  // radians to degrees conversion
        private const Double rads = Math.PI / 180.0;  // degrees to radians conversion

        //h = 0 degrees: Center of Sun's disk touches a mathematical horizon
        //h = -0.25 degrees: Sun's upper limb touches a mathematical horizon
        //h = -0.583 degrees: Center of Sun's disk touches the horizon; atmospheric refraction accounted for
        //h = -0.833 degrees: Sun's supper limb touches the horizon; atmospheric refraction accounted for
        //h = -6 degrees: Civil twilight (one can no longer read outside without artificial illumination)
        //h = -12 degrees: Nautical twilight (navigation using a sea horizon no longer possible)
        //h = -15 degrees: Amateur astronomical twilight (the sky is dark enough for most astronomical observations)
        //h = -18 degrees: Astronomical twilight (the sky is completely dark)
        private const Double RiseSetAngle = -0.833 * rads;
        private const Double CivilTwilightAngle = -6.0 * rads;
        private const Double NautiTwilightAngle = -12.0 * rads;
        private const Double AstroTwilightAngle = -18.0 * rads;

        // The mean sidereal day is 23 hours, 56 minutes, 4.09054 seconds long
        private const Double MSD = 1.00273790935; //one mean solar day = 1.00273790935 mean sidereal days.
        private const Double CalculationPrecision = 1.0 / 3600.0; // 1 sec precision of the computed times

        private const double AirRefr = 34.0 / 60.0; // athmospheric refraction in degrees 
        #endregion

        private Double civil_twam, civil_twpm, nauti_twam, nauti_twpm, astro_twam, astro_twpm; // civil, nautical, astronomical twilight morning and evening times
        private Double altmax;       // sun altitude at noon time
        private Double altmin;       // sun altitude at midnight time
        private Double daylen;       // daylength
        private Double riseTime;     // sunrise Time
        private Double setTime;      // sunset Time
        private Double noonTime;     // noon time


        /// <summary>
        /// program calculating the sunrise and sunset for
        /// the current date and a fixed location(latitude,longitude)
        /// 
        /// Adapted from Jarmo Lammi (http://personal.inet.fi/cool/jjlammi/) and Paul Schlyter (http://stjarnhimlen.se/)
        /// </summary>
        /// <param name="Aujourdhui">The date for which the times are calculated</param>
        public SunriseSunset(DateTime Aujourdhui)
        {

            // calculates latitude and longitude ======================================================================================
            Double latitude = SID_monitor.Properties.Settings.Default.LatDeg + SID_monitor.Properties.Settings.Default.LatMins / 60 + SID_monitor.Properties.Settings.Default.LatSecs / 3600;
            if (SID_monitor.Properties.Settings.Default.LatHemi == "S")
            {
                latitude = -latitude;
            }
            latitude = latitude * rads; // convert to radians


            Double longitude = SID_monitor.Properties.Settings.Default.LongDeg + SID_monitor.Properties.Settings.Default.LongMins / 60 + SID_monitor.Properties.Settings.Default.LongSecs / 3600;
            if (SID_monitor.Properties.Settings.Default.LongHemi == "W")
            {
                longitude = -longitude;
            }
            longitude = longitude * rads; // convert to radians


            //Orbital elements of the Sun:
            //N = 0.0                           // longitude of the ascending node
            //i = 0.0                           // inclination to the ecliptic (plane of the Earth's orbit)
            //w = 282.9404 + 4.70935E-5 * d     // argument of perihelion
            //a = 1.000000  (AU)                // semi-major axis, or mean distance from Sun
            //e = 0.016709 - 1.151E-9 * d       // eccentricity (0=circle, 0-1=ellipse, 1=parabola)
            //M = 356.0470 + 0.9856002585 * d   // mean anomaly (0 at perihelion; increases uniformly with time)
            //ecl = 23.4393 - 3.563E-7 * d      // Obliquity of the ecliptic

            //Related orbital elements are: 
            //w1 = N + w   = longitude of perihelion
            //L  = M + w1  = mean longitude
            //q  = a*(1-e) = perihelion distance
            //Q  = a*(1+e) = aphelion distance
            //P  = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)
            //T  = Epoch_of_M - (M(deg)/360_deg) / P  = time of perihelion
            //v  = true anomaly (angle between position and perihelion)
            //E  = eccentric anomaly

            #region Algorithm without iteration (less precise)
            //Double RA = Double.NaN;
            //Double Dec = Double.NaN;
            //Double GMST0 = Double.NaN;

            //DateTime currentTime = Aujourdhui.Date.AddHours(12); // we start by computing the Sun's RA and Dec, and sidereal time at noon local time

            //CalculateSunPosition(currentTime, out RA, out Dec, out GMST0);
            ////UT_Sun_in_south = ( RA - GMST0 - long ) / 15.0 in hours
            //Double UT_Sun_in_south = FNrange(RA - GMST0 - longitude); // in radians

            //// Noon Time
            //noonTime = UT_Sun_in_south / (Math.PI / 12.0);

            ////Now we're going to compute the Sun's Local Hour Angle (LHA) at rise/set (or at twilight, if we've decided to compute the time of twilight).
            //Double cosLHA0 = getCosLHA(latitude, Dec, RiseSetAngle);
            //Double cosLHA6 = getCosLHA(latitude, Dec, CivilTwilightAngle);
            //Double cosLHA12 = getCosLHA(latitude, Dec, NautiTwilightAngle);
            //Double cosLHA18 = getCosLHA(latitude, Dec, AstroTwilightAngle);


            //// Sunrise and sunset times
            //if ((cosLHA0 >= -1.0) && (cosLHA0 <= 1.0))  //If cos(LHA) is between +1.0 and -1.0, then we take the arccos to find LHA. Convert from degrees to hours by dividing by 15.0
            //{
            //    //Now, if we add LHA to UT_Sun_in_south, we get the time of sunset. If we subtract LHA from UT_Sun_in_south, we get the time of sunrise.
            //    riseTime = (UT_Sun_in_south - Math.Acos(cosLHA0)) / (Math.PI / 12.0);
            //    setTime = (UT_Sun_in_south + Math.Acos(cosLHA0)) / (Math.PI / 12.0);
            //    daylen = setTime - riseTime;
            //}
            //else if (cosLHA0 < -1.0)  //If cos(LHA) is less than -1.0, then the Sun is always below our altitude limit. One example is midwinter in the arctics, when the Sun never gets above the horizon.
            //{
            //    riseTime = Double.NaN;
            //    setTime = Double.NaN;
            //    daylen = 0.0;
            //}
            //else //If cos(LHA) is greater than +1.0, then the Sun is always above our altitude limit. If we were computing rise/set times, the Sun is then aboute the horizon continuously; we have Midnight Sun. Or, if we computed a twilight, then the sky never gets dark (a good example is Stockholm, Sweden, at midsummer midnight: the Sun then only reaches about 7 degrees below the horizon: there will be civil twilight, but never nautical or astronomical twilight, at midsummer in Stockholm).
            //{
            //    riseTime = Double.NaN;
            //    setTime = Double.NaN;
            //    daylen = 24.0;
            //}


            //// sun altitude at noon time
            //altmax = 90.0 + Dec * degs - latitude * degs;
            //if (latitude < Dec)
            //{
            //    altmax = 180.0 - altmax;
            //}

            //// civil twilight times
            //if ((cosLHA6 >= -1.0) && (cosLHA6 <= 1.0))
            //{
            //    civil_twam = (UT_Sun_in_south - Math.Acos(cosLHA6)) / (Math.PI / 12.0);      // morning civil twilight begin
            //    civil_twpm = (UT_Sun_in_south + Math.Acos(cosLHA6)) / (Math.PI / 12.0);      // evening civil twilight end
            //}
            //else
            //{
            //    civil_twam = Double.NaN;
            //    civil_twpm = Double.NaN;
            //}

            //// nautical twilight times
            //if ((cosLHA12 >= -1.0) && (cosLHA12 <= 1.0))
            //{
            //    nauti_twam = (UT_Sun_in_south - Math.Acos(cosLHA12)) / (Math.PI / 12.0);     // morning nautical twilight begin
            //    nauti_twpm = (UT_Sun_in_south + Math.Acos(cosLHA12)) / (Math.PI / 12.0);     // evening nautical twilight end
            //}
            //else
            //{
            //    nauti_twam = Double.NaN;
            //    nauti_twpm = Double.NaN;
            //}

            //// astronomical twilight times
            //if ((cosLHA18 >= -1.0) && (cosLHA18 <= 1.0))
            //{
            //    astro_twam = (UT_Sun_in_south - Math.Acos(cosLHA18)) / (Math.PI / 12.0);     // morning astronomical twilight begin
            //    astro_twpm = (UT_Sun_in_south + Math.Acos(cosLHA18)) / (Math.PI / 12.0);     // evening astronomical twilight end
            //}
            //else
            //{
            //    astro_twam = Double.NaN;
            //    astro_twpm = Double.NaN;
            //}
            #endregion

            // Noon Time and altitude calculations ====================================================================================
            DateTime currentTime = Aujourdhui.Date.AddHours(12); // we start by computing the Sun's RA and Dec, and sidereal time at noon local time
            Double RA = Double.NaN;
            Double Dec = Double.NaN;
            Double GMST0 = Double.NaN;

            CalculateSunPosition(currentTime, out RA, out Dec, out GMST0);
            //UT_Sun_in_south = ( RA - GMST0 - long ) / 15.0 in hours
            Double UT_Sun_in_south = FNrange(RA - GMST0 - longitude); // in radians

            // Noon Time
            noonTime = UT_Sun_in_south / (Math.PI / 12.0);

            // sun altitudes
            altmax = Math.Asin(Math.Sin(latitude) * Math.Sin(Dec) + Math.Cos(latitude) * Math.Cos(Dec)) * degs; // sun altitude at noon time
            altmin = Math.Asin(Math.Sin(latitude) * Math.Sin(Dec) - Math.Cos(latitude) * Math.Cos(Dec)) * degs; // sun altitude at midnight time
            if ((altmax > 0.0) && (altmax < 30.0)) { altmax += AirRefr; }          // Air refraction included at small altitudes
            if ((altmin > 0.0) && (altmin < 30.0)) { altmin += AirRefr; }          // Air refraction included at small altitudes

            // Sunrise and sunset times ===============================================================================================
            riseTime = CalculateSunTime(Aujourdhui, longitude, latitude, RiseSetAngle, true);
            setTime = CalculateSunTime(Aujourdhui, longitude, latitude, RiseSetAngle, false);

            if (!(Double.IsNaN(riseTime) || Double.IsInfinity(riseTime)) &&
                !(Double.IsNaN(setTime) || Double.IsInfinity(setTime)))
            {
                daylen = setTime - riseTime;
            }
            else
            {
                if (Double.IsPositiveInfinity(setTime))
                {
                    daylen = 0.0;
                }
                if (Double.IsNegativeInfinity(setTime))
                {
                    daylen = 24.0;
                }
                riseTime = Double.NaN;
                setTime = Double.NaN;
            }


            // civil twilight times ===================================================================================================
            civil_twam = CalculateSunTime(Aujourdhui, longitude, latitude, CivilTwilightAngle, true);
            civil_twpm = CalculateSunTime(Aujourdhui, longitude, latitude, CivilTwilightAngle, false);

            if (Double.IsInfinity(civil_twam) || Double.IsInfinity(civil_twpm))
            {
                civil_twam = Double.NaN;
                civil_twpm = Double.NaN;
            }

            // nautical twilight times ================================================================================================
            nauti_twam = CalculateSunTime(Aujourdhui, longitude, latitude, NautiTwilightAngle, true);
            nauti_twpm = CalculateSunTime(Aujourdhui, longitude, latitude, NautiTwilightAngle, false);

            if (Double.IsInfinity(nauti_twam) || Double.IsInfinity(nauti_twpm))
            {
                nauti_twam = Double.NaN;
                nauti_twpm = Double.NaN;
            }

            // astronomical twilight times ============================================================================================
            astro_twam = CalculateSunTime(Aujourdhui, longitude, latitude, AstroTwilightAngle, true);
            astro_twpm = CalculateSunTime(Aujourdhui, longitude, latitude, AstroTwilightAngle, false);

            if (Double.IsInfinity(astro_twam) || Double.IsInfinity(astro_twpm))
            {
                astro_twam = Double.NaN;
                astro_twpm = Double.NaN;
            }
        }

        /// <summary>
        /// Returns the time that correspond to the Sun's position at the given Angle altitude
        /// </summary>
        /// <param name="Aujourdhui">The date at which the calculation is performed</param>
        /// <param name="longitude">observer's longitude</param>
        /// <param name="latitude">observer's latitude</param>
        /// <param name="Angle">The altitude of the sun for which the time is calculated</param>
        /// <param name="rise">true for a sunrise calculation, false for a sunset calculation</param>
        /// <returns>the time in the 0-24h range. -inf if no sunrise. +inf if no sunset.</returns>
        private Double CalculateSunTime(DateTime Aujourdhui, Double longitude, Double latitude, Double Angle, bool rise)
        {

            int i = 0; // iteration counter

            Double currentTime = 12.0; // we start by computing the Sun's RA / Dec and sidereal time at noon local time
            DateTime oldTime;

            Double UT_Sun_in_south;
            Double cosLHA;

            Double RA;  // Suns's RA
            Double Dec; // Suns's declination
            Double GMST0; // sidereal time

            do // we iterate till we get the desired precision
            {
                if (i++ > 10) // we limit the number of iterations
                {
                    Program.MainForm.outputTextBoxDockablePanel.AddOutputTextWarningMessage("Bad precision of the sunrise time calculation.\n");
                    break;
                }

                oldTime = Aujourdhui.Date.AddHours(currentTime);
                CalculateSunPosition(oldTime, out RA, out Dec, out GMST0);

                //UT_Sun_in_south = ( RA - GMST0 - long ) / 15.0 in hours
                UT_Sun_in_south = FNrange(RA - GMST0 - longitude); // in radians

                //Now we're going to compute the Sun's Local Hour Angle (LHA) at rise/set (or at twilight, if we've decided to compute the time of twilight).
                cosLHA = getCosLHA(latitude, Dec, Angle);

                if ((cosLHA >= -1.0) && (cosLHA <= 1.0))  //If cos(LHA) is between +1.0 and -1.0, then we take the arccos to find LHA. Convert from degrees to hours by dividing by 15.0
                {
                    //Now, if we add LHA to UT_Sun_in_south, we get the time of sunset. If we subtract LHA from UT_Sun_in_south, we get the time of sunrise.
                    if (rise) // rise time asked
                    {
                        currentTime = (UT_Sun_in_south - Math.Acos(cosLHA)) / (Math.PI / 12.0) / MSD;
                    }
                    else
                    {
                        currentTime = (UT_Sun_in_south + Math.Acos(cosLHA)) / (Math.PI / 12.0) / MSD;
                    }
                    if (currentTime >= 24.0) { currentTime -= 24.0; }
                    if (currentTime < 0.0) { currentTime += 24.0; }
                }
                else if (cosLHA < -1.0)  //If cos(LHA) is less than -1.0, then the Sun is always below our altitude limit. One example is midwinter in the arctics, when the Sun never gets above the horizon.
                // or the opposite?
                {
                    currentTime = Double.NegativeInfinity;
                }
                else //If cos(LHA) is greater than +1.0, then the Sun is always above our altitude limit. If we were computing rise/set times, the Sun is then aboute the horizon continuously; we have Midnight Sun. Or, if we computed a twilight, then the sky never gets dark (a good example is Stockholm, Sweden, at midsummer midnight: the Sun then only reaches about 7 degrees below the horizon: there will be civil twilight, but never nautical or astronomical twilight, at midsummer in Stockholm).
                // or the opposite?
                {
                    currentTime = Double.PositiveInfinity;
                }
            }
            while ((!(Double.IsNaN(currentTime) || Double.IsInfinity(currentTime)))
                   && (Math.Abs(currentTime - oldTime.TimeOfDay.TotalHours) > CalculationPrecision));

            return currentTime;
        }

        /// <summary>
        /// Returns the cosine of the Sun's Local Hour Angle
        /// </summary>
        /// <param name="lat">observer's latitude in radians</param>
        /// <param name="Decl">Sun's declination in radians</param>
        /// <param name="h">given sun specific altitude in radians</param>
        private Double getCosLHA(Double lat, Double Decl, double h)
        {
            // This is the angle the Earth must turn from sunrise to noon, or from noon to sunset:
            //            sin(h) - sin(lat)*sin(Decl)
            //cos(LHA) = -----------------------------
            //            cos(lat) * cos(Decl)

            return (Math.Sin(h) - Math.Sin(lat) * Math.Sin(Decl)) / (Math.Cos(lat) * Math.Cos(Decl));
        }

        /// <summary>
        /// Returns sun coordinates (RA and Dec) and the sidereal time at Greenwich at 00:00UT
        /// </summary>
        /// <param name="date">the date at which the calculation is done</param>
        /// <param name="RA">returns the sun RA</param>
        /// <param name="Dec">returns the sun Dec</param>
        /// <param name="GMST0">returns the sun sidereal time</param>
        private void CalculateSunPosition(DateTime date, out Double RA, out Double Dec, out Double GMST0)
        {

            // days since J2000
            Double d = date.Year * 367 - 7 * (date.Year + (date.Month + 9) / 12) / 4 + 275 * date.Month / 9 + date.Day - 730531.5 + date.TimeOfDay.Hours / 24.0;

            //First, compute the eccentric anomaly E from the mean anomaly M and from the eccentricity e (E and M in degrees):
            //E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )
            //or (if E and M are expressed in radians):
            //E = M + e * sin(M) * ( 1.0 + e * cos(M) )


            Double M = FNrange((356.0470 + 0.9856002585 * d) * rads);  // Mean Anomaly
            Double ecc = (0.016709 - 1.151E-9 * d);           // Eccentricity
            Double E = FNrange(M + ecc * Math.Sin(M) * (1.0 + ecc * Math.Cos(M))); // Eccentricity Anomaly

            //Then compute the Sun's distance r and its true anomaly v from:
            //xv = r * cos(v) = cos(E) - e
            //yv = r * sin(v) = sqrt(1.0 - e*e) * sin(E)
            //v = atan2( yv, xv )
            //r = sqrt( xv*xv + yv*yv )

            Double xv = Math.Cos(E) - ecc;
            Double yv = Math.Sqrt(1.0 - ecc * ecc) * Math.Sin(E);
            Double v = Math.Atan2(yv, xv);
            Double r = Math.Sqrt(xv * xv + yv * yv);


            //Now, compute the Sun's true longitude:
            //lonsun = v + w
            Double w = (282.9404 + 4.70935E-5 * d) * rads;
            Double lonsun = FNrange(v + w);

            //Convert lonsun,r to ecliptic rectangular geocentric coordinates xs,ys:
            //xs = r * cos(lonsun)
            //ys = r * sin(lonsun)
            //(since the Sun always is in the ecliptic plane, zs is of course zero). xs,ys is the Sun's position in a coordinate system in the plane of the ecliptic.

            Double xs = r * Math.Cos(lonsun);
            Double ys = r * Math.Sin(lonsun);

            //To convert this to equatorial, rectangular, geocentric coordinates, compute:
            //xe = xs
            //ye = ys * cos(ecl)
            //ze = ys * sin(ecl)

            Double ecl = (23.4393 - 3.563E-7 * d) * rads;
            Double xe = xs;
            Double ye = ys * Math.Cos(ecl);
            Double ze = ys * Math.Sin(ecl);


            //Finally, compute the Sun's Right Ascension (RA) and Declination (Dec):
            //RA  = atan2( ye, xe )
            //Dec = atan2( ze, sqrt(xe*xe+ye*ye) )

            RA = Math.Atan2(ye, xe);
            Dec = Math.Atan2(ze, Math.Sqrt(xe * xe + ye * ye));

            //We start by computing the sidereal time at Greenwich at 00:00 Universal Time, let's call this quantity GMST0:
            //GMST0 = L + 180
            //L is the Sun's mean longitude, which we compute as: 
            //L = M + w

            GMST0 = FNrange(M + w + Math.PI);

        }

        #region Properties
        public String SunsetTime
        {
            get
            {
                if (Double.IsNaN(this.setTime))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.setTime*60.0)/60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String SunriseTime
        {
            get
            {
                if (Double.IsNaN(this.riseTime))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.riseTime*60.0)/60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String MorningCivilTwilight
        {
            get
            {
                if (Double.IsNaN(this.civil_twam))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.civil_twam * 60.0) / 60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String MorningNauticalTwilight
        {
            get
            {
                if (Double.IsNaN(this.nauti_twam))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.nauti_twam * 60.0) / 60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String MorningAstronomicalTwilight
        {
            get
            {
                if (Double.IsNaN(this.astro_twam))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.astro_twam * 60.0) / 60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String EveningCivilTwilight
        {
            get
            {
                if (Double.IsNaN(this.civil_twpm))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.civil_twpm * 60.0) / 60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String EveningNauticalTwilight
        {
            get
            {
                if (Double.IsNaN(this.nauti_twpm))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.nauti_twpm * 60.0) / 60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String EveningAstronomicalTwilight
        {
            get
            {
                if (Double.IsNaN(this.astro_twpm))
                {
                    return "-----";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.astro_twpm * 60.0) / 60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5) + " UT";
                    }
                    else
                    {
                        return outString + " UT";
                    }
                }
            }
        }

        public String SunMaxAlt
        {
            get
            {
                return this.altmax.ToString("F2") + "°";
            }
        }

        public String SunMinAlt
        {
            get
            {
                return this.altmin.ToString("F2") + "°";
            }
        }

        public String NoonTime
        {
            get
            {
                TimeSpan val = new TimeSpan();
                val = TimeSpan.FromHours(Math.Round(this.noonTime * 60.0) / 60.0);
                String outString = val.ToString();
                if (outString.Length > 5)
                {
                    return outString.Remove(5) + " UT";
                }
                else
                {
                    return outString + " UT";
                }
            }
        }

        public String DayLength
        {
            get
            {
                if (this.daylen == 24.0)
                {
                    return "24:00";
                }
                else
                {
                    TimeSpan val = new TimeSpan();
                    val = TimeSpan.FromHours(Math.Round(this.daylen*60.0)/60.0);
                    String outString = val.ToString();
                    if (outString.Length > 5)
                    {
                        return outString.Remove(5);
                    }
                    else
                    {
                        return outString;
                    }
                }
            }
        }

        /// <summary>
        ///  returns
        /// 0 for day
        /// 1 for evening civil twilight
        /// 2 for evening nautical twilight
        /// 3 for evening astronomical twilight
        /// 4 for night
        /// 5 for morning astronomical twilight
        /// 6 for morning nautical twilight
        /// 7 for morning civil twilight
        /// </summary>
        /// <param name="current_time"></param>
        /// <returns></returns>
        public int DayNightStatus(double current_time)
        {
            int day_night_status;

            if (current_time <= this.astro_twam)
            {
                day_night_status = 4; // night
            }
            else if (current_time <= this.nauti_twam)
            {
                day_night_status = 5; // astronomical morning twilight
            }
            else if (current_time <= this.civil_twam)
            {
                day_night_status = 6; // nautical morning twilight
            }
            else if (current_time <= this.riseTime)
            {
                day_night_status = 7; // civil morning twilight
            }
            else if (current_time <= this.setTime)
            {
                day_night_status = 0; // day
            }
            else if (current_time <= this.civil_twpm)
            {
                day_night_status = 1; // civil evening twilight
            }
            else if (current_time <= this.nauti_twpm)
            {
                day_night_status = 2; // nautical evening twilight
            }
            else if (current_time <= this.astro_twpm)
            {
                day_night_status = 3; // astronomical evening twilight
            }
            else
            {
                day_night_status = 4; // night
            }

            return day_night_status;

        }
        #endregion


        /// <summary>
        /// the function below returns an angle in the range 0 to 2*pi
        /// </summary>
        private Double FNrange(Double x)
        {
            Double b = 0.5 * x / Math.PI;
            Double a = 2.0 * Math.PI * (b - (Int64)(b));
            if (a < 0) a = 2.0 * Math.PI + a;
            return a;
        }

    }
}
