From 5715977e62e4aa9c864b8c4518866ea0e2e3e774 Mon Sep 17 00:00:00 2001 From: Pim van Pelt Date: Wed, 2 Oct 2019 23:55:32 +0200 Subject: [PATCH] Simple sunset/sunrise algorithm --- include/sun.h | 35 +++++++++++ src/sun.c | 163 ++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 include/sun.h create mode 100644 src/sun.c diff --git a/include/sun.h b/include/sun.h new file mode 100644 index 0000000..c7404ba --- /dev/null +++ b/include/sun.h @@ -0,0 +1,35 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +// See http://edwilliams.org/sunrise_sunset_algorithm.htm +#define ZENITH_OFFICIAL 90.83 +#define ZENITH_CIVIL 96.0 +#define ZENITH_NAUTICAL 102.0 +#define ZENITH_ASTRONOMICAL 108.0 + +/* Given the latlong, returns the seconds in the day, in UTC, that the sun rises + * on a given year-month-day, or 0 if it never rises. + */ +uint32_t calcSunRise(int year, int month, int day, float latitude, float longitude); + +/* Given the latlong, returns the seconds in the day, in UTC, that the sun sets + * on a given year-month-day, or 0 if it never sets. + */ +uint32_t calcSunSet(int year, int month, int day, float latitude, float longitude); + +/* Given a month and hemisphere, a sunrise and sunset (output from + * calcSun{Rise,Set}()), returns the amount of seconds of the day where + * daylight is observed. + * + * Is reasonably aware of polar regions (ie permanent day/night). Exception is the + * first full day/night, which are truncated to 86400 respectively 0 seconds. This is + * only a concern for any latitude > 66.6 or latitude < -66.6. + */ +uint32_t calcDaylight(int month, float latitude, uint32_t rise, uint32_t set); diff --git a/src/sun.c b/src/sun.c new file mode 100644 index 0000000..b0a1ddb --- /dev/null +++ b/src/sun.c @@ -0,0 +1,163 @@ +#include "sun.h" + +// See http://edwilliams.org/sunrise_sunset_algorithm.htm +#define ZENITH_OFFICIAL 90.83 +#define ZENITH_CIVIL 96.0 +#define ZENITH_NAUTICAL 102.0 +#define ZENITH_ASTRONOMICAL 108.0 + +static float calcSun(int year, int month, int day, float latitude, + float longitude, bool sunset, float zenith) { + int N1 = floor(275 * month / 9); + int N2 = floor((month + 9) / 12); + int N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3)); + int N = N1 - (N2 * N3) + day - 30; + + float lngHour = longitude / 15.; + + float t; + if (sunset == false) { + // if rising time is desired: + t = N + ((6 - lngHour) / 24); + } else { + // if setting time is desired: + t = N + ((18 - lngHour) / 24); + } + + float M = (0.9856 * t) - 3.289; + + // calculate the Sun's true longitude + float L = M + (1.916 * sin((M_PI / 180.0f) * M)) + + (0.020 * sin((M_PI / 180.0f) * 2 * M)) + 282.634; + if (L < 0) L += 360.0f; + if (L > 360) L -= 360.0f; + + // 5a. calculate the Sun's right ascension + float RA = (180.0f / M_PI) * atan(0.91764 * tan((M_PI / 180.0f) * L)); + if (RA < 0) RA += 360; + if (RA > 360) RA -= 360; + + // 5b. right ascension value needs to be in the same quadrant as L + float Lquadrant = (floor(L / 90)) * 90; + float RAquadrant = (floor(RA / 90)) * 90; + RA = RA + (Lquadrant - RAquadrant); + + // 5c. right ascension value needs to be converted into hours + RA = RA / 15.; + + // 6. calculate the Sun's declination + float sinDec = 0.39782 * sin((M_PI / 180.0f) * L); + float cosDec = cos(asin(sinDec)); + + // 7a. calculate the Sun's local hour angle + float cosH = (cos((M_PI / 180.0f) * zenith) - + (sinDec * sin((M_PI / 180.0f) * latitude))) / + (cosDec * cos((M_PI / 180.0f) * latitude)); + + if (cosH > 1) { + return 0; + } else if (cosH < -1) { + return 0; + } + + // 7b. finish calculating H and convert into hours + float H; + if (sunset == false) { + // if rising time is desired: + H = 360 - (180.0f / M_PI) * acos(cosH); + } else { + // if setting time is desired: + H = (180.0f / M_PI) * acos(cosH); + } + H = H / 15.; + + // 8. calculate local mean time of rising/setting + float T = H + RA - (0.06571 * t) - 6.622; + + // 9. adjust back to UTC + float UT = T - lngHour; + if (UT < 0) UT += 24; + if (UT > 24) UT -= 24; + + return UT; +} + +uint32_t calcSunRise(int year, int month, int day, float latitude, + float longitude) { + return (uint32_t)(3600 * calcSun(year, month, day, latitude, longitude, false, + ZENITH_OFFICIAL)) % + 86400; +} + +uint32_t calcSunSet(int year, int month, int day, float latitude, + float longitude) { + return (uint32_t)(3600 * calcSun(year, month, day, latitude, longitude, true, + ZENITH_OFFICIAL)) % + 86400; +} + +uint32_t calcDaylight(int month, float latitude, uint32_t rise, uint32_t set) { + if (set == 0 || rise == 0) { + if (month < 4 || month > 9) { + if (latitude < 0) return 86400; + else return 0; + } else { + if (latitude > 0) return 86400; + else return 0; + } + } + + if (set < rise) return (86400 - (rise - set)) % 86400; + return (set - rise) % 86400; +} + +// str needs to be at least 20 characters +/* +static void clk(uint32_t seconds, char *str) { + if (!str) return; + snprintf(str, 20, "%u %02d:%02d:%02d", seconds, seconds / 3600, + (seconds % 3600) / 60, seconds % 60); +} + +int main(int argc, char **argv) { + int year = 2004, month = 8, day = 21; + double latitude = 39.95; // convert to just degrees. No min/sec + double longitude = 75.15; + uint32_t set, rise, daylight; + + if (argc == 6) { + year = atoi(argv[1]); + month = atoi(argv[2]); + day = atoi(argv[3]); + latitude = atof(argv[4]); + longitude = atof(argv[5]); + + } else { + printf("\nUsage: (note convert latitude and longitude to degrees) \n"); + printf("./sunrise year month date latitude longitude \n\n"); + printf( + " Latitude and longitude must be in degrees and or fraction of " + "degrees. Not min sec\n"); + printf("\n US Listings Of Latitude and Longitude\n"); + printf(" http://www.census.gov/geo/www/tiger/latlng.txt\n"); + + printf("\n Example: (Just outside Philadelphia PA, USA)\n\n"); + printf("./sunrise 2004 8 21 39.95 -75.15 \n"); + printf("./sunrise `date \"+%cY %cm %cd\"` 39.95 -75.15\n\n", '%', '%', '%'); + } + + rise = calcSunRise(year, month, day, latitude, longitude); + set = calcSunSet(year, month, day, latitude, longitude); + daylight = calcDaylight(month, latitude, rise, set); + + char riseTime[20], setTime[20], daylightTime[20]; + + clk(rise, riseTime); + clk(set, setTime); + clk(daylight, daylightTime); + printf("%d-%02d-%02d UTC Sunrise=%s Sunset=%s Daylight=%s\n", year, month, + day, riseTime, setTime, daylightTime); + + return 0; +} +*/