5
$\begingroup$

Is there a simple method of determining, given a UTC date/time, whether it is day or night at a given lat/long coordinate?

I am currently using a formula based on a Sunrise/Sunset Algorithm from the US Naval Observatory, but unfortunately I am having trouble applying to arrive at a simple night/day indicator.

Should I build upon this formula or is there some more efficient way of doing this that I am missing?

Disclaimer: I am a programmer, not an astronomer, so I find the formulas for astronomical calculations a bit confusing.

Extra stuff: Ideally I would like to artificially shorten the period of "night" by two hours, for the purposes of eliminating some false positives from sensor data that occurs during twilight/dawn. In other words, if sunset occurs at 7 pm at a given location, I would like the option of padding the value by 60 minutes such that the output remains "day" until 8 pm for that location. (And apply the same padding prior to sunrise.)

Update:

I wound up calculating the sun's altitude for the location using a UTC timestamp and the location's lat/long coordinates. Using figures for various twilight amounts (6, 12, 18 or degrees below horizon) I am able to automatically select the portion of "night" that is sufficiently dark. (e.g. 6 degrees below horizon is not as dark as 12 degrees.) This approach was much easier than attempting to determine sunrise and sunset times for the locale's time zone (and all of the associated time zone headaches).

$\endgroup$
5
  • $\begingroup$ The link in your question is dead. It seems to have moved to edwilliams.org/sunrise_sunset_algorithm.htm as on 2019/08 $\endgroup$
    – Waslap
    Aug 31, 2019 at 4:34
  • $\begingroup$ @Waslap Thanks, I've updated the link. $\endgroup$
    – JYelton
    Sep 4, 2019 at 22:17
  • $\begingroup$ Can you share your code? $\endgroup$
    – Ori
    Feb 20, 2022 at 14:50
  • $\begingroup$ @Ori This was 10 years ago. I don't have access to the code base any more. However if I am able to find my notes I will update the question. $\endgroup$
    – JYelton
    Feb 21, 2022 at 15:06
  • $\begingroup$ @Ori I've found some of the original code, and posted it as an additional answer. I do not guarantee its functionality nor completeness, but it may be of some help to others looking for ideas. $\endgroup$
    – JYelton
    Feb 22, 2022 at 22:37

2 Answers 2

7
$\begingroup$

Formulas on that page are about as simple as you'll find and should meet your needs. If you're confused by something specific, just ask!

For shortening the period of night, you should go by the distance (in degrees, say) the Sun is below the horizion, rather than any fixed time. That is the setting of the "zenith" variable. You might like http://en.wikipedia.org/wiki/Twilight for background, then the values "official", "civil", "nautical", and "astronomical" on the page you referenced should make sense. Without even knowing the use for your sensor, I would suggest starting with the nautical twilight value, that is set the variable zenith to 102. (Nautical twilight is pretty dark!)

UPD: The math on the page looks good. I tried coding it and it gave the right answer (for my location today anyway) to within about a minute.

$\endgroup$
1
  • $\begingroup$ Thanks for pointing out the usage of zenith, this may actually be more useful than creating a buffer, as it essentially means "at what point do we consider the sun to have set/risen." I have been using official (90°) thus far, and that's probably the reason I encounter problems -- it's not really dark yet! $\endgroup$
    – JYelton
    Apr 5, 2012 at 21:52
2
$\begingroup$

As requested by @Ori, I found the code that was used for this purpose. I believe a newer version of this code eliminated a lot of extra information, since only the Altitude was ultimately needed to determine the number of degrees above/below the horizon that the sun was for a given location.

There are a lot of comments which should hopefully help anyone who is new to coding. This was written in C# (.NET) but should be relatively easy to adapt to other languages. I believe this code can be greatly simplified and reduced, but I am sharing it in its original form because it may be helpful to adapt to other purposes.

Finally, credit given where due, to Paul Schlyter whose web page was very helpful, and referenced as a link in the comments as well.

public static class SunCalc
{
    public struct SunData
    {
        public double Azimuth;
        public double Altitude;
        public double RightAscension;
        public double Declination;
    }

    public static SunData GetSunPosition(DateTime utcTime, double locationLatitude, double locationLongitude)
    {
        var dayNumber = GetDayNumber(utcTime);
        var argumentOfPerihelion = Sun_ArgumentOfPerihelion(dayNumber);
        var eclipticObliquity = EclipticObliquity(dayNumber);
        /*
         * First, compute the eccentric anomaly E from the mean anomaly M and from the eccentricity e (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) )
         *      
         */
        var meanAnomaly = Sun_MeanAnomaly(dayNumber);
        var eccentricity = Sun_Eccentricity(dayNumber);
        var eccentricAnomaly = meanAnomaly + (180 / Math.PI) * eccentricity * Math.Sin(Deg2Rad(meanAnomaly)) * (1.0 + eccentricity * Math.Cos(Deg2Rad(meanAnomaly)));

        /*
         * 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 )
         * 
         * (note that the r computed here is later used as rs)
         */
        var xv = Math.Cos(Deg2Rad(eccentricAnomaly)) - eccentricity;
        var yv = Math.Sqrt(1.0 - eccentricity * eccentricity) * Math.Sin(Deg2Rad(eccentricAnomaly));
        var v = Rad2Deg(Math.Atan2(yv, xv));
        var r = Math.Sqrt(xv * xv + yv * yv);

        /*
         * Now, compute the Sun's true longitude:
         * 
         *      lonsun = v + w
         */
        var sunTrueLongitude = Rev(v + argumentOfPerihelion);

        /*
         * Convert lonsun, r to ecliptic rectangular geocentric coordinates xs,ys:
         * 
         *     xs = r * cos(lonsun)
         *     ys = r * sin(lonsun)
         */
        var xs = r * Math.Cos(Deg2Rad(sunTrueLongitude));
        var ys = r * Math.Sin(Deg2Rad(sunTrueLongitude));

        /*
         * (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.
         * To convert this to equatorial, rectangular, geocentric coordinates, compute:
         * 
         *     xe = xs
         *     ye = ys * cos(ecl)
         *     ze = ys * sin(ecl)
         */
        var xe = xs;
        var ye = ys * Math.Cos(Deg2Rad(eclipticObliquity));
        var ze = ys * Math.Sin(Deg2Rad(eclipticObliquity));

        /*
         * Finally, compute the Sun's Right Ascension (RA) and Declination (Dec):
         *     RA  = atan2( ye, xe )
         *     Dec = atan2( ze, sqrt(xe*xe+ye*ye) )
         */
        var rightAscension = Rad2Deg(Math.Atan2(ye, xe));
        var declination = Rad2Deg(Math.Atan2(ze, Math.Sqrt(xe * xe + ye * ye)));

        /*
         * Calculate Greenwich Sidereal Time, Sidereal Time and the Sun's Hour Angle
         */
        var sunMeanLongitude = Sun_MeanLongitude(dayNumber);
        var gmst0 = sunMeanLongitude / 15 + 12;

        var siderealTime = RevTime(gmst0 + utcTime.Hour + (utcTime.Minute / 60F) + locationLongitude / 15);
        var hourAngle = RevTime(siderealTime - rightAscension / 15);

        /*
         * Convert the Sun's Hour Angle and Declination to a rectangular coordinate system where the X
         * axis points to the celestial equator in the south, the Y axis to the horizon in the west,
         * and the Z axis to the north celestial pole.
         */
        var x = Math.Cos(Deg2Rad(hourAngle * 15)) * Math.Cos(Deg2Rad(declination));
        var y = Math.Sin(Deg2Rad(hourAngle * 15)) * Math.Cos(Deg2Rad(declination));
        var z = Math.Sin(Deg2Rad(declination));

        /*
         * Rotate this x,y,z axis system along an axis going east-west (Y axis) in such a way that the
         * Z axis will point to the zenith. At the North Pole, the angle of rotation will be zero since
         * there the north celestial pole already is in the zenith. At other latitudes the angle of
         * rotation becomes 90 - latitude.
         */
        var xhor = x * Math.Sin(Deg2Rad(locationLatitude)) - z * Math.Cos(Deg2Rad(locationLatitude));
        var yhor = y;
        var zhor = x * Math.Cos(Deg2Rad(locationLatitude)) + z * Math.Sin(Deg2Rad(locationLatitude));

        /*
         * Compute azimuth and altitude.
         */
        var azimuth = Rad2Deg(Math.Atan2(yhor, xhor)) + 180;
        var altitude = Rad2Deg(Math.Asin(zhor));

        return new SunData
        {
            RightAscension = rightAscension,
            Declination = declination,
            Azimuth = azimuth,
            Altitude = altitude
        };
    }

    private static double GetDayNumber(DateTime dt)
    {
        // http://www.stjarnhimlen.se/comp/ppcomp.html#5
        /*
         * The time scale in these formulae are counted in days. Hours, minutes, seconds are expressed as fractions of a day.
         * Day 0.0 occurs at 2000 Jan 0.0 UT (or 1999 Dec 31, 0:00 UT).
         * This "day number" d is computed as follows (y=year, m=month, D=date, UT=UT in hours+decimals):
         * 
         *      d = 367*y - 7 * ( y + (m+9)/12 ) / 4 + 275*m/9 + D - 730530
         *      
         * Note that ALL divisions here should be INTEGER divisions.
         * Finally, include the time of the day, by adding:
         * 
         *      d = d + UT/24.0        (this is a floating-point division)
         *      
         */

        var d = 367 * dt.Year - 7 * (dt.Year + (dt.Month + 9) / 12) / 4 + 275 * dt.Month / 9 + dt.Day - 730530;
        double hm = dt.Hour + (dt.Minute / 60F);
        return d + hm / 24;
    }

    /// <summary>
    /// Longitude of the Ascending Node (N)
    /// </summary>
    private static double Sun_LongitudeOfAscendingNode()
    {
        return 0.0D;
    }

    /// <summary>
    /// Inclination to the Ecliptic (i) (plane of the Earth's orbit)
    /// </summary>
    private static double Sun_InclinationToEcliptic()
    {
        return 0.0D;
    }

    /// <summary>
    /// Argument of Perihelion (w)
    /// </summary>
    private static double Sun_ArgumentOfPerihelion(double dayNumber)
    {
        return 282.9404 + 4.70935e-5 * dayNumber;
    }

    /// <summary>
    /// Semi-major Axis (a) or mean distance from sun, 1.000000 (AU)
    /// </summary>
    private static double Sun_SemiMajorAxis()
    {
        return 1.0D;
    }

    /// <summary>
    /// Eccentricity (e) where 0 = circle, 0-1 = ellipse, and 1 = parabola
    /// </summary>
    private static double Sun_Eccentricity(double dayNumber)
    {
        return 0.016709 - 1.151e-9 * dayNumber;
    }

    /// <summary>
    /// Mean anomaly (M) (0 at perihelion; increases uniformly with time)
    /// </summary>
    private static double Sun_MeanAnomaly(double dayNumber)
    {
        return Rev(356.0470 + 0.9856002585 * dayNumber);
    }

    /// <summary>
    /// Ecliptic Obliquity (ecl)
    /// </summary>
    private static double EclipticObliquity(double dayNumber)
    {
        return 23.4393 - 3.563e-7 * dayNumber;
    }

    /// <summary>
    /// Longitude of Perihelion (w1)
    /// </summary>
    private static double Sun_LongitudeOfPerihelion(double dayNumber)
    {
        var n = Sun_LongitudeOfAscendingNode();
        var w = Sun_ArgumentOfPerihelion(dayNumber);
        return Rev(n + w);
    }

    /// <summary>
    /// Mean Longitude (L)
    /// </summary>
    private static double Sun_MeanLongitude(double dayNumber)
    {
        var m = Sun_MeanAnomaly(dayNumber);
        var w1 = Sun_LongitudeOfPerihelion(dayNumber);
        return Rev(m + w1);
    }

    /// <summary>
    /// Perihelion Distance (q)
    /// </summary>
    private static double Sun_PerihelionDistance(double dayNumber)
    {
        var a = Sun_SemiMajorAxis();
        var e = Sun_Eccentricity(dayNumber);
        return a * (1 - e);
    }

    /// <summary>
    /// Aphelion Distance (Q)
    /// </summary>
    private static double Sun_AphelionDistance(double dayNumber)
    {
        var a = Sun_SemiMajorAxis();
        var e = Sun_Eccentricity(dayNumber);
        return a * (1 + e);
    }

    /// <summary>
    /// Convert degrees to radians
    /// </summary>
    private static double Deg2Rad(double angleDegrees)
    {
        return angleDegrees * (Math.PI / 180.0);
    }

    /// <summary>
    /// Convert radians to degrees
    /// </summary>
    private static double Rad2Deg(double angleRadians)
    {
        return angleRadians * (180.0 / Math.PI);
    }

    /// <summary>
    /// Revolution function, normalizes an angle to between 0 and 360 degrees by adding or subtracting even multiples of 360.
    /// </summary>
    private static double Rev(double x)
    {
        return x - Math.Floor(x / 360.0) * 360.0;
    }

    /// <summary>
    /// Revolution function, normalizes a time value (in hours) to between 0 and 24 by adding or subtracting even multiples of 24.
    /// </summary>
    private static double RevTime(double x)
    {
        return x - Math.Floor(x / 24.0) * 24.0;
    }
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.