ported to use libnova for calculations

This commit is contained in:
Steffen Vogel 2013-01-29 13:05:11 +01:00
parent 4bac41cf0b
commit 58463813cb
7 changed files with 115 additions and 308 deletions

View file

@ -1,7 +1,7 @@
bin_PROGRAMS = sun
sun_SOURCES = sun_main.c sun.c
sun_LDADD = -lm
sun_SOURCES = sun_main.c
sun_LDADD = -lm -lnova
if GEONAMES_SUPPORT
noinst_PROGRAMS = geonames

View file

@ -71,9 +71,9 @@ am__geonames_SOURCES_DIST = geonames_main.c geonames.c
geonames_OBJECTS = $(am_geonames_OBJECTS)
am__DEPENDENCIES_1 =
@GEONAMES_SUPPORT_TRUE@geonames_DEPENDENCIES = $(am__DEPENDENCIES_1)
am__sun_SOURCES_DIST = sun_main.c sun.c geonames.c
am__sun_SOURCES_DIST = sun_main.c geonames.c
@GEONAMES_SUPPORT_TRUE@am__objects_1 = geonames.$(OBJEXT)
am_sun_OBJECTS = sun_main.$(OBJEXT) sun.$(OBJEXT) $(am__objects_1)
am_sun_OBJECTS = sun_main.$(OBJEXT) $(am__objects_1)
sun_OBJECTS = $(am_sun_OBJECTS)
@GEONAMES_SUPPORT_TRUE@am__DEPENDENCIES_2 = $(am__DEPENDENCIES_1)
sun_DEPENDENCIES = $(am__DEPENDENCIES_2)
@ -183,8 +183,8 @@ target_alias = @target_alias@
top_build_prefix = @top_build_prefix@
top_builddir = @top_builddir@
top_srcdir = @top_srcdir@
sun_SOURCES = sun_main.c sun.c $(am__append_1)
sun_LDADD = -lm $(am__append_2)
sun_SOURCES = sun_main.c $(am__append_1)
sun_LDADD = -lm -lnova $(am__append_2)
@GEONAMES_SUPPORT_TRUE@geonames_SOURCES = geonames_main.c geonames.c
@GEONAMES_SUPPORT_TRUE@geonames_LDADD = $(DEPS_GEONAMES_LIBS)
@GEONAMES_SUPPORT_TRUE@AM_CFLAGS = $(DEPS_GEONAMES_CFLAGS)
@ -280,7 +280,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geonames.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geonames_main.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/sun.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/sun_main.Po@am__quote@
.c.o:

View file

@ -73,7 +73,7 @@ int geonames_lookup(const char *place, struct coords *result, char *name, int n)
#ifdef GEONAMES_CACHE_SUPPORT
if (geonames_cache_lookup(place, result, name, n) == EXIT_SUCCESS) {
#ifdef DEBUG
printf("using cached entry\n");
printf("using cached geonames entry\n");
#endif
return EXIT_SUCCESS;
}
@ -143,7 +143,7 @@ int geonames_parse(struct json_object *jobj, struct coords *result, char *name,
struct json_object *jobj_place = json_object_array_get_idx(json_object_object_get(jobj, "geonames"), 0);
result->lat = json_object_get_double(json_object_object_get(jobj_place, "lat"));
result->lon = json_object_get_double(json_object_object_get(jobj_place, "lng"));
result->lng = json_object_get_double(json_object_object_get(jobj_place, "lng"));
if (name && n > 0) {
strncpy(name, json_object_get_string(json_object_object_get(jobj_place, "name")), n);
@ -178,7 +178,7 @@ int geonames_cache_lookup(const char *place, struct coords *result, char *name,
for (col = 0, tok = strtok(line, "\t"); tok != NULL; tok = strtok(NULL, "\t")) {
switch (col) {
case 0:
if (strcmp(tok, place) != 0) {
if (strcasecmp(tok, place) != 0) {
continue; /* skip row */
}
break;
@ -188,7 +188,7 @@ int geonames_cache_lookup(const char *place, struct coords *result, char *name,
break;
case 2:
result->lon = strtod(tok, NULL);
result->lng = strtod(tok, NULL);
break;
case 3:
@ -216,7 +216,7 @@ int geonames_cache_store(const char *place, struct coords *result, char *name, i
/* build cache entry */
char line[256];
snprintf(line, sizeof(line), "%s\t%.5f\t%.5f\t%s\n", place, result->lat, result->lon, name);
snprintf(line, sizeof(line), "%s\t%.5f\t%.5f\t%s\n", place, result->lat, result->lng, name);
if (fputs(line, file) == EOF) {
fclose(file);

View file

@ -32,8 +32,8 @@
#define GEONAMES_CACHE_FILE ".geonames.cache" /* in users home dir */
struct coords {
double lng;
double lat;
double lon;
};
int geonames_lookup(const char *place, struct coords *coords, char *name, int n);

160
src/sun.c
View file

@ -1,160 +0,0 @@
/**
* Calculation routines
*
* @copyright 2012 Steffen Vogel
* @license http://www.gnu.org/licenses/gpl.txt GNU Public License
* @author Steffen Vogel <post@steffenvogel.de>
* @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 <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <stdio.h>
#include "../config.h"
#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;
}

View file

@ -1,53 +0,0 @@
/**
* Header file for calculation
*
* @copyright 2012 Steffen Vogel
* @license http://www.gnu.org/licenses/gpl.txt GNU Public License
* @author Steffen Vogel <post@steffenvogel.de>
* @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 <http://www.gnu.org/licenses/>.
*/
#include <time.h>
#ifdef GEONAMES_SUPPORT
#include "geonames.h"
#else
struct coords { double lat, lon; };
#endif
struct sun_coords { double dk, ra; };
#define M_2PI (M_PI * 2)
enum mode { INVALID, 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);

View file

@ -25,6 +25,9 @@
* along with sun. If not, see <http://www.gnu.org/licenses/>.
*/
#define _XOPEN_SOURCE 700
#define EXIT_CIRCUMPOLAR 2
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
@ -32,13 +35,21 @@
#include <getopt.h>
#include <float.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <libnova/libnova.h>
#include <libnova/solar.h>
#include "../config.h"
#include "sun.h"
enum mode { INVALID, RISE, SET, TRANSIT, DAYTIME, NIGHTTIME };
extern long timezone;
static struct option long_options[] = {
{"twilight", required_argument, 0, 't'},
{"horizon", required_argument, 0, 't'},
{"date", required_argument, 0, 'd'},
{"format", required_argument, 0, 'f'},
{"lat", required_argument, 0, 'a'},
@ -53,7 +64,7 @@ static struct option long_options[] = {
};
static char *long_options_descs[] = {
"use special twilight (nautic|civil|astro)",
"use special twilight (nautic|civil|astronomical)",
"calculcate for specified date (eg. 2011-12-25)",
"output format (eg. %H:%M:%S)",
"geographical latitude (-90° to 90°)",
@ -68,11 +79,12 @@ static char *long_options_descs[] = {
void version () {
printf("%s %s\n", PACKAGE_NAME, PACKAGE_VERSION);
printf("libnova %s\n", LIBNOVA_VERSION);
}
void usage() {
printf("Usage:\n sun mode [options]\n\n");
printf(" mode is one of: rise, set, noon, daytime, nighttime\n\n");
printf(" mode is one of: rise, set, transit, daytime, nighttime\n\n");
printf("Options:\n");
struct option *op = long_options;
@ -116,35 +128,32 @@ char * strreplace(char *subject, char *search, char *replace) {
int main(int argc, char *argv[]) {
/* default options */
double twilight = -50.0 / 60; /* 50 Bogenminuten; no twilight, normal sunset/rise */
char *format = "%H:%M";
double horizon = LN_SOLAR_STANDART_HORIZON; /* 50 Bogenminuten; no twilight, normal sunset/rise */
char *format = "%H:%M:%S";
char *query = NULL;
bool error = false;
int timezone = 0;
enum mode mode = INVALID;
struct tm date;
struct coords pos = { INFINITY, INFINITY };
/* default date: now */
time_t t;
time(&t);
localtime_r(&t, &date);
double julian;
struct tm date = { 0 };
struct ln_rst_time rst;
struct ln_lnlat_posn observer = { DBL_MAX, DBL_MAX };
/* default timezone: system */
struct timezone tz;
if (gettimeofday(NULL, &tz) == 0) {
timezone = -tz.tz_minuteswest / 60.0;
}
tzset();
/* default time: now */
julian = ln_get_julian_from_sys();
/* parse mode */
if (argc > 1 && argv[1][0] != '-') {
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], "transit") == 0) mode = TRANSIT;
else if (strcmp(argv[1], "daytime") == 0) mode = DAYTIME;
else if (strcmp(argv[1], "nighttime") == 0) mode = NIGHTTIME;
/* skip mode parameter for following getopt() */
argc--;
argv++;
}
@ -160,13 +169,13 @@ int main(int argc, char *argv[]) {
switch (c) {
case 't':
if (strcmp(optarg, "civil") == 0) {
twilight = -6.0;
horizon = LN_SOLAR_CIVIL_HORIZON;
}
else if (strcmp(optarg, "nautic") == 0) {
twilight = -12.0;
horizon = LN_SOLAR_NAUTIC_HORIZON;
}
else if (strcmp(optarg, "astro") == 0) {
twilight = -18.0;
else if (strcmp(optarg, "astronomical") == 0) {
horizon = LN_SOLAR_ASTRONOMICAL_HORIZON;
}
else {
fprintf(stderr, "invalid twilight: %s\n", optarg);
@ -175,10 +184,21 @@ int main(int argc, char *argv[]) {
break;
case 'd':
if (strptime(optarg, "%Y-%m-%d", &date) == 0) {
if (strptime(optarg, "%Y-%m-%d", &date) == NULL) {
fprintf(stderr, "invalid date: %s\n", optarg);
error = true;
}
else {
time_t t = mktime(&date);
julian = ln_get_julian_from_timet(&t);
#ifdef DEBUG
char date_str[64];
strftime(date_str, 64, "%Y-%m-%d", &date);
printf("parsed date: %s\n", date_str);
#endif
}
break;
case 'f':
@ -186,11 +206,11 @@ int main(int argc, char *argv[]) {
break;
case 'a':
pos.lat = strtod(optarg, NULL);
observer.lat = strtod(optarg, NULL);
break;
case 'o':
pos.lon = strtod(optarg, NULL);
observer.lng = strtod(optarg, NULL);
break;
#ifdef GEONAMES_SUPPORT
case 'q':
@ -199,7 +219,7 @@ int main(int argc, char *argv[]) {
#endif
case 'z':
timezone = atoi(optarg);
timezone = -3600 * atoi(optarg);
break;
case 'v':
@ -225,33 +245,19 @@ int main(int argc, char *argv[]) {
#ifdef GEONAMES_SUPPORT
/* lookup place at http://geonames.org */
if (query && geonames_lookup(query, &pos, NULL, 0) != 0) {
if (query && geonames_lookup(query, (struct pos *) &observer, NULL, 0) != 0) {
fprintf(stderr, "failed to lookup location: %s\n", query);
error = true;
}
#endif
/* validate coordinates */
if (pos.lat == INFINITY) {
fprintf(stderr, "latitude is missing\n");
if (fabs(observer.lat) > 90) {
fprintf(stderr, "invalid latitude\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, "longitude is missing\n");
error = true;
}
else if (fabs(pos.lon) > 180) {
fprintf(stderr, "invalid longitude: %.4f\n", pos.lon);
if (fabs(observer.lng) > 180) {
fprintf(stderr, "invalid longitude\n");
error = true;
}
@ -264,48 +270,63 @@ int main(int argc, char *argv[]) {
#ifdef DEBUG
char date_str[64];
strftime(date_str, 64, "%Y-%m-%d", &date);
time_t t;
ln_get_timet_from_julian(julian, &t);
strftime(date_str, 64, "%Y-%m-%d", gmtime(&t));
printf("calculate for: %s\n", date_str);
printf("for position: %f, %f\n", pos.lat, pos.lon);
printf("with twilight: %f\n", twilight);
printf("for position: %f, %f\n", observer.lat, observer.lng);
printf("with horizon: %f\n", horizon);
printf("with timezone: UTC +%dh\n", timezone / -3600);
#endif
/* start the calculation */
struct tm result = sun(pos, date, mode, twilight, timezone);
if (ln_get_solar_rst_horizon(julian, &observer, horizon, &rst) == 1) {
fprintf(stderr, "sun is circumpolar\n");
return EXIT_CIRCUMPOLAR;
}
else {
double result_jd;
char result_str[64];
struct tm result_date;
struct ln_date result_ln;
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
switch (mode) {
case RISE: result_jd = rst.rise; break;
case SET: result_jd = rst.set; break;
case TRANSIT: result_jd = rst.transit; break;
case DAYTIME: result_jd = rst.set - rst.rise; break;
case NIGHTTIME: result_jd = rst.set - rst.rise; break;
}
result.tm_year = -1900;
result.tm_mon = -1;
result.tm_mday = 0;
}
if (mode == DAYTIME || mode == NIGHTTIME) {
ln_get_date(result_jd - 0.5, &result_ln);
strftime(result_str, 64, format, &result);
printf("%s\n", result_str);
if (strstr(format, "%s") != NULL) {
char timestamp_str[16];
int seconds = round(result_jd * 86400);
snprintf(timestamp_str, sizeof(timestamp_str), "%lu", seconds);
format = strreplace(format, "%s", timestamp_str);
}
return EXIT_SUCCESS;
result_date.tm_year = -1900;
result_date.tm_mon = -1;
result_date.tm_mday = 0;
}
else {
ln_get_date(result_jd - timezone / 86400.0, &result_ln);
result_date.tm_year = result_ln.years - 1900;
result_date.tm_mon = result_ln.months - 1;
result_date.tm_mday = result_ln.days;
}
result_date.tm_hour = result_ln.hours;
result_date.tm_min = result_ln.minutes;
result_date.tm_sec = result_ln.seconds;
strftime(result_str, 64, format, &result_date);
printf("%s\n", result_str);
return EXIT_SUCCESS;
}
}