Move source and headers to src dir.

7 files changed, 973 insertions, 0 deletions

diff --git a/src/colorramp.c b/src/colorramp.c
new file mode 100644
index 0000000..2743036
--- /dev/null
+++ b/src/colorramp.c
@@ -0,0 +1,150 @@
+/* colorramp.c -- color temperature calculation source
+   This file is part of Redshift.
+
+   Redshift 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
+   (at your option) any later version.
+
+   Redshift 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 Redshift.  If not, see <http://www.gnu.org/licenses/>.
+
+   Copyright (c) 2009  Jon Lund Steffensen <jonlst@gmail.com>
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+
+/* Source: http://www.vendian.org/mncharity/dir3/blackbody/
+   Rescaled to make exactly 6500K equal to full intensity in all channels. */
+static const float blackbody_color[] = {
+	1.0000, 0.0425, 0.0000, /* 1000K */
+	1.0000, 0.0668, 0.0000, /* 1100K */
+	1.0000, 0.0911, 0.0000, /* 1200K */
+	1.0000, 0.1149, 0.0000, /* ... */
+	1.0000, 0.1380, 0.0000,
+	1.0000, 0.1604, 0.0000,
+	1.0000, 0.1819, 0.0000,
+	1.0000, 0.2024, 0.0000,
+	1.0000, 0.2220, 0.0000,
+	1.0000, 0.2406, 0.0000,
+	1.0000, 0.2630, 0.0062,
+	1.0000, 0.2868, 0.0155,
+	1.0000, 0.3102, 0.0261,
+	1.0000, 0.3334, 0.0379,
+	1.0000, 0.3562, 0.0508,
+	1.0000, 0.3787, 0.0650,
+	1.0000, 0.4008, 0.0802,
+	1.0000, 0.4227, 0.0964,
+	1.0000, 0.4442, 0.1136,
+	1.0000, 0.4652, 0.1316,
+	1.0000, 0.4859, 0.1505,
+	1.0000, 0.5062, 0.1702,
+	1.0000, 0.5262, 0.1907,
+	1.0000, 0.5458, 0.2118,
+	1.0000, 0.5650, 0.2335,
+	1.0000, 0.5839, 0.2558,
+	1.0000, 0.6023, 0.2786,
+	1.0000, 0.6204, 0.3018,
+	1.0000, 0.6382, 0.3255,
+	1.0000, 0.6557, 0.3495,
+	1.0000, 0.6727, 0.3739,
+	1.0000, 0.6894, 0.3986,
+	1.0000, 0.7058, 0.4234,
+	1.0000, 0.7218, 0.4485,
+	1.0000, 0.7375, 0.4738,
+	1.0000, 0.7529, 0.4992,
+	1.0000, 0.7679, 0.5247,
+	1.0000, 0.7826, 0.5503,
+	1.0000, 0.7970, 0.5760,
+	1.0000, 0.8111, 0.6016,
+	1.0000, 0.8250, 0.6272,
+	1.0000, 0.8384, 0.6529,
+	1.0000, 0.8517, 0.6785,
+	1.0000, 0.8647, 0.7040,
+	1.0000, 0.8773, 0.7294,
+	1.0000, 0.8897, 0.7548,
+	1.0000, 0.9019, 0.7801,
+	1.0000, 0.9137, 0.8051,
+	1.0000, 0.9254, 0.8301,
+	1.0000, 0.9367, 0.8550,
+	1.0000, 0.9478, 0.8795,
+	1.0000, 0.9587, 0.9040,
+	1.0000, 0.9694, 0.9283,
+	1.0000, 0.9798, 0.9524,
+	1.0000, 0.9900, 0.9763,
+	1.0000, 1.0000, 1.0000, /* 6500K */
+	0.9771, 0.9867, 1.0000,
+	0.9554, 0.9740, 1.0000,
+	0.9349, 0.9618, 1.0000,
+	0.9154, 0.9500, 1.0000,
+	0.8968, 0.9389, 1.0000,
+	0.8792, 0.9282, 1.0000,
+	0.8624, 0.9179, 1.0000,
+	0.8465, 0.9080, 1.0000,
+	0.8313, 0.8986, 1.0000,
+	0.8167, 0.8895, 1.0000,
+	0.8029, 0.8808, 1.0000,
+	0.7896, 0.8724, 1.0000,
+	0.7769, 0.8643, 1.0000,
+	0.7648, 0.8565, 1.0000,
+	0.7532, 0.8490, 1.0000,
+	0.7420, 0.8418, 1.0000,
+	0.7314, 0.8348, 1.0000,
+	0.7212, 0.8281, 1.0000,
+	0.7113, 0.8216, 1.0000,
+	0.7018, 0.8153, 1.0000,
+	0.6927, 0.8092, 1.0000,
+	0.6839, 0.8032, 1.0000,
+	0.6755, 0.7975, 1.0000,
+	0.6674, 0.7921, 1.0000,
+	0.6595, 0.7867, 1.0000,
+	0.6520, 0.7816, 1.0000,
+	0.6447, 0.7765, 1.0000,
+	0.6376, 0.7717, 1.0000,
+	0.6308, 0.7670, 1.0000,
+	0.6242, 0.7623, 1.0000,
+	0.6179, 0.7579, 1.0000,
+	0.6117, 0.7536, 1.0000,
+	0.6058, 0.7493, 1.0000,
+	0.6000, 0.7453, 1.0000,
+	0.5944, 0.7414, 1.0000 /* 10000K */
+};
+
+
+static void
+interpolate_color(float a, const float *c1, const float *c2, float *c)
+{
+	c[0] = (1.0-a)*c1[0] + a*c2[0];
+	c[1] = (1.0-a)*c1[1] + a*c2[1];
+	c[2] = (1.0-a)*c1[2] + a*c2[2];
+}
+
+void
+colorramp_fill(uint16_t *gamma_r, uint16_t *gamma_g, uint16_t *gamma_b,
+	       int size, int temp, float gamma[3])
+{
+	/* Calculate white point */
+	float white_point[3];
+	float alpha = (temp % 100) / 100.0;
+	int temp_index = ((temp - 1000) / 100)*3;
+	interpolate_color(alpha, &blackbody_color[temp_index],
+			  &blackbody_color[temp_index+3], white_point);
+
+	for (int i = 0; i < size; i++) {
+		gamma_r[i] = pow((float)i/size, 1.0/gamma[0]) *
+			UINT16_MAX * white_point[0];
+		gamma_g[i] = pow((float)i/size, 1.0/gamma[1]) *
+			UINT16_MAX * white_point[1];
+		gamma_b[i] = pow((float)i/size, 1.0/gamma[2]) *
+			UINT16_MAX * white_point[2];
+	}
+}
diff --git a/src/colorramp.h b/src/colorramp.h
new file mode 100644
index 0000000..e7f5554
--- /dev/null
+++ b/src/colorramp.h
@@ -0,0 +1,28 @@
+/* colorramp.h -- color temperature calculation header
+   This file is part of Redshift.
+
+   Redshift 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
+   (at your option) any later version.
+
+   Redshift 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 Redshift.  If not, see <http://www.gnu.org/licenses/>.
+
+   Copyright (c) 2009  Jon Lund Steffensen <jonlst@gmail.com>
+*/
+
+#ifndef _REDSHIFT_COLORRAMP_H
+#define _REDSHIFT_COLORRAMP_H
+
+#include <stdint.h>
+
+void colorramp_fill(uint16_t *gamma_r, uint16_t *gamma_g, uint16_t *gamma_b,
+		    int size, int temp, float gamma[3]);
+
+#endif /* ! _REDSHIFT_COLORRAMP_H */
diff --git a/src/randr.c b/src/randr.c
new file mode 100644
index 0000000..5226b81
--- /dev/null
+++ b/src/randr.c
@@ -0,0 +1,155 @@
+/* randr.c -- X RandR gamma adjustment source
+   This file is part of Redshift.
+
+   Redshift 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
+   (at your option) any later version.
+
+   Redshift 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 Redshift.  If not, see <http://www.gnu.org/licenses/>.
+
+   Copyright (c) 2009  Jon Lund Steffensen <jonlst@gmail.com>
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <xcb/xcb.h>
+#include <xcb/randr.h>
+
+#include "colorramp.h"
+
+
+int
+randr_check_extension()
+{
+	xcb_generic_error_t *error;
+
+	/* Open X server connection */
+	xcb_connection_t *conn = xcb_connect(NULL, NULL);
+
+	/* Query RandR version */
+	xcb_randr_query_version_cookie_t ver_cookie =
+		xcb_randr_query_version(conn, 1, 3);
+	xcb_randr_query_version_reply_t *ver_reply =
+		xcb_randr_query_version_reply(conn, ver_cookie, &error);
+
+	if (error) {
+		fprintf(stderr, "RANDR Query Version, error: %d\n",
+			error->error_code);
+		xcb_disconnect(conn);
+		return -1;
+	}
+
+	if (ver_reply->major_version < 1 || ver_reply->minor_version < 3) {
+		fprintf(stderr, "Unsupported RANDR version (%u.%u)\n",
+			ver_reply->major_version, ver_reply->minor_version);
+		free(ver_reply);
+		xcb_disconnect(conn);
+		return -1;
+	}
+
+	free(ver_reply);
+
+	/* Close connection */
+	xcb_disconnect(conn);
+
+	return 0;
+}
+
+int
+randr_set_temperature(int temp, float gamma[3])
+{
+	xcb_generic_error_t *error;
+
+	/* Open X server connection */
+	xcb_connection_t *conn = xcb_connect(NULL, NULL);
+
+	/* Get first screen */
+	const xcb_setup_t *setup = xcb_get_setup(conn);
+	xcb_screen_iterator_t iter = xcb_setup_roots_iterator(setup);
+	xcb_screen_t *screen = iter.data;
+
+	/* Get list of CRTCs for the screen */
+	xcb_randr_get_screen_resources_current_cookie_t res_cookie =
+		xcb_randr_get_screen_resources_current(conn, screen->root);
+	xcb_randr_get_screen_resources_current_reply_t *res_reply =
+		xcb_randr_get_screen_resources_current_reply(conn, res_cookie,
+							     &error);
+
+	if (error) {
+		fprintf(stderr, "RANDR Get Screen Resources Current,"
+			" error: %d\n", error->error_code);
+		xcb_disconnect(conn);
+		return -1;
+	}
+
+	xcb_randr_crtc_t *crtcs =
+		xcb_randr_get_screen_resources_current_crtcs(res_reply);
+	xcb_randr_crtc_t crtc = crtcs[0];
+
+	free(res_reply);
+
+	/* Request size of gamma ramps */
+	xcb_randr_get_crtc_gamma_size_cookie_t gamma_size_cookie =
+		xcb_randr_get_crtc_gamma_size(conn, crtc);
+	xcb_randr_get_crtc_gamma_size_reply_t *gamma_size_reply =
+		xcb_randr_get_crtc_gamma_size_reply(conn, gamma_size_cookie,
+						    &error);
+
+	if (error) {
+		fprintf(stderr, "RANDR Get CRTC Gamma Size, error: %d\n",
+			error->error_code);
+		xcb_disconnect(conn);
+		return -1;
+	}
+
+	int gamma_ramp_size = gamma_size_reply->size;
+
+	free(gamma_size_reply);
+
+	if (gamma_ramp_size == 0) {
+		fprintf(stderr, "Error: Gamma ramp size too small, %i\n",
+			gamma_ramp_size);
+		xcb_disconnect(conn);
+		return -1;
+	}
+
+	/* Create new gamma ramps */
+	uint16_t *gamma_ramps = malloc(3*gamma_ramp_size*sizeof(uint16_t));
+	if (gamma_ramps == NULL) abort();
+
+	uint16_t *gamma_r = &gamma_ramps[0*gamma_ramp_size];
+	uint16_t *gamma_g = &gamma_ramps[1*gamma_ramp_size];
+	uint16_t *gamma_b = &gamma_ramps[2*gamma_ramp_size];
+
+	colorramp_fill(gamma_r, gamma_g, gamma_b, gamma_ramp_size,
+		       temp, gamma);
+
+	/* Set new gamma ramps */
+	xcb_void_cookie_t gamma_set_cookie =
+		xcb_randr_set_crtc_gamma_checked(conn, crtc, gamma_ramp_size,
+						 gamma_r, gamma_g, gamma_b);
+	error = xcb_request_check(conn, gamma_set_cookie);
+
+	if (error) {
+		fprintf(stderr, "RANDR Set CRTC Gamma, error: %d\n",
+			error->error_code);
+		free(gamma_ramps);
+		xcb_disconnect(conn);
+		return -1;
+	}
+
+	free(gamma_ramps);
+
+	/* Close connection */
+	xcb_disconnect(conn);
+
+	return 0;
+}
diff --git a/src/randr.h b/src/randr.h
new file mode 100644
index 0000000..7efa94c
--- /dev/null
+++ b/src/randr.h
@@ -0,0 +1,26 @@
+/* randr.h -- X RandR gamma adjustment header
+   This file is part of Redshift.
+
+   Redshift 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
+   (at your option) any later version.
+
+   Redshift 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 Redshift.  If not, see <http://www.gnu.org/licenses/>.
+
+   Copyright (c) 2009  Jon Lund Steffensen <jonlst@gmail.com>
+*/
+
+#ifndef _REDSHIFT_RANDR_H
+#define _REDSHIFT_RANDR_H
+
+int randr_check_extension();
+int randr_set_temperature(int temp, float gamma[3]);
+
+#endif /* ! _REDSHIFT_RANDR_H */
diff --git a/src/redshift.c b/src/redshift.c
new file mode 100644
index 0000000..70bc4fa
--- /dev/null
+++ b/src/redshift.c
@@ -0,0 +1,246 @@
+/* redshift.c -- Main program source
+   This file is part of Redshift.
+
+   Redshift 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
+   (at your option) any later version.
+
+   Redshift 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 Redshift.  If not, see <http://www.gnu.org/licenses/>.
+
+   Copyright (c) 2009  Jon Lund Steffensen <jonlst@gmail.com>
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+#include <locale.h>
+
+#include "solar.h"
+#include "randr.h"
+
+
+/* Bounds for parameters. */
+#define MIN_LAT   -90.0
+#define MAX_LAT    90.0
+#define MIN_LON  -180.0
+#define MAX_LON   180.0
+#define MIN_TEMP   1000
+#define MAX_TEMP  10000
+#define MIN_GAMMA   0.1
+#define MAX_GAMMA  10.0
+
+/* Default values for parameters. */
+#define DEFAULT_DAY_TEMP    5500
+#define DEFAULT_NIGHT_TEMP  3700
+#define DEFAULT_GAMMA        1.0
+
+/* Angular elevation of the sun at which the color temperature
+   transition period starts and ends (in degress).
+   Transition during twilight, and while the sun is lower than
+   3.0 degrees above the horizon. */
+#define TRANSITION_LOW     SOLAR_CIVIL_TWILIGHT_ELEV
+#define TRANSITION_HIGH    3.0
+
+
+#define USAGE  \
+	"Usage: %s -l LAT:LON -t DAY:NIGHT [OPTIONS...]\n"
+#define HELP \
+	USAGE \
+	" Set color temperature of display according to time of day.\n" \
+	"  -g R:G:B\tAdditional gamma correction to apply\n" \
+	"  -h\t\tDisplay this help message\n" \
+	"  -l LAT:LON\tYour current location\n" \
+	"  -t DAY:NIGHT\tColor temperature to set at daytime/night\n" \
+	"  -v\t\tVerbose output\n"
+
+/* DEGREE SIGN is Unicode U+00b0 */
+#define DEG_CHAR  0xb0
+
+
+int
+main(int argc, char *argv[])
+{
+	/* Check extensions needed for color temperature adjustment. */
+	int r = randr_check_extension();
+	if (r < 0) {
+		fprintf(stderr, "Unable to set color temperature:"
+			" Needed extension is missing.\n");
+		exit(EXIT_FAILURE);
+	}
+
+	/* Init locale for degree symbol. */
+	char *loc = setlocale(LC_CTYPE, "");
+	if (loc == NULL) {
+		fprintf(stderr, "Unable to set locale.\n");
+		exit(EXIT_FAILURE);
+	}
+
+	/* Parse arguments. */
+	float lat = NAN;
+	float lon = NAN;
+	int temp_day = DEFAULT_DAY_TEMP;
+	int temp_night = DEFAULT_NIGHT_TEMP;
+	float gamma[3] = { DEFAULT_GAMMA, DEFAULT_GAMMA, DEFAULT_GAMMA };
+	int verbose = 0;
+	char *s;
+
+	int opt;
+	while ((opt = getopt(argc, argv, "g:hl:t:v")) != -1) {
+		switch (opt) {
+		case 'g':
+			s = strchr(optarg, ':');
+			if (s == NULL) {
+				/* Use value for all channels */
+				float g = atof(optarg);
+				gamma[0] = gamma[1] = gamma[2] = g;
+			} else {
+				/* Parse separate value for each channel */
+				*(s++) = '\0';
+				gamma[0] = atof(optarg); /* Red */
+				char *g_s = s;
+				s = strchr(s, ':');
+				if (s == NULL) {
+					fprintf(stderr, USAGE, argv[0]);
+					exit(EXIT_FAILURE);
+				}
+
+				*(s++) = '\0';
+				gamma[1] = atof(g_s); /* Blue */
+				gamma[2] = atof(s); /* Green */
+			}
+			break;
+		case 'h':
+			printf(HELP, argv[0]);
+			exit(EXIT_SUCCESS);
+			break;
+		case 'l':
+			s = strchr(optarg, ':');
+			if (s == NULL) {
+				fprintf(stderr, USAGE, argv[0]);
+				exit(EXIT_FAILURE);
+			}
+			*(s++) = '\0';
+			lat = atof(optarg);
+			lon = atof(s);
+			break;
+		case 't':
+			s = strchr(optarg, ':');
+			if (s == NULL) {
+				fprintf(stderr, USAGE, argv[0]);
+				exit(EXIT_FAILURE);
+			}
+			*(s++) = '\0';
+			temp_day = atoi(optarg);
+			temp_night = atoi(s);
+			break;
+		case 'v':
+			verbose = 1;
+			break;
+		default:
+			fprintf(stderr, USAGE, argv[0]);
+			exit(EXIT_FAILURE);
+		}
+	}
+
+	/* Latitude and longitude must be set */
+	if (isnan(lat) || isnan(lon)) {
+		fprintf(stderr, USAGE, argv[0]);
+		fprintf(stderr, "Latitude and longitude must be set.\n");
+		exit(EXIT_FAILURE);
+	}
+
+	if (verbose) {
+		printf("Location: %f%lc, %f%lc\n",
+		       lat, DEG_CHAR, lon, DEG_CHAR);
+	}
+
+	/* Latitude */
+	if (lat < MIN_LAT || lat > MAX_LAT) {
+		fprintf(stderr,
+			"Latitude must be between %.1f%lc and %.1f%lc.\n",
+			MIN_LAT, DEG_CHAR, MAX_LAT, DEG_CHAR);
+		exit(EXIT_FAILURE);
+	}
+
+	/* Longitude */
+	if (lon < MIN_LON || lon > MAX_LON) {
+		fprintf(stderr,
+			"Longitude must be between %.1f%lc and %.1f%lc.\n",
+			MIN_LON, DEG_CHAR, MAX_LON, DEG_CHAR);
+		exit(EXIT_FAILURE);
+	}
+
+	/* Color temperature at daytime */
+	if (temp_day < MIN_TEMP || temp_day >= MAX_TEMP) {
+		fprintf(stderr, "Temperature must be between %uK and %uK.\n",
+			MIN_TEMP, MAX_TEMP);
+		exit(EXIT_FAILURE);
+	}
+
+	/* Color temperature at night */
+	if (temp_night < MIN_TEMP || temp_night >= MAX_TEMP) {
+		fprintf(stderr, "Temperature must be between %uK and %uK.\n",
+			MIN_TEMP, MAX_TEMP);
+		exit(EXIT_FAILURE);
+	}
+
+	/* Gamma */
+	if (gamma[0] < MIN_GAMMA || gamma[0] > MAX_GAMMA ||
+	    gamma[1] < MIN_GAMMA || gamma[1] > MAX_GAMMA ||
+	    gamma[2] < MIN_GAMMA || gamma[2] > MAX_GAMMA) {
+		fprintf(stderr, "Gamma value must be between %.1f and %.1f.\n",
+			MIN_GAMMA, MAX_GAMMA);
+		exit(EXIT_FAILURE);
+	}
+
+	/* Current angular elevation of the sun */
+	time_t now = time(NULL);
+	double elevation = solar_elevation(now, lat, lon);
+
+	if (verbose) {
+		printf("Solar elevation: %f%lc\n", elevation, DEG_CHAR);
+	}
+
+	/* Use elevation of sun to set color temperature */
+	int temp = 0;
+	if (elevation < TRANSITION_LOW) {
+		temp = temp_night;
+		if (verbose) printf("Period: Night\n");
+	} else if (elevation < TRANSITION_HIGH) {
+		/* Transition period: interpolate */
+		float a = (TRANSITION_LOW - elevation) /
+			(TRANSITION_LOW - TRANSITION_HIGH);
+		temp = (1.0-a)*temp_night + a*temp_day;
+		if (verbose) {
+			printf("Period: Transition (%.2f%% day)\n", a*100);
+		}
+	} else {
+		temp = temp_day;
+		if (verbose) printf("Period: Daytime\n");
+	}
+
+	if (verbose) {
+		printf("Color temperature: %uK\n", temp);
+		printf("Gamma: %.3f, %.3f, %.3f\n",
+		       gamma[0], gamma[1], gamma[2]);
+	}
+
+	/* Set color temperature */
+	r = randr_set_temperature(temp, gamma);
+	if (r < 0) {
+		fprintf(stderr, "Unable to set color temperature.\n");
+		exit(EXIT_FAILURE);
+	}
+
+	return EXIT_SUCCESS;
+}
diff --git a/src/solar.c b/src/solar.c
new file mode 100644
index 0000000..2d66fff
--- /dev/null
+++ b/src/solar.c
@@ -0,0 +1,318 @@
+/* solar.c -- Solar position source
+   This file is part of Redshift.
+
+   Redshift 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
+   (at your option) any later version.
+
+   Redshift 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 Redshift.  If not, see <http://www.gnu.org/licenses/>.
+
+   Copyright (c) 2009  Jon Lund Steffensen <jonlst@gmail.com>
+*/
+
+/* Ported from javascript code by U.S. Department of Commerce,
+   National Oceanic & Atmospheric Administration:
+   http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
+   It is based on equations from "Astronomical Algorithms" by
+   Jean Meeus. */
+
+#include <math.h>
+#include <time.h>
+
+#include "solar.h"
+
+#define RAD(x)  ((x)*(M_PI/180))
+#define DEG(x)  ((x)*(180/M_PI))
+
+
+/* Angels of various times of day. */
+static const double time_angle[] = {
+	[SOLAR_TIME_ASTRO_DAWN] = RAD(-90.0 + SOLAR_ASTRO_TWILIGHT_ELEV),
+	[SOLAR_TIME_NAUT_DAWN] = RAD(-90.0 + SOLAR_NAUT_TWILIGHT_ELEV),
+	[SOLAR_TIME_CIVIL_DAWN] = RAD(-90.0 + SOLAR_CIVIL_TWILIGHT_ELEV),
+	[SOLAR_TIME_SUNRISE] = RAD(-90.0 + SOLAR_DAYTIME_ELEV),
+	[SOLAR_TIME_NOON] = RAD(0.0),
+	[SOLAR_TIME_SUNSET] = RAD(90.0 - SOLAR_DAYTIME_ELEV),
+	[SOLAR_TIME_CIVIL_DUSK] = RAD(90.0 - SOLAR_CIVIL_TWILIGHT_ELEV),
+	[SOLAR_TIME_NAUT_DUSK] = RAD(90.0 - SOLAR_NAUT_TWILIGHT_ELEV),
+	[SOLAR_TIME_ASTRO_DUSK] = RAD(90.0 - SOLAR_ASTRO_TWILIGHT_ELEV)
+};
+
+
+/* Unix time from Julian day */
+static time_t
+unix_time_from_jd(double jd)
+{
+	return 86400.0*(jd - 2440587.5);
+}
+
+/* Julian day from unix time */
+static double
+jd_from_unix_time(time_t t)
+{
+	return (t / 86400.0) + 2440587.5;
+}
+
+/* Julian centuries since J2000.0 from Julian day */
+static double
+jcent_from_jd(double jd)
+{
+	return (jd - 2451545.0) / 36525.0;
+}
+
+/* Julian day from Julian centuries since J2000.0 */
+static double
+jd_from_jcent(double t)
+{
+	return 36525.0*t + 2451545.0;
+}
+
+/* Geometric mean longitude of the sun.
+   t: Julian centuries since J2000.0
+   Return: Geometric mean logitude in radians. */
+static double
+sun_geom_mean_lon(double t)
+{
+	/* FIXME returned value should always be positive */
+	return RAD(fmod(280.46646 + t*(36000.76983 + t*0.0003032), 360));
+}
+
+/* Geometric mean anomaly of the sun.
+   t: Julian centuries since J2000.0
+   Return: Geometric mean anomaly in radians. */
+static double
+sun_geom_mean_anomaly(double t)
+{
+	return RAD(357.52911 + t*(35999.05029 - t*0.0001537));
+}
+
+/* Eccentricity of earth orbit.
+   t: Julian centuries since J2000.0
+   Return: Eccentricity (unitless). */
+static double
+earth_orbit_eccentricity(double t)
+{
+	return 0.016708634 - t*(0.000042037 + t*0.0000001267);
+}
+
+/* Equation of center of the sun.
+   t: Julian centuries since J2000.0
+   Return: Center(?) in radians */
+static double
+sun_equation_of_center(double t)
+{
+	/* Use the first three terms of the equation. */
+	double m = sun_geom_mean_anomaly(t);
+	double c = sin(m)*(1.914602 - t*(0.004817 + 0.000014*t)) +
+		sin(2*m)*(0.019993 - 0.000101*t) +
+		sin(3*m)*0.000289;
+	return RAD(c);
+}
+
+/* True longitude of the sun.
+   t: Julian centuries since J2000.0
+   Return: True longitude in radians */
+static double
+sun_true_lon(double t)
+{
+	double l_0 = sun_geom_mean_lon(t);
+	double c = sun_equation_of_center(t);
+	return l_0 + c;
+}
+
+/* Apparent longitude of the sun. (Right ascension).
+   t: Julian centuries since J2000.0
+   Return: Apparent longitude in radians */
+static double
+sun_apparent_lon(double t)
+{
+	double o = sun_true_lon(t);
+	return RAD(DEG(o) - 0.00569 - 0.00478*sin(RAD(125.04 - 1934.136*t)));
+}
+
+/* Mean obliquity of the ecliptic
+   t: Julian centuries since J2000.0
+   Return: Mean obliquity in radians */
+static double
+mean_ecliptic_obliquity(double t)
+{
+	double sec = 21.448 - t*(46.815 + t*(0.00059 - t*0.001813));
+	return RAD(23.0 + (26.0 + (sec/60.0))/60.0);
+}
+
+/* Corrected obliquity of the ecliptic.
+   t: Julian centuries since J2000.0
+   Return: Currected obliquity in radians */
+static double
+obliquity_corr(double t)
+{
+	double e_0 = mean_ecliptic_obliquity(t);
+	double omega = 125.04 - t*1934.136;
+	return RAD(DEG(e_0) + 0.00256*cos(RAD(omega)));
+}
+
+/* Declination of the sun.
+   t: Julian centuries since J2000.0
+   Return: Declination in radians */
+static double
+solar_declination(double t)
+{
+	double e = obliquity_corr(t);
+	double lambda = sun_apparent_lon(t);
+	return asin(sin(e)*sin(lambda));
+}
+
+/* Difference between true solar time and mean solar time.
+   t: Julian centuries since J2000.0
+   Return: Difference in minutes */
+static double
+equation_of_time(double t)
+{
+	double epsilon = obliquity_corr(t);
+	double l_0 = sun_geom_mean_lon(t);
+	double e = earth_orbit_eccentricity(t);
+	double m = sun_geom_mean_anomaly(t);
+	double y = pow(tan(epsilon/2.0), 2.0);
+
+	double eq_time = y*sin(2*l_0) - 2*e*sin(m) +
+		4*e*y*sin(m)*cos(2*l_0) -
+		0.5*y*y*sin(4*l_0) -
+		1.25*e*e*sin(2*m);
+	return 4*DEG(eq_time);
+}
+
+/* Hour angle at the location for the given angular elevation.
+   lat: Latitude of location in degrees
+   decl: Declination in radians
+   elev: Angular elevation angle in radians
+   Return: Hour angle in radians */
+static double
+hour_angle_from_elevation(double lat, double decl, double elev)
+{
+	double omega = acos((cos(fabs(elev)) - sin(RAD(lat))*sin(decl))/
+			    (cos(RAD(lat))*cos(decl)));
+	return copysign(omega, -elev);
+}
+
+/* Angular elevation at the location for the given hour angle.
+   lat: Latitude of location in degrees
+   decl: Declination in radians
+   ha: Hour angle in radians
+   Return: Angular elevation in radians */
+static double
+elevation_from_hour_angle(double lat, double decl, double ha)
+{
+	return asin(cos(ha)*cos(RAD(lat))*cos(decl) +
+		    sin(RAD(lat))*sin(decl));
+}
+
+/* Time of apparent solar noon of location on earth.
+   t: Julian centuries since J2000.0
+   lon: Longitude of location in degrees
+   Return: Time difference from mean solar midnigth in minutes */
+static double
+time_of_solar_noon(double t, double lon)
+{
+	/* First pass uses approximate solar noon to
+	   calculate equation of time. */
+	double t_noon = jcent_from_jd(jd_from_jcent(t) - lon/360.0);
+	double eq_time = equation_of_time(t_noon);
+	double sol_noon = 720 - 4*lon - eq_time;
+
+	/* Recalculate using new solar noon. */
+	t_noon = jcent_from_jd(jd_from_jcent(t) - 0.5 + sol_noon/1440.0);
+	eq_time = equation_of_time(t_noon);
+	sol_noon = 720 - 4*lon - eq_time;
+
+	/* No need to do more iterations */
+	return sol_noon;
+}
+
+/* Time of given apparent solar angular elevation of location on earth.
+   t: Julian centuries since J2000.0
+   t_noon: Apparent solar noon in Julian centuries since J2000.0
+   lat: Latitude of location in degrees
+   lon: Longtitude of location in degrees
+   elev: Solar angular elevation in radians
+   Return: Time difference from mean solar midnight in minutes */
+static double
+time_of_solar_elevation(double t, double t_noon,
+			double lat, double lon, double elev)
+{
+	/* First pass uses approximate sunrise to
+	   calculate equation of time. */
+	double eq_time = equation_of_time(t_noon);
+	double sol_decl = solar_declination(t_noon);
+	double ha = hour_angle_from_elevation(lat, sol_decl, elev);
+	double sol_offset = 720 - 4*(lon + DEG(ha)) - eq_time;
+
+	/* Recalculate using new sunrise. */
+	double t_rise = jcent_from_jd(jd_from_jcent(t) + sol_offset/1440.0);
+	eq_time = equation_of_time(t_rise);
+	sol_decl = solar_declination(t_rise);
+	ha = hour_angle_from_elevation(lat, sol_decl, elev);
+	sol_offset = 720 - 4*(lon + DEG(ha)) - eq_time;
+
+	/* No need to do more iterations */
+	return sol_offset;
+}
+
+/* Solar angular elevation at the given location and time.
+   t: Julian centuries since J2000.0
+   lat: Latitude of location
+   lon: Longitude of location
+   Return: Solar angular elevation in radians */
+static double
+solar_elevation_from_time(double t, double lat, double lon)
+{
+	/* Minutes from midnight */
+	double jd = jd_from_jcent(t);
+	double offset = (jd - round(jd) - 0.5)*1440.0;
+
+	double eq_time = equation_of_time(t);
+	double ha = RAD((720 - offset - eq_time)/4 - lon);
+	double decl = solar_declination(t);
+	return elevation_from_hour_angle(lat, decl, ha);
+}
+
+double
+solar_elevation(time_t date, double lat, double lon)
+{
+	double jd = jd_from_unix_time(date);
+	return DEG(solar_elevation_from_time(jcent_from_jd(jd), lat, lon));
+}
+
+void
+solar_table_fill(time_t date, double lat, double lon, time_t *table)
+{
+	/* Calculate Julian day */
+	double jd = jd_from_unix_time(date);
+
+	/* Calculate Julian day number */
+	double jdn = round(jd);
+	double t = jcent_from_jd(jdn);
+
+	/* Calculate apparent solar noon */
+	double sol_noon = time_of_solar_noon(t, lon);
+	double j_noon = jdn - 0.5 + sol_noon/1440.0;
+	double t_noon = jcent_from_jd(j_noon);
+	table[SOLAR_TIME_NOON] = unix_time_from_jd(j_noon);
+
+	/* Calculate solar midnight */
+	table[SOLAR_TIME_MIDNIGHT] = unix_time_from_jd(j_noon + 0.5);
+
+	/* Calulate absoute time of other phenomena */
+	for (int i = 2; i < SOLAR_TIME_MAX; i++) {
+		double angle = time_angle[i];
+		double offset =
+			time_of_solar_elevation(t, t_noon, lat, lon, angle);
+		table[i] = unix_time_from_jd(jdn - 0.5 + offset/1440.0);
+	}
+}
diff --git a/src/solar.h b/src/solar.h
new file mode 100644
index 0000000..5897373
--- /dev/null
+++ b/src/solar.h
@@ -0,0 +1,50 @@
+/* solar.h -- Solar position header
+   This file is part of Redshift.
+
+   Redshift 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
+   (at your option) any later version.
+
+   Redshift 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 Redshift.  If not, see <http://www.gnu.org/licenses/>.
+
+   Copyright (c) 2009  Jon Lund Steffensen <jonlst@gmail.com>
+*/
+
+#ifndef _SOLAR_H
+#define _SOLAR_H
+
+#include <time.h>
+
+/* Model of atmospheric refraction near horizon (in degrees). */
+#define SOLAR_ATM_REFRAC  0.833
+
+#define SOLAR_ASTRO_TWILIGHT_ELEV   -18.0
+#define SOLAR_NAUT_TWILIGHT_ELEV    -12.0
+#define SOLAR_CIVIL_TWILIGHT_ELEV    -6.0
+#define SOLAR_DAYTIME_ELEV           (0.0 - SOLAR_ATM_REFRAC)
+
+typedef enum {
+	SOLAR_TIME_NOON = 0,
+	SOLAR_TIME_MIDNIGHT,
+	SOLAR_TIME_ASTRO_DAWN,
+	SOLAR_TIME_NAUT_DAWN,
+	SOLAR_TIME_CIVIL_DAWN,
+	SOLAR_TIME_SUNRISE,
+	SOLAR_TIME_SUNSET,
+	SOLAR_TIME_CIVIL_DUSK,
+	SOLAR_TIME_NAUT_DUSK,
+	SOLAR_TIME_ASTRO_DUSK,
+	SOLAR_TIME_MAX
+} solar_time_t;
+
+double solar_elevation(time_t t, double lat, double lon);
+void solar_table_fill(time_t date, double lat, double lon, time_t *table);
+
+#endif /* ! _SOLAR_H */