CSharp: SunCalc for calculating

 

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

namespace WebAppPdfDemo
{

    /// <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 SunMoonCalcs
    {
        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);
        }
        /// <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 SunCalc
        {
            private const double J0 = 0.0009;

            public class AzAlt
            {
                public double azimuth = 0;
                public double altitude = 0;
            }
            /// <summary>
            /// 
            /// </summary>
            /// <param name="dt"></param>
            /// <param name="lat"></param>
            /// <param name="lng"></param>
            /// <returns></returns>
            public static AzAlt GetPosition(DateTime dt, double lat, double lng)
            {
                double lw = rad * -lng;
                double phi = rad * lat;
                double d = JulianDays(dt);
                RaDec c = SunCoords(d);
                double H = SiderealTime(d, lw) - c.ra;
                return new AzAlt { azimuth = Azimuth(H, phi, c.dec), altitude = Altitude(H, phi, c.dec) };
            }
            /// <summary>
            /// 
            /// </summary>
            /// <param name="d"></param>
            /// <param name="lw"></param>
            /// <returns></returns>
            public static double JulianCycle(double d, double lw)
            {
                return Math.Round(d - J0 - lw / (2 * PI));
            }
            /// <summary>
            /// 
            /// </summary>
            /// <param name="Ht"></param>
            /// <param name="lw"></param>
            /// <param name="n"></param>
            /// <returns></returns>
            public static double ApproxTransit(double Ht, double lw, double n)
            {
                return J0 + (Ht + lw) / (2 * PI) + n;
            }
            /// <summary>
            /// 
            /// </summary>
            /// <param name="ds"></param>
            /// <param name="M"></param>
            /// <param name="L"></param>
            /// <returns></returns>
            public static double SolarTransitJ(double ds, double M, double L)
            {
                return J2000 + ds + 0.0053 * Math.Sin(M) - 0.0069 * Math.Sin(2.0 * L);
            }
            /// <summary>
            /// 
            /// </summary>
            /// <param name="h"></param>
            /// <param name="phi"></param>
            /// <param name="d"></param>
            /// <returns></returns>
            public static double HourAngle(double h, double phi, double d)
            {
                return Math.Acos((Math.Sin(h) - Math.Sin(phi) * Math.Sin(d)) / (Math.Cos(phi) * Math.Cos(d)));
            }

            /// <summary>
            /// returns set time for the given sun altitude
            /// </summary>
            /// <param name="h"></param>
            /// <param name="lw"></param>
            /// <param name="phi"></param>
            /// <param name="dec"></param>
            /// <param name="n"></param>
            /// <param name="M"></param>
            /// <param name="L"></param>
            /// <returns></returns>
            public static double GetSetJ(double h, double lw, double phi, double dec, double n, double M, double L)
            {
                double w = HourAngle(h, phi, dec);
                double a = ApproxTransit(w, lw, n);
                return SolarTransitJ(a, M, L);
            }

            /// <summary>
            /// solar disc diameter 
            /// </summary>
            /// <param name="dt"></param>
            /// <param name="lat"></param>
            /// <param name="lng"></param>
            /// <param name="rise"></param>
            /// <param name="set"></param>
            /// <param name="angle"></param>
 
            public static void GetTimes(DateTime dt, double lat, double lng, out DateTime rise, out DateTime set,
                double angle = -0.833)
            {
                double lw = rad * -lng;
                double phi = rad * lat;
                double d = JulianDays(dt);
                double n = JulianCycle(d, lw);
                double ds = ApproxTransit(0, lw, n);

                double M = SolarMeanAnomaly(ds);
                double L = EclipticLongitude(M);
                double dec = Declination(L, 0);

                double Jnoon = SolarTransitJ(ds, M, L);
                double Jset = GetSetJ(angle * rad, lw, phi, dec, n, M, L);
                double Jrise = Jnoon - (Jset - Jnoon);

                rise = double.IsNaN(Jrise) ? DateTime.MinValue : FromJulianDate(Jrise);
                set = double.IsNaN(Jset) ? DateTime.MinValue : FromJulianDate(Jset);
            }
            /// <summary>
            /// 
            /// </summary>
            /// <param name="dt"></param>
            /// <param name="lat"></param>
            /// <param name="lng"></param>
            /// <returns></returns>
            public static List<SunTimeRiseSet> GetTimes(DateTime dt, double lat, double lng)
            {
                List<SunTimeRiseSet> retval = new List<SunTimeRiseSet>();
                DateTime rise, set;
                foreach (SunTime st in SunTimes)
                {
                    GetTimes(dt, lat, lng, out rise, out set, st.Angle);
                    retval.Add(new SunTimeRiseSet
                    {
                        Angle = st.Angle,
                        MorningName = st.MorningName,
                        EveningName = st.EveningName,
                        RiseTime = rise,
                        SetTime = set
                    });
                }

                return retval;
            }
        }
        /// <summary>
        /// 
        /// </summary>
        public class MoonCalc
        {


            /// <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 GetMoonTimes(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;
                    }
                }
            }
        }
    }




}