diff --git a/src/sun.c b/src/sun.c new file mode 100644 index 0000000..5876b26 --- /dev/null +++ b/src/sun.c @@ -0,0 +1,159 @@ +/** + * Calculation routines + * + * @copyright 2012 Steffen Vogel + * @license http://www.gnu.org/licenses/gpl.txt GNU Public License + * @author Steffen Vogel + * @link http://www.steffenvogel.de/2012/03/14/cron-jobs-fur-sonnenauf-untergang/ + */ +/* + * This file is part of sun + * + * sun is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * any later version. + * + * sun is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with sun. If not, see . + */ + +#include +#include + +#include "sun.h" + +double deg2rad(double deg) { + return M_PI * deg / 180.0; +} + +double rad2deg(double rad) { + return 180 * rad / M_PI; +} + +int julian_date(struct tm date) { + int year = date.tm_year + 1900; + int month = date.tm_mon + 1; + int day = date.tm_mday; + + int a = (14 - month) / 12; + int y = year + 4800 - a; + int m = month + 12 * a - 3; + int jdn = day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045; + + return jdn; +} + +double in_pi(double x) { + x -= floor(x / M_2PI) * M_2PI; + + return (x < 0) ? x + M_2PI : x; +} + +double sun_eps(double days) { + double poly; + + poly = -46.8150 * days - 0.00059 * powf(days, 2) + 0.001813 * powf(days, 3); + poly /= 3600.0; + + return deg2rad(23.43929111 + poly); +} + +double sun_eof(double days, double ra) { + double ra_mean = 18.71506921 + 2400.0513369 * days + (2.5862e-5 - 1.72e-9 * days) * powf(days, 2); + ra_mean = 24.0 * in_pi(M_2PI * ra_mean / 24.0) / M_2PI; + double ra_diff = ra_mean - ra; + + if (ra_diff < -12) ra_diff += 24; + if (ra_diff > 12) ra_diff -= 24; + + ra_diff *= 1.0027379; + +#ifdef DEBUG + printf("ra_mean = %.20f\n", ra_mean); + printf("ra_diff = %.20f\n", ra_diff); +#endif + + return ra_diff; +} + +struct sun_coords sun_pos(double days) { + double m = in_pi(M_2PI * (0.993133 + 99.997361 * days)); + double l = in_pi(M_2PI * (0.7859453 + m / M_2PI + (6893 * sin(m) + 72.0 * sin(2 * m) + 6191.2 * days) / 1296.0e3)); + double e = sun_eps(days); + double ra = atan(tan(l)*cos(e)); + double dk = asin(sin(e) * sin(l)); + + if (ra < 0) ra += M_PI; + if (l > M_PI) ra += M_PI; + + ra *= 12 / M_PI; + +#ifdef DEBUG + printf("m = %.20f\n", m); + printf("l = %.20f\n", l); + printf("e = %.20f\n", e); + printf("dk = %.20f\n", dk); + printf("ra = %.20f\n", ra); +#endif + + struct sun_coords spos = { dk, ra }; + return spos; +} + +struct tm sun(struct coords pos, struct tm date, enum mode mode, double twilight, int timezone) { + /* calculate sun position */ + double jd = julian_date(date); + double days = (jd - 2451545) / 36525; + struct sun_coords spos = sun_pos(days); + + /* use equation of time */ + double h = deg2rad(twilight); + double time = sun_eof(days, spos.ra); + double time_diff = 12.0 * acos((sin(h) - sin(deg2rad(pos.lat)) * sin(spos.dk)) / (cos(deg2rad(pos.lat)) * cos(spos.dk))) / M_PI; + +#ifdef DEBUG + printf("jd = %.20f\n", jd); + printf("days = %.20f\n", days); + printf("delta_days = %.20f\n", jd - 2451545); + printf("h = %.20f\n", h); + printf("time = %.20f\n", time); + printf("time_diff = %.20f\n", time_diff); +#endif + + double local = 12 - time; + double global = local - pos.lon / 15 + timezone; + + switch (mode) { + case NOON: /* do nothing */ break; + case RISE: global -= time_diff; break; + case SET: global += time_diff; break; + case DAYTIME: global = 2 * time_diff; break; + case NIGHTTIME: global = 24 - 2 * time_diff; break; + } + + + if (global < 0) { + global += 24; + } + else if (global >= 24) { + global -= 24; + } + +#ifdef DEBUG + printf("local = %f\n", local); + printf("global = %f\n", global); +#endif + + date.tm_hour = (int) global; + date.tm_min = 60 * (global - date.tm_hour) + 0.5; /* math rounding */ + date.tm_sec = 0; /* below our accuracy */ + //date.tm_sec = 3600 * (global - date.tm_hour - date.tm_min / 60.0); + + return date; +} diff --git a/src/sun.h b/src/sun.h new file mode 100644 index 0000000..36132dd --- /dev/null +++ b/src/sun.h @@ -0,0 +1,55 @@ +/** + * Header file for calculation + * + * @copyright 2012 Steffen Vogel + * @license http://www.gnu.org/licenses/gpl.txt GNU Public License + * @author Steffen Vogel + * @link http://www.steffenvogel.de/2012/03/14/cron-jobs-fur-sonnenauf-untergang/ + */ +/* + * This file is part of sun + * + * sun is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * any later version. + * + * sun is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with sun. If not, see . + */ + + +#include + +#ifdef GEONAMES_SUPPORT +#include "geonames.h" +#else +struct coords { double lat, lon; }; +#endif + +struct sun_coords { double dk, ra; }; + +//#define DEBUG 1 +#define GEONAMES_SUPPORT 1 +#define M_2PI (M_PI * 2) + +enum mode { RISE, SET, NOON, DAYTIME, NIGHTTIME }; + +/* helper function */ +void usage(); +double deg2rad(double deg); +double rad2deg(double rad); +double in_pi(double x); +int julian_date(struct tm date); + +/* calculation functions */ +struct tm sun(struct coords pos, struct tm date, enum mode mode, double twilight, int timezone); +struct sun_coords sun_pos(double days); +double sun_eof(double days, double ra); +double sun_eps(double days); + diff --git a/src/sun_main.c b/src/sun_main.c new file mode 100644 index 0000000..6f39a36 --- /dev/null +++ b/src/sun_main.c @@ -0,0 +1,286 @@ +/** + * Main routine + * + * Does parsing of command line options and start calculation + * + * @copyright 2012 Steffen Vogel + * @license http://www.gnu.org/licenses/gpl.txt GNU Public License + * @author Steffen Vogel + * @link http://www.steffenvogel.de/2012/03/14/cron-jobs-fur-sonnenauf-untergang/ + */ +/* + * This file is part of sun + * + * sun is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * any later version. + * + * sun is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with sun. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include + +#include "sun.h" + +static struct option long_options[] = { + {"twilight", required_argument, 0, 't'}, + {"date", required_argument, 0, 'd'}, + {"format", required_argument, 0, 'f'}, + {"lat", required_argument, 0, 'a'}, + {"lon", required_argument, 0, 'o'}, +#ifdef GEONAMES_SUPPORT + {"query", required_argument, 0, 'q'}, +#endif + {"zone", required_argument, 0, 'z'}, + {"help", no_argument, 0, 'h'}, + {0} +}; + +static char *long_options_descs[] = { + "use special twilight (nautic|civil|astro)", + "calculcate for specified date (eg. 2011-12-25)", + "output format (eg. %H:%M:%S)", + "geographical latitude (-90° to 90°)", + "geographical longitude (-180° to 180°)", +#ifdef GEONAMES_SUPPORT + "query geonames.org for geographical position", +#endif + "use timezone for output", + "show this help" +}; + +void usage() { + printf("usage: sun mode [options]\n\n"); + printf(" mode is one of: rise, set, noon, daytime, nighttime\n\n"); + printf(" following OPTIONS are available\n"); + + struct option *op = long_options; + char **desc = long_options_descs; + while (op->name && desc) { + printf("\t-%c, --%s\t%s\n", op->val, op->name, *desc); + op++; + desc++; + } + + printf("\nA combination of --lat, --lon or --query is required!\n"); + printf("For bugs or feature requests please contact Steffen Vogel \n"); +} + +char * strreplace(char *subject, char *search, char *replace) { + int new_len = strlen(subject); + int search_len = strlen(search); + int replace_len = strlen(replace); + char *tmp; + + for (tmp = strstr(subject, search); tmp != NULL; tmp = strstr(tmp + search_len, search)) { + new_len += replace_len - search_len; + } + + char *old = subject; + char *new = malloc(new_len); + + new[0] = '\0'; /* empty string */ + for (tmp = strstr(subject, search); tmp != NULL; tmp = strstr(tmp + search_len, search)) { + new_len = strlen(new); + + strncpy(new + new_len, old, tmp - old); + strcpy(new + new_len + (tmp - old), replace); + old = tmp + search_len; + } + + strcpy(new + strlen(new), old); + + return new; +} + +int main(int argc, char *argv[]) { + /* default options */ + double twilight = -50.0 / 60; /* 50 Bogenminuten; no twilight, normal sunset/rise */ + int timezone = 0; /* UTC */ + char *format = "%H:%M"; + char *query = NULL; + bool error = false; + + enum mode mode; + struct coords pos = { INFINITY, INFINITY }; + struct tm date; + + /* default date: now */ + time_t t; + time(&t); + localtime_r(&t, &date); + + /* parse command */ + if (argc <= 1) { + fprintf(stderr, "mode required\n\n"); + error = true; + } + else if (strcmp(argv[1], "rise") == 0) mode = RISE; + else if (strcmp(argv[1], "set") == 0) mode = SET; + else if (strcmp(argv[1], "noon") == 0) mode = NOON; + else if (strcmp(argv[1], "daytime") == 0) mode = DAYTIME; + else if (strcmp(argv[1], "nighttime") == 0) mode = NIGHTTIME; + else { + fprintf(stderr, "invalid mode: %s\n", argv[1]); + error = true; + } + + /* parse command line arguments */ + while (1) { + int optidx; + int c = getopt_long(argc-1, argv+1, "ht:d:f:a:o:q:z:", long_options, &optidx); + + /* detect the end of the options. */ + if (c == -1) break; + + switch (c) { + case 't': + if (strcmp(optarg, "civil") == 0) { + twilight = -6.0; + } + else if (strcmp(optarg, "nautic") == 0) { + twilight = -12.0; + } + else if (strcmp(optarg, "astro") == 0) { + twilight = -18.0; + } + else { + fprintf(stderr, "invalid twilight: %s\n", optarg); + error = true; + } + break; + + case 'd': + if (strptime(optarg, "%Y-%m-%d", &date) == 0) { + fprintf(stderr, "invalid date: %s\n", optarg); + error = true; + } + break; + + case 'f': + format = strdup(optarg); + break; + + case 'a': + pos.lat = strtod(optarg, NULL); + break; + + case 'o': + pos.lon = strtod(optarg, NULL); + break; +#ifdef GEONAMES_SUPPORT + case 'q': + query = strdup(optarg); + break; +#endif + + case 'z': + timezone = atoi(optarg); + break; + + case 'h': + usage(); + return EXIT_SUCCESS; + + case '?': + default: + error = true; + } + } + +#ifdef GEONAMES_SUPPORT + if (query && geonames_lookup(query, &pos, NULL, 0) != 0) { + fprintf(stderr, "failed to lookup location: %s\n", query); + error = true; + } +#endif + + if (pos.lat == INFINITY) { + fprintf(stderr, "please provide a complete position: latitude is missing\n"); + error = true; + } + else { + if (fabs(pos.lat) > 90) { + fprintf(stderr, "invalid latitude: %.4f\n", pos.lat); + error = true; + } + else if (fabs(pos.lat) > 65) { + fprintf(stderr, "attention: results may be inaccurate! (delta_T > 5min)\n"); + } + } + + if (pos.lon == INFINITY) { + fprintf(stderr, "please provide a complete position: longitude is missing\n"); + error = true; + } + else if (fabs(pos.lon) > 180) { + fprintf(stderr, "invalid longitude: %.4f\n", pos.lon); + error = true; + } + + if (error) { + printf("\n"); + usage(); + return EXIT_FAILURE; + } + +#ifdef DEBUG + char date_str[64]; + strftime(date_str, 64, "%Y-%m-%d", &date); + printf("calculate for %s\n", date_str); + printf("for position: %f, %f\n", pos.lat, pos.lon); + printf("with twilight: %f\n", twilight); +#endif + + /* start the calculation */ + struct tm result = sun(pos, date, mode, twilight, timezone); + + char result_str[64]; + + /* workaround for unix timestamps and 'struct tm's limitation on years > 1990 */ + if (mode == DAYTIME || mode == NIGHTTIME) { + if (strstr(format, "%s") != NULL) { + result.tm_year = 70; + result.tm_mon = 0; + result.tm_mday = 1; + result.tm_hour++; + + char buffer[16]; + sprintf(buffer, "%lu", mktime(&result)); + +#ifdef DEBUG + printf("buffer: %s\n", buffer); +#endif + + format = strreplace(format, "%s", buffer); + + result.tm_hour--; + +#ifdef DEBUG + printf("format: %s\n", format); +#endif + } + + result.tm_year = -1900; + result.tm_mon = -1; + result.tm_mday = 0; + } + + strftime(result_str, 64, format, &result); + printf("%s\n", result_str); + + return EXIT_SUCCESS; +}