CSharp: SunTimeCalculator

using System;
using System.Collections.Generic;
using System.Linq;
using System.Web;

namespace WebAppPdfDemo
{



    /// <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);
        }

        /// <summary>
        /// 计算月亮升起和降落时间
        /// </summary>
        /// <param name="date"></param>
        /// 
        /// 
        /// <param name="longitude"></param>
        /// <param name="latitude"></param>
        /// <returns></returns>
        public static MoonTimeResult GetMoonTime(DateTime date, double longitude, double latitude)
        {
            double moonrise = 0;
            double moonset = 0;
            MoonRiset(date.Year, date.Month, date.Day, longitude, latitude, -35.0 / 60.0, 1, ref moonrise, ref moonset);
            DateTime moonriseTime = ToLocalTime(date, moonrise);
            DateTime moonsetTime = ToLocalTime(date, 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; }

        public MoonTimeResult(DateTime moonrise, DateTime moonset)
        {
            MoonriseTime = moonrise;
            MoonsetTime = moonset;
        }
    }




}