Initial import.

6 files changed, 735 insertions, 0 deletions

diff --git a/SConstruct b/SConstruct
new file mode 100644
index 0000000..17f2140
--- /dev/null
+++ b/SConstruct
@@ -0,0 +1,8 @@
+
+sources = ['redshift.c', 'solar.c', 'colortemp.c']
+
+env = Environment()
+env.ParseConfig('pkg-config --cflags --libs xcb xcb-randr')
+env.Program('redshift', sources,
+            CFLAGS='-std=c99 -D_BSD_SOURCE',
+            LINKFLAGS='-lm')
diff --git a/colortemp.c b/colortemp.c
new file mode 100644
index 0000000..5b40df4
--- /dev/null
+++ b/colortemp.c
@@ -0,0 +1,260 @@
+/* colortemp.c */
+/* Copyright (c) 2009, Jon Lund Steffensen <jonlst@gmail.com> */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+#include <xcb/xcb.h>
+#include <xcb/randr.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.9917, 1.0014, 1.0149,
+	0.9696, 0.9885, 1.0149,
+	0.9488, 0.9761, 1.0149,
+	0.9290, 0.9642, 1.0149,
+	0.9102, 0.9529, 1.0149,
+	0.8923, 0.9420, 1.0149,
+	0.8753, 0.9316, 1.0149,
+	0.8591, 0.9215, 1.0149,
+	0.8437, 0.9120, 1.0149,
+	0.8289, 0.9028, 1.0149,
+	0.8149, 0.8939, 1.0149,
+	0.8014, 0.8854, 1.0149,
+	0.7885, 0.8772, 1.0149,
+	0.7762, 0.8693, 1.0149,
+	0.7644, 0.8617, 1.0149,
+	0.7531, 0.8543, 1.0149,
+	0.7423, 0.8472, 1.0149,
+	0.7319, 0.8404, 1.0149,
+	0.7219, 0.8338, 1.0149,
+	0.7123, 0.8274, 1.0149,
+	0.7030, 0.8213, 1.0149,
+	0.6941, 0.8152, 1.0149,
+	0.6856, 0.8094, 1.0149,
+	0.6773, 0.8039, 1.0149,
+	0.6693, 0.7984, 1.0149,
+	0.6617, 0.7932, 1.0149,
+	0.6543, 0.7881, 1.0149,
+	0.6471, 0.7832, 1.0149,
+	0.6402, 0.7784, 1.0149,
+	0.6335, 0.7737, 1.0149,
+	0.6271, 0.7692, 1.0149,
+	0.6208, 0.7648, 1.0149,
+	0.6148, 0.7605, 1.0149,
+	0.6089, 0.7564, 1.0149,
+	0.6033, 0.7524, 1.0149  /* 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];
+}
+
+int
+colortemp_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
+colortemp_set_temperature(int temp, float gamma)
+{
+	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;
+	}
+
+	/* 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);
+
+	printf("White point: %f, %f, %f\n",
+	       white_point[0], white_point[1], white_point[2]);
+
+	/* 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];
+
+	for (int i = 0; i < gamma_ramp_size; i++) {
+		gamma_r[i] = pow((float)i/gamma_ramp_size, 1.0/gamma) *
+			UINT16_MAX * white_point[0];
+		gamma_g[i] = pow((float)i/gamma_ramp_size, 1.0/gamma) *
+			UINT16_MAX * white_point[1];
+		gamma_b[i] = pow((float)i/gamma_ramp_size, 1.0/gamma) *
+			UINT16_MAX * white_point[2];
+	}
+
+	/* 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/colortemp.h b/colortemp.h
new file mode 100644
index 0000000..d2e9630
--- /dev/null
+++ b/colortemp.h
@@ -0,0 +1,9 @@
+/* colortemp.h */
+
+#ifndef _COLORTEMP_H
+#define _COLORTEMP_H
+
+int colortemp_check_extension();
+int colortemp_set_temperature(int temp, float gamma);
+
+#endif /* ! _COLORTEMP_H */
diff --git a/redshift.c b/redshift.c
new file mode 100644
index 0000000..ef71609
--- /dev/null
+++ b/redshift.c
@@ -0,0 +1,123 @@
+/* redshift.c */
+/* Copyright (c) 2009, Jon Lund Steffensen <jonlst@gmail.com> */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <locale.h>
+
+#include "solar.h"
+#include "colortemp.h"
+
+#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 USAGE  \
+	"Usage: %s LAT LON DAY-TEMP NIGHT-TEMP [GAMMA]\n"
+
+#define DEG_CHAR  0xb0
+
+
+static void
+printtime(time_t time)
+{
+	char s[64];
+	strftime(s, 64, "%a, %d %b %Y %H:%M:%S %z", localtime(&time));
+	printf("%s\n", s);
+}
+
+int
+main(int argc, char *argv[])
+{
+	/* Check extensions needed for color temperature adjustment. */
+	int r = colortemp_check_extension();
+	if (r < 0) {
+		fprintf(stderr, "Unable to set color temperature.\n");
+		exit(EXIT_FAILURE);
+	}
+
+	/* Init locale for special symbols */
+	char *loc = setlocale(LC_CTYPE, "");
+	if (loc == NULL) {
+		fprintf(stderr, "Unable to set locale.\n");
+		exit(EXIT_FAILURE);
+	}
+
+	/* Load arguments */
+	if (argc < 5) {
+		printf(USAGE, argv[0]);
+		exit(EXIT_FAILURE);
+	}
+
+	/* Latitude */
+	float lat = atof(argv[1]);
+	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 */
+	float lon = atof(argv[2]);
+	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 */
+	int temp_day = atoi(argv[3]);
+	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 */
+	int temp_night = atoi(argv[4]);
+	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 value */
+	float gamma = 1.0;
+	if (argc > 5) gamma = atof(argv[5]);
+
+	/* Current angular elevation of the sun */
+	time_t now = time(NULL);
+	double elevation = solar_elevation(now, lat, lon);
+	printf("Solar elevation is %f\n", elevation);
+
+	/* Use elevation of sun to set color temperature */
+	int temp = 0;
+	if (elevation < SOLAR_CIVIL_TWILIGHT_ELEV) {
+		temp = temp_night;
+		printf("Set color temperature to %uK\n", temp_night);
+	} else if (elevation < SOLAR_DAYTIME_ELEV) {
+		float a = (SOLAR_DAYTIME_ELEV - elevation) /
+			(SOLAR_DAYTIME_ELEV - SOLAR_CIVIL_TWILIGHT_ELEV);
+		temp = (1.0-a)*temp_day + a*temp_night;
+		printf("Interpolated (%f) temperature is %uK\n", a, temp);
+	} else {
+		temp = temp_day;
+		printf("Set color temperature to %uK\n", temp_day);
+	}
+
+	/* Set color temperature */
+	r = colortemp_set_temperature(temp, gamma);
+	if (r < 0) {
+		fprintf(stderr, "Unable to set color temperature.\n");
+		exit(EXIT_FAILURE);
+	}
+
+	return EXIT_SUCCESS;
+}
diff --git a/solar.c b/solar.c
new file mode 100644
index 0000000..5a9704b
--- /dev/null
+++ b/solar.c
@@ -0,0 +1,302 @@
+/* solar.c */
+/* 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/solar.h b/solar.h
new file mode 100644
index 0000000..dc6e513
--- /dev/null
+++ b/solar.h
@@ -0,0 +1,33 @@
+/* solar.h */
+
+#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 */