CSharp: SunTimeCalculator

 

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

namespace GeovinDu.Ticket.Common
{
    /// <summary>
    /// 
    /// </summary>
    public static class DateTimeJavaScriptExt
    {

        /// <summary>
        /// Based on a JavaScript library SunCalc for calculating sun/moon position and light phases.
        ///https://github.com/mourner/suncalc
        /// </summary>
        /// <param name="dt"></param>
        /// <returns></returns>
        public static double ValueOf(this DateTime dt) // JavaScript Date.valueOf()
        {
            dt = dt.Kind == DateTimeKind.Local ? dt.ToUniversalTime() : dt;
            return (dt - new DateTime(1970, 1, 1, 0, 0, 0, DateTimeKind.Utc)).TotalMilliseconds;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="dt"></param>
        /// <param name="ms"></param>
        /// <returns></returns>
        public static DateTime FromJScriptValue(this DateTime dt, double ms)
        {
            return new DateTime(1970, 1, 1, 0, 0, 0, DateTimeKind.Utc).AddMilliseconds(ms);
        }
    }
    /// <summary>
    /// 计算日出日落时间，月升月落时间
    /// 
    /// </summary>
    public class SunTimeCalculator
    {


        #region 辅助函数
        /// <summary>
        /// 历元2000.0，即以2000年第一天开端为计日起始（天文学以第一天为0日而非1日）。
        /// 它与UT（就是世界时，格林尼治平均太阳时）1999年末重合。 
        /// </summary>
        /// <param name="y"></param>
        /// <param name="m"></param>
        /// <param name="d"></param>
        /// <returns></returns>
        private static long Days_since_2000_Jan_0(int y, int m, int d)
        {
            return (367L * (y) - ((7 * ((y) + (((m) + 9) / 12))) / 4) + ((275 * (m)) / 9) + (d) - 730530L);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Revolution(double x)
        {
            return (x - 360.0 * Math.Floor(x * Inv360));
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Rev180(double x)
        {
            return (x - 360.0 * Math.Floor(x * Inv360 + 0.5));
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        private static double GMST0(double d)
        {
            double sidtim0;
            sidtim0 = Revolution((180.0 + 356.0470 + 282.9404) +
                (0.9856002585 + 4.70935E-5) * d);
            return sidtim0;
        }

        /// <summary>
        /// 
        /// </summary>
        private static double Inv360 = 1.0 / 360.0;
        #endregion

        #region 度与弧度转换系数，为球面三角计算作准备
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Sind(double x)
        {
            return Math.Sin(x * Degrad);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Cosd(double x)
        {
            return Math.Cos(x * Degrad);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Tand(double x)
        {
            return Math.Tan(x * Degrad);

        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Atand(double x)
        {
            return Radge * Math.Atan(x);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Asind(double x)
        {
            return Radge * Math.Asin(x);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Acosd(double x)
        {
            return Radge * Math.Acos(x);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="y"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double Atan2d(double y, double x)
        {
            return Radge * Math.Atan2(y, x);

        }
        /// <summary>
        /// 
        /// </summary>
        private static double Radge = 180.0 / Math.PI;
        /// <summary>
        /// 
        /// </summary>
        private static double Degrad = Math.PI / 180.0;

        #endregion

        #region 与日出日落时间相关计算
        /// <summary>
        /// 
        /// </summary>
        /// <param name="year"></param>
        /// <param name="month"></param>
        /// <param name="day"></param>
        /// <param name="lon"></param>
        /// <param name="lat"></param>
        /// <param name="altit"></param>
        /// <param name="upper_limb"></param>
        /// <returns></returns>
        private static double DayLen(int year, int month, int day, double lon, double lat,
            double altit, int upper_limb)
        {
            double d,  /* Days since 2000 Jan 0.0 (negative before) */
                obl_ecl,    /* Obliquity (inclination) of Earth's axis */
                //黄赤交角，在2000.0历元下国际规定为23度26分21.448秒，但有很小的时间演化。

                sr,         /* Solar distance, astronomical units */
                slon,       /* True solar longitude */
                sin_sdecl,  /* Sine of Sun's declination */
                //太阳赤纬的正弦值。
                cos_sdecl,  /* Cosine of Sun's declination */
                sradius,    /* Sun's apparent radius */
                t;          /* Diurnal arc */

            /* Compute d of 12h local mean solar time */
            d = Days_since_2000_Jan_0(year, month, day) + 0.5 - lon / 360.0;

            /* Compute obliquity of ecliptic (inclination of Earth's axis) */
            obl_ecl = 23.4393 - 3.563E-7 * d;
            //这个黄赤交角时变公式来历复杂，很大程度是经验性的，不必追究。

            /* Compute Sun's position */
            slon = 0.0;
            sr = 0.0;
            Sunpos(d, ref slon, ref sr);

            /* Compute sine and cosine of Sun's declination */
            sin_sdecl = Sind(obl_ecl) * Sind(slon);
            cos_sdecl = Math.Sqrt(1.0 - sin_sdecl * sin_sdecl);
            //用球面三角学公式计算太阳赤纬。

            /* Compute the Sun's apparent radius, degrees */
            sradius = 0.2666 / sr;
            //视半径，同前。

            /* Do correction to upper limb, if necessary */
            if (upper_limb != 0)
                altit -= sradius;

            /* Compute the diurnal arc that the Sun traverses to reach */
            /* the specified altitide altit: */
            //根据设定的地平高度判据计算周日弧长。
            double cost;
            cost = (Sind(altit) - Sind(lat) * sin_sdecl) /
                (Cosd(lat) * cos_sdecl);
            if (cost >= 1.0)
                t = 0.0;                      /* Sun always below altit */
            //极夜。
            else if (cost <= -1.0)
                t = 24.0;                     /* Sun always above altit */
            //极昼。
            else t = (2.0 / 15.0) * Acosd(cost); /* The diurnal arc, hours */
            //周日弧换算成小时计。
            return t;

        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="d"></param>
        /// <param name="lon"></param>
        /// <param name="r"></param>
        private static void Sunpos(double d, ref double lon, ref double r)
        {
            double M,//太阳的平均近点角，从太阳观察到的地球（=从地球看到太阳的）距近日点（近地点）的角度。
                w, //近日点的平均黄道经度。
                e, //地球椭圆公转轨道离心率。
                E, //太阳的偏近点角。计算公式见下面。

                x, y,
                v;  //真近点角，太阳在任意时刻的真实近点角。


            M = Revolution(356.0470 + 0.9856002585 * d);//自变量的组成：2000.0时刻太阳黄经为356.0470度,此后每天约推进一度（360度/365天
            w = 282.9404 + 4.70935E-5 * d;//近日点的平均黄经。

            e = 0.016709 - 1.151E-9 * d;//地球公转椭圆轨道离心率的时间演化。以上公式和黄赤交角公式一样，不必深究。

            E = M + e * Radge * Sind(M) * (1.0 + e * Cosd(M));
            x = Cosd(E) - e;
            y = Math.Sqrt(1.0 - e * e) * Sind(E);
            r = Math.Sqrt(x * x + y * y);
            v = Atan2d(y, x);
            lon = v + w;
            if (lon >= 360.0)
                lon -= 360.0;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="d"></param>
        /// <param name="RA"></param>
        /// <param name="dec"></param>
        /// <param name="r"></param>
        private static void Sun_RA_dec(double d, ref double RA, ref double dec, ref double r)
        {
            double lon, obl_ecl, x, y, z;
            lon = 0.0;

            Sunpos(d, ref lon, ref r);
            //计算太阳的黄道坐标。

            x = r * Cosd(lon);
            y = r * Sind(lon);
            //计算太阳的直角坐标。

            obl_ecl = 23.4393 - 3.563E-7 * d;
            //黄赤交角，同前。

            z = y * Sind(obl_ecl);
            y = y * Cosd(obl_ecl);
            //把太阳的黄道坐标转换成赤道坐标（暂改用直角坐标）。

            RA = Atan2d(y, x);
            dec = Atan2d(z, Math.Sqrt(x * x + y * y));
            //最后转成赤道坐标。显然太阳的位置是由黄道坐标方便地直接确定的，但必须转换到赤
            //道坐标里才能结合地球的自转确定我们需要的白昼长度。

        }
        /// <summary>
        /// 日出没时刻计算
        /// </summary>
        /// <param name="year">年</param>
        /// <param name="month">月</param>
        /// <param name="day">日</param>
        /// <param name="lon">经度</param>
        /// <param name="lat">纬度</param>
        /// <param name="altit"></param>
        /// <param name="upper_limb"></param>
        /// <param name="trise">日出时刻</param>
        /// <param name="tset">日没时刻</param>
        /// <returns>太阳有出没现象，返回0 极昼，返回+1 极夜，返回-1</returns>
        private static int SunRiset(int year, int month, int day, double lon, double lat,
            double altit, int upper_limb, ref double trise, ref double tset)
        {
            double d,  /* Days since 2000 Jan 0.0 (negative before) */
                //以历元2000.0起算的日数。

                sr,         /* Solar distance, astronomical units */
                //太阳距离，以天文单位计算（约1.5亿公里）。      

                sRA,        /* Sun's Right Ascension */
                //同前，太阳赤经。

                sdec,       /* Sun's declination */
                //太阳赤纬。

                sradius,    /* Sun's apparent radius */
                //太阳视半径，约16分（受日地距离、大气折射等诸多影响）

                t,          /* Diurnal arc */
                //周日弧，太阳一天在天上走过的弧长。

                tsouth,     /* Time when Sun is at south */
                sidtime;    /* Local sidereal time */
            //当地恒星时，即地球的真实自转周期。比平均太阳日（日常时间）长3分56秒。      

            int rc = 0; /* Return cde from function - usually 0 */

            /* Compute d of 12h local mean solar time */
            d = Days_since_2000_Jan_0(year, month, day) + 0.5 - lon / 360.0;
            //计算观测地当日中午时刻对应2000.0起算的日数。

            /* Compute local sideral time of this moment */
            sidtime = Revolution(GMST0(d) + 180.0 + lon);
            //计算同时刻的当地恒星时（以角度为单位）。以格林尼治为基准，用经度差校正。

            /* Compute Sun's RA + Decl at this moment */
            sRA = 0.0;
            sdec = 0.0;
            sr = 0.0;
            Sun_RA_dec(d, ref sRA, ref sdec, ref sr);
            //计算同时刻太阳赤经赤纬。

            /* Compute time when Sun is at south - in hours UT */
            tsouth = 12.0 - Rev180(sidtime - sRA) / 15.0;
            //计算太阳日的正午时刻，以世界时（格林尼治平太阳时）的小时计。

            /* Compute the Sun's apparent radius, degrees */
            sradius = 0.2666 / sr;
            //太阳视半径。0.2666是一天文单位处的太阳视半径（角度）。

            /* Do correction to upper limb, if necessary */
            if (upper_limb != 0)
                altit -= sradius;
            //如果要用上边缘，就要扣除一个视半径。

            /* Compute the diurnal arc that the Sun traverses to reach */
            //计算周日弧。直接利用球面三角公式。如果碰到极昼极夜问题，同前处理。
            /* the specified altitide altit: */

            double cost;
            cost = (Sind(altit) - Sind(lat) * Sind(sdec)) /
                (Cosd(lat) * Cosd(sdec));
            if (cost >= 1.0)
            {
                rc = -1;
                t = 0.0;
            }
            else
            {
                if (cost <= -1.0)
                {
                    rc = +1;
                    t = 12.0;      /* Sun always above altit */
                }
                else
                    t = Acosd(cost) / 15.0;   /* The diurnal arc, hours */
            }


            /* Store rise and set times - in hours UT */
            trise = tsouth - t;
            tset = tsouth + t;

            return rc;
        }
        #endregion

        /// <summary>
        /// 计算日出日没时间
        /// </summary>
        /// <param name="date"></param>
        /// <param name="longitude">经度</param>
        /// <param name="latitude">纬度</param>
        /// <returns></returns>
        public static SunTimeResult GetSunTime(DateTime date, double longitude, double latitude)
        {
            double start = 0;
            double end = 0;
            SunRiset(date.Year, date.Month, date.Day, longitude, latitude, -35.0 / 60.0, 1, ref start, ref end);
            DateTime sunrise = ToLocalTime(date, start);
            DateTime sunset = ToLocalTime(date, end);
            return new SunTimeResult(sunrise, sunset);
        }

        #region moon 月亮计算


        private const double dayMs = 86400000;
        private const double J1970 = 2440588;
        private const double J2000 = 2451545;
        private const double PI = Math.PI;
        private const double rad = Math.PI / 180.0;
        private const double e = rad * 23.4397; // obliquity of the Earth


        /// <summary>
        /// 
        /// </summary>
        public class SunTime
        {
            public double Angle { get; set; }
            public string MorningName { get; set; }
            public string EveningName { get; set; }
        }
        /// <summary>
        /// 
        /// </summary>
        public class SunTimeRiseSet : SunTime
        {
            public DateTime RiseTime { get; set; }
            public DateTime SetTime { get; set; }
        }

        /// <summary>
        /// sun times configuration (angle, morning name, evening name)
        /// </summary>

        public static List<SunTime> SunTimes = new List<SunTime>(new SunTime[]
        {
            new SunTime { Angle = -0.833, MorningName = "sunrise", EveningName = "sunset" },
            new SunTime { Angle = -0.3, MorningName = "sunriseEnd", EveningName = "sunsetStart" },
            new SunTime { Angle = -6, MorningName = "dawn", EveningName = "dusk" },
            new SunTime { Angle = -12, MorningName = "nauticalDawn", EveningName = "nauticalDusk" },
            new SunTime { Angle = -18, MorningName = "nightEnd", EveningName = "night" },
            new SunTime { Angle = 6, MorningName = "goldenHourEnd", EveningName = "goldenHour" }
        });

        /// <summary>
        /// adds a custom time to the times config
        /// </summary>
        /// <param name="sunTime"></param>

        public static void AddTime(SunTime sunTime)
        {
            SunTimes.Add(sunTime);
        }
        /// <summary>
        /// 
        /// </summary>
        public class RaDec
        {
            public double ra = 0;
            public double dec = 0;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="dt"></param>
        /// <returns></returns>
        public static double ToJulianDate(DateTime dt)
        {
            dt = dt.Kind == DateTimeKind.Local ? dt.ToUniversalTime() : dt;
            return (dt - new DateTime(1970, 1, 1, 0, 0, 0, DateTimeKind.Utc)).TotalMilliseconds / dayMs - 0.5 + J1970;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="jd"></param>
        /// <returns></returns>
        public static DateTime FromJulianDate(double jd)
        {
            return double.IsNaN(jd)
                ? DateTime.MinValue
                : new DateTime(1970, 1, 1, 0, 0, 0, DateTimeKind.Utc).AddMilliseconds((jd + 0.5 - J1970) * dayMs);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="dt"></param>
        /// <returns></returns>
        public static double JulianDays(DateTime dt)
        {
            dt = dt.Kind == DateTimeKind.Local ? dt.ToUniversalTime() : dt;
            return ToJulianDate(dt) - J2000;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="l"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        public static double RightAscension(double l, double b)
        {
            return Math.Atan2(Math.Sin(l) * Math.Cos(e) - Math.Tan(b) * Math.Sin(e), Math.Cos(l));
        }

        public static double Declination(double l, double b)
        {
            return Math.Asin(Math.Sin(b) * Math.Cos(e) + Math.Cos(b) * Math.Sin(e) * Math.Sin(l));
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="H"></param>
        /// <param name="phi"></param>
        /// <param name="dec"></param>
        /// <returns></returns>
        public static double Azimuth(double H, double phi, double dec)
        {
            return Math.Atan2(Math.Sin(H), Math.Cos(H) * Math.Sin(phi) - Math.Tan(dec) * Math.Cos(phi));
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="H"></param>
        /// <param name="phi"></param>
        /// <param name="dec"></param>
        /// <returns></returns>
        public static double Altitude(double H, double phi, double dec)
        {
            return Math.Asin(Math.Sin(phi) * Math.Sin(dec) + Math.Cos(phi) * Math.Cos(dec) * Math.Cos(H));
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="d"></param>
        /// <param name="lw"></param>
        /// <returns></returns>
        public static double SiderealTime(double d, double lw)
        {
            return rad * (280.16 + 360.9856235 * d) - lw;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="h"></param>
        /// <returns></returns>
        public static double AstroRefraction(double h)
        {
            if (h < 0) // the following formula works for positive altitudes only.
            {
                h = 0; // if h = -0.08901179 a div/0 would occur.
            }

            // formula 16.4 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
            // 1.02 / tan(h + 10.26 / (h + 5.10)) h in degrees, result in arc minutes -> converted to rad:
            return 0.0002967 / Math.Tan(h + 0.00312536 / (h + 0.08901179));
        }

        /// <summary>
        /// general sun calculations
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>

        public static double SolarMeanAnomaly(double d)
        {
            return rad * (357.5291 + 0.98560028 * d);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="dt"></param>
        /// <param name="h"></param>
        /// <returns></returns>
        public static DateTime HoursLater(DateTime dt, double h)
        {
            return new DateTime(1970, 1, 1, 0, 0, 0, DateTimeKind.Utc).AddMilliseconds(dt.ValueOf() + h * dayMs / 24); //ValueOf
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="M"></param>
        /// <returns></returns>
        public static double EclipticLongitude(double M)
        {
            double C = rad * (1.9148 * Math.Sin(M) + 0.02 * Math.Sin(2 * M) +
                              0.0003 * Math.Sin(3 * M)); // equation of center
            double P = rad * 102.9372; // perihelion of the Earth
            return M + C + P + PI;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        public static RaDec SunCoords(double d)
        {
            double M = SolarMeanAnomaly(d);
            double L = EclipticLongitude(M);
            return new RaDec { dec = Declination(L, 0), ra = RightAscension(L, 0) };
        }
        /// <summary>
        /// 
        /// </summary>
        public class MoonRaDecDist
        {
            public double ra = 0;
            public double dec = 0;
            public double dist = 0;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        public static MoonRaDecDist MoonCoords(double d) // geocentric ecliptic coordinates of the moon
        {
            double L = rad * (218.316 + 13.176396 * d); // ecliptic longitude
            double M = rad * (134.963 + 13.064993 * d); // mean anomaly
            double F = rad * (93.272 + 13.229350 * d); // mean distance

            double l = L + rad * 6.289 * Math.Sin(M); // longitude
            double b = rad * 5.128 * Math.Sin(F); // latitude
            double dt = 385001 - 20905 * Math.Cos(M); // distance to the moon in km

            return new MoonRaDecDist { ra = RightAscension(l, b), dec = Declination(l, b), dist = dt };
        }
        /// <summary>
        /// 
        /// </summary>
        public class MoonAzAltDistPa
        {
            public double azimuth = 0;
            public double altitude = 0;
            public double distance = 0;
            public double parallacticAngle = 0;
        }
        /// <summary>
        /// 
        /// </summary>
        public class MoonFracPhaseAngle
        {
            public double fraction = 0;
            public double phase = 0;
            public double angle = 0;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="dt"></param>
        /// <param name="lat"></param>
        /// <param name="lng"></param>
        /// <returns></returns>
        public static MoonAzAltDistPa GetMoonPosition(DateTime dt, double lat, double lng)
        {
            double lw = rad * -lng;
            double phi = rad * lat;
            double d = JulianDays(dt);

            MoonRaDecDist c = MoonCoords(d);
            double H = SiderealTime(d, lw) - c.ra;
            double h = Altitude(H, phi, c.dec);
            // formula 14.1 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
            double pa = Math.Atan2(Math.Sin(H), Math.Tan(phi) * Math.Cos(c.dec) - Math.Sin(c.dec) * Math.Cos(H));

            h += AstroRefraction(h); // altitude correction for refraction
            return new MoonAzAltDistPa { azimuth = Azimuth(H, phi, c.dec), altitude = h, distance = c.dist, parallacticAngle = pa };
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="dt"></param>
        /// <returns></returns>
        public static MoonFracPhaseAngle GetMoonIllumination(DateTime dt)
        {
            double d = JulianDays(dt);
            RaDec s = SunCoords(d);
            MoonRaDecDist m = MoonCoords(d);
            double sdist = 149598000; // distance from Earth to Sun in km
            double phi = Math.Acos(Math.Sin(s.dec) * Math.Sin(m.dec) +
                                   Math.Cos(s.dec) * Math.Cos(m.dec) * Math.Cos(s.ra - m.ra));
            double inc = Math.Atan2(sdist * Math.Sin(phi), m.dist - sdist * Math.Cos(phi));
            double angle = Math.Atan2(Math.Cos(s.dec) * Math.Sin(s.ra - m.ra), Math.Sin(s.dec) * Math.Cos(m.dec) -
                Math.Cos(s.dec) * Math.Sin(m.dec) * Math.Cos(s.ra - m.ra));
            return new MoonFracPhaseAngle
            {
                fraction = (1 + Math.Cos(inc)) / 2,
                phase = 0.5 + 0.5 * inc * (angle < 0 ? -1 : 1) / PI,
                angle = angle
            };
        }

        /// <summary>
        /// DateTime.Max = always up, DateTime.Min = always down
        /// </summary>
        /// <param name="dt"></param>
        /// <param name="lat">纬度</param>
        /// <param name="lng">经度</param>
        /// <param name="risem"></param>
        /// <param name="setm"></param>
        /// <param name="alwaysUp"></param>
        /// <param name="alwaysDown"></param>
        public static void MoonRiset(DateTime dt, double lat, double lng, out DateTime risem, out DateTime setm,
            out bool? alwaysUp, out bool? alwaysDown)
        {
            dt = new DateTime(dt.Year, dt.Month, dt.Day, 0, 0, 0, DateTimeKind.Utc);
            DateTime t = dt;

            double hc = 0.133 * rad;
            double h0 = GetMoonPosition(t, lat, lng).altitude - hc;
            double h1, h2, rise = 0, set = 0, a, b, xe, ye = 0, d, x1, x2, dx;
            int roots;

            for (double i = 1.0; i <= 24.0; i += 2.0)
            {
                h1 = GetMoonPosition(HoursLater(t, i), lat, lng).altitude - hc;
                h2 = GetMoonPosition(HoursLater(t, i + 1), lat, lng).altitude - hc;

                a = (h0 + h2) / 2 - h1;
                b = (h2 - h0) / 2;
                xe = -b / (2 * a);
                ye = (a * xe + b) * xe + h1;
                d = b * b - 4 * a * h1;
                roots = 0;

                if (d >= 0)
                {
                    dx = Math.Sqrt(d) / (Math.Abs(a) * 2);
                    x1 = xe - dx;
                    x2 = xe + dx;
                    if (Math.Abs(x1) <= 1)
                    {
                        roots++;
                    }

                    if (Math.Abs(x2) <= 1)
                    {
                        roots++;
                    }

                    if (x1 < -1)
                    {
                        x1 = x2;
                    }

                    if (roots == 1)
                    {
                        if (h0 < 0)
                        {
                            rise = i + x1;
                        }
                        else
                        {
                            set = i + x1;
                        }
                    }
                    else if (roots == 2)
                    {
                        rise = i + (ye < 0 ? x2 : x1);
                        set = i + (ye < 0 ? x1 : x2);
                    }

                    if (rise > 0 && set > 0)
                    {
                        break;
                    }

                    h0 = h2;
                }
            }

            risem = DateTime.MinValue;
            setm = DateTime.MinValue;

            if (rise > 0)
            {
                risem = HoursLater(t, rise);
            }

            if (set > 0)
            {
                setm = HoursLater(t, set);
            }

            alwaysUp = null;
            alwaysDown = null;

            if (rise < 0 && set < 0)
            {
                if (ye > 0)
                {
                    alwaysUp = true;
                    alwaysDown = false;
                    risem = DateTime.MaxValue;
                    setm = DateTime.MaxValue;
                }
                else
                {
                    alwaysDown = true;
                    alwaysUp = false;
                }
            }
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="year"></param>
        /// <param name="month"></param>
        /// <param name="day"></param>
        /// <param name="lat">纬度</param>
        /// <param name="lng">经度</param>
        /// <param name="risem"></param>
        /// <param name="setm"></param>
        /// <param name="alwaysUp"></param>
        /// <param name="alwaysDown"></param>
        public static void MoonRisetInt(int year, int month, int day, double lat, double lng, out DateTime risem, out DateTime setm,
                out bool? alwaysUp, out bool? alwaysDown)
        {
            DateTime dt = new DateTime(year, month, day, 0, 0, 0, DateTimeKind.Utc);

            DateTime t = dt;

            double hc = 0.133 * rad;
            double h0 = GetMoonPosition(t, lat, lng).altitude - hc;
            double h1, h2, rise = 0, set = 0, a, b, xe, ye = 0, d, x1, x2, dx;
            int roots;

            for (double i = 1.0; i <= 24.0; i += 2.0)
            {
                h1 = GetMoonPosition(HoursLater(t, i), lat, lng).altitude - hc;
                h2 = GetMoonPosition(HoursLater(t, i + 1), lat, lng).altitude - hc;

                a = (h0 + h2) / 2 - h1;
                b = (h2 - h0) / 2;
                xe = -b / (2 * a);
                ye = (a * xe + b) * xe + h1;
                d = b * b - 4 * a * h1;
                roots = 0;

                if (d >= 0)
                {
                    dx = Math.Sqrt(d) / (Math.Abs(a) * 2);
                    x1 = xe - dx;
                    x2 = xe + dx;
                    if (Math.Abs(x1) <= 1)
                    {
                        roots++;
                    }

                    if (Math.Abs(x2) <= 1)
                    {
                        roots++;
                    }

                    if (x1 < -1)
                    {
                        x1 = x2;
                    }

                    if (roots == 1)
                    {
                        if (h0 < 0)
                        {
                            rise = i + x1;
                        }
                        else
                        {
                            set = i + x1;
                        }
                    }
                    else if (roots == 2)
                    {
                        rise = i + (ye < 0 ? x2 : x1);
                        set = i + (ye < 0 ? x1 : x2);
                    }

                    if (rise > 0 && set > 0)
                    {
                        break;
                    }

                    h0 = h2;
                }
            }

            risem = DateTime.MinValue;
            setm = DateTime.MinValue;

            if (rise > 0)
            {
                risem = HoursLater(t, rise);
            }

            if (set > 0)
            {
                setm = HoursLater(t, set);
            }

            alwaysUp = null; //false;//
            alwaysDown = null; //false;//

            if (rise < 0 && set < 0)
            {
                if (ye > 0)
                {
                    alwaysUp = true;
                    alwaysDown = false;
                    risem = DateTime.MaxValue;
                    setm = DateTime.MaxValue;
                }
                else
                {
                    alwaysDown = true;
                    alwaysUp = false;
                }
            }
        }
        #endregion

        /// <summary>
        /// 计算月亮升起和降落时间
        /// </summary>
        /// <param name="date"></param>
        /// <param name="latitude">纬度</param>
        /// <param name="longitude">经度</param>
        /// <returns></returns>
        public static MoonTimeResult GetMoonTime(DateTime date, double latitude, double longitude)
        {
            double moonrise = 0;
            double moonset = 0;
            DateTime dt = new DateTime(date.Year, date.Month, date.Day, 0, 0, 0, DateTimeKind.Utc);
            bool? up = null;
            bool? down = null;
            DateTime rise;
            DateTime set;
            MoonRiset(dt, latitude, longitude, out rise, out set, out up, out down); //, out up, out down
            DateTime moonriseTime = rise;// ToLocalTime(rise, moonrise);// 
            DateTime moonsetTime = set;//ToLocalTime(set, moonset); // 
            return new MoonTimeResult(moonriseTime, moonsetTime);
        }

        /// <summary>
        /// 私有方法：将 UTC 时间转换为本地时间
        /// </summary>
        /// <param name="time"></param>
        /// <param name="utTime"></param>
        /// <returns></returns>

        private static DateTime ToLocalTime(DateTime time, double utTime)
        {
            int hour = Convert.ToInt32(Math.Floor(utTime));
            double temp = utTime - hour;
            hour += 8; // 转换为东8区北京时间
            temp = temp * 60;
            int minute = Convert.ToInt32(Math.Floor(temp));
            try
            {
                return new DateTime(time.Year, time.Month, time.Day, hour, minute, 0);
            }
            catch
            {
                return new DateTime(time.Year, time.Month, time.Day, 0, 0, 0);
            }
        }




    }


    /// <summary>
    /// 日出日落时间结果类
    /// </summary>
    public class SunTimeResult
    {
        /// <summary>
        /// 
        /// </summary>
        public DateTime SunriseTime { get; set; }
        /// <summary>
        /// 
        /// </summary>
        public DateTime SunsetTime { get; set; }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="sunrise"></param>
        /// <param name="sunset"></param>
        public SunTimeResult(DateTime sunrise, DateTime sunset)
        {
            SunriseTime = sunrise;
            SunsetTime = sunset;
        }
    }

    /// <summary>
    /// 月亮升起和降落时间结果类
    /// </summary>
    public class MoonTimeResult
    {

        /// <summary>
        /// 
        /// </summary>
        public DateTime MoonriseTime { get; set; }
        /// <summary>
        /// 
        /// </summary>
        public DateTime MoonsetTime { get; set; }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="moonrise"></param>
        /// <param name="moonset"></param>
        public MoonTimeResult(DateTime moonrise, DateTime moonset)
        {
            MoonriseTime = moonrise;
            MoonsetTime = moonset;
        }
    }

}