Simple sunset/sunrise algorithm

include/sun.h

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <math.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/time.h>
+#include <time.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
+
+/* 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);

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;
+}
+*/