split out from radharc

Signed-off-by: Mattias Andrée <maandree@member.fsf.org>

6 files changed, 1087 insertions, 0 deletions

diff --git a/src/10deg b/src/10deg
new file mode 100644
index 0000000..d054a4c
--- /dev/null
+++ b/src/10deg
@@ -0,0 +1,391 @@
+0.6472 0.3506
+0.6348 0.3612
+0.6225 0.3710
+0.6104 0.3797
+0.5983 0.3874
+0.5864 0.3940
+0.5747 0.3996
+0.5632 0.4041
+0.5519 0.4077
+0.5408 0.4104
+0.5300 0.4123
+0.5195 0.4134
+0.5093 0.4139
+0.4995 0.4138
+0.4900 0.4132
+0.4808 0.4122
+0.4720 0.4108
+0.4636 0.4091
+0.4555 0.4071
+0.4478 0.4050
+0.4404 0.4026
+0.4333 0.4002
+0.4266 0.3976
+0.4202 0.3950
+0.4140 0.3923
+0.4082 0.3895
+0.4026 0.3868
+0.3973 0.3841
+0.3922 0.3813
+0.3873 0.3786
+0.3827 0.3759
+0.3783 0.3733
+0.3741 0.3707
+0.3701 0.3681
+0.3662 0.3656
+0.3625 0.3631
+0.3590 0.3607
+0.3557 0.3583
+0.3524 0.3560
+0.3494 0.3538
+0.3464 0.3516
+0.3436 0.3494
+0.3409 0.3473
+0.3383 0.3453
+0.3358 0.3433
+0.3334 0.3413
+0.3311 0.3394
+0.3289 0.3376
+0.3268 0.3358
+0.3247 0.3341
+0.3228 0.3324
+0.3209 0.3307
+0.3190 0.3291
+0.3173 0.3275
+0.3156 0.3260
+0.3140 0.3245
+0.3124 0.3231
+0.3109 0.3217
+0.3094 0.3203
+0.3080 0.3190
+0.3066 0.3177
+0.3053 0.3164
+0.3040 0.3152
+0.3028 0.3140
+0.3016 0.3128
+0.3004 0.3117
+0.2993 0.3106
+0.2982 0.3095
+0.2971 0.3084
+0.2961 0.3074
+0.2951 0.3064
+0.2942 0.3054
+0.2932 0.3044
+0.2923 0.3035
+0.2914 0.3026
+0.2906 0.3017
+0.2897 0.3008
+0.2889 0.3000
+0.2881 0.2992
+0.2874 0.2983
+0.2866 0.2976
+0.2859 0.2968
+0.2852 0.2960
+0.2845 0.2953
+0.2839 0.2946
+0.2832 0.2939
+0.2826 0.2932
+0.2820 0.2925
+0.2813 0.2918
+0.2808 0.2912
+0.2802 0.2905
+0.2796 0.2899
+0.2791 0.2893
+0.2785 0.2887
+0.2780 0.2881
+0.2775 0.2876
+0.2770 0.2870
+0.2765 0.2865
+0.2761 0.2859
+0.2756 0.2854
+0.2751 0.2849
+0.2747 0.2844
+0.2743 0.2839
+0.2738 0.2834
+0.2734 0.2829
+0.2730 0.2825
+0.2726 0.2820
+0.2722 0.2816
+0.2719 0.2811
+0.2715 0.2807
+0.2711 0.2803
+0.2708 0.2799
+0.2704 0.2794
+0.2701 0.2790
+0.2697 0.2786
+0.2694 0.2783
+0.2691 0.2779
+0.2688 0.2775
+0.2684 0.2771
+0.2681 0.2768
+0.2678 0.2764
+0.2675 0.2761
+0.2673 0.2757
+0.2670 0.2754
+0.2667 0.2751
+0.2664 0.2747
+0.2662 0.2744
+0.2659 0.2741
+0.2656 0.2738
+0.2654 0.2735
+0.2651 0.2732
+0.2649 0.2729
+0.2646 0.2726
+0.2644 0.2723
+0.2642 0.2720
+0.2639 0.2718
+0.2637 0.2715
+0.2635 0.2712
+0.2633 0.2709
+0.2631 0.2707
+0.2629 0.2704
+0.2626 0.2702
+0.2624 0.2699
+0.2622 0.2697
+0.2620 0.2694
+0.2618 0.2692
+0.2617 0.2690
+0.2615 0.2687
+0.2613 0.2685
+0.2611 0.2683
+0.2609 0.2681
+0.2607 0.2678
+0.2606 0.2676
+0.2604 0.2674
+0.2602 0.2672
+0.2601 0.2670
+0.2599 0.2668
+0.2597 0.2666
+0.2596 0.2664
+0.2594 0.2662
+0.2593 0.2660
+0.2591 0.2658
+0.2590 0.2656
+0.2588 0.2654
+0.2587 0.2653
+0.2585 0.2651
+0.2584 0.2649
+0.2582 0.2647
+0.2581 0.2645
+0.2580 0.2644
+0.2578 0.2642
+0.2577 0.2640
+0.2575 0.2639
+0.2574 0.2637
+0.2573 0.2635
+0.2572 0.2634
+0.2570 0.2632
+0.2569 0.2631
+0.2568 0.2629
+0.2567 0.2628
+0.2566 0.2626
+0.2564 0.2625
+0.2563 0.2623
+0.2562 0.2622
+0.2561 0.2620
+0.2560 0.2619
+0.2559 0.2618
+0.2558 0.2616
+0.2557 0.2615
+0.2555 0.2613
+0.2554 0.2612
+0.2553 0.2611
+0.2552 0.2609
+0.2551 0.2608
+0.2550 0.2607
+0.2549 0.2606
+0.2548 0.2604
+0.2547 0.2603
+0.2546 0.2602
+0.2546 0.2601
+0.2545 0.2600
+0.2544 0.2598
+0.2543 0.2597
+0.2542 0.2596
+0.2541 0.2595
+0.2540 0.2594
+0.2539 0.2593
+0.2538 0.2591
+0.2537 0.2590
+0.2537 0.2589
+0.2536 0.2588
+0.2535 0.2587
+0.2534 0.2586
+0.2533 0.2585
+0.2533 0.2584
+0.2532 0.2583
+0.2531 0.2582
+0.2530 0.2581
+0.2529 0.2580
+0.2529 0.2579
+0.2528 0.2578
+0.2527 0.2577
+0.2526 0.2576
+0.2526 0.2575
+0.2525 0.2574
+0.2524 0.2573
+0.2524 0.2572
+0.2523 0.2572
+0.2522 0.2571
+0.2521 0.2570
+0.2521 0.2569
+0.2520 0.2568
+0.2519 0.2567
+0.2519 0.2566
+0.2518 0.2565
+0.2517 0.2565
+0.2517 0.2564
+0.2516 0.2563
+0.2516 0.2562
+0.2515 0.2561
+0.2514 0.2560
+0.2514 0.2560
+0.2513 0.2559
+0.2512 0.2558
+0.2512 0.2557
+0.2511 0.2557
+0.2511 0.2556
+0.2510 0.2555
+0.2510 0.2554
+0.2509 0.2554
+0.2508 0.2553
+0.2508 0.2552
+0.2507 0.2551
+0.2507 0.2551
+0.2506 0.2550
+0.2506 0.2549
+0.2505 0.2548
+0.2505 0.2548
+0.2504 0.2547
+0.2503 0.2546
+0.2503 0.2546
+0.2502 0.2545
+0.2502 0.2544
+0.2501 0.2544
+0.2501 0.2543
+0.2500 0.2542
+0.2500 0.2542
+0.2499 0.2541
+0.2499 0.2541
+0.2499 0.2540
+0.2498 0.2539
+0.2498 0.2539
+0.2497 0.2538
+0.2497 0.2537
+0.2496 0.2537
+0.2496 0.2536
+0.2495 0.2536
+0.2495 0.2535
+0.2494 0.2534
+0.2494 0.2534
+0.2493 0.2533
+0.2493 0.2533
+0.2493 0.2532
+0.2492 0.2532
+0.2492 0.2531
+0.2491 0.2530
+0.2491 0.2530
+0.2490 0.2529
+0.2490 0.2529
+0.2490 0.2528
+0.2489 0.2528
+0.2489 0.2527
+0.2488 0.2527
+0.2488 0.2526
+0.2488 0.2526
+0.2487 0.2525
+0.2487 0.2525
+0.2487 0.2524
+0.2486 0.2524
+0.2486 0.2523
+0.2485 0.2523
+0.2485 0.2522
+0.2485 0.2522
+0.2484 0.2521
+0.2484 0.2521
+0.2484 0.2520
+0.2483 0.2520
+0.2483 0.2519
+0.2482 0.2519
+0.2482 0.2518
+0.2482 0.2518
+0.2481 0.2517
+0.2481 0.2517
+0.2481 0.2516
+0.2480 0.2516
+0.2480 0.2516
+0.2480 0.2515
+0.2479 0.2515
+0.2479 0.2514
+0.2479 0.2514
+0.2478 0.2513
+0.2478 0.2513
+0.2478 0.2513
+0.2477 0.2512
+0.2477 0.2512
+0.2477 0.2511
+0.2476 0.2511
+0.2476 0.2510
+0.2476 0.2510
+0.2476 0.2510
+0.2475 0.2509
+0.2475 0.2509
+0.2475 0.2508
+0.2474 0.2508
+0.2474 0.2508
+0.2474 0.2507
+0.2473 0.2507
+0.2473 0.2506
+0.2473 0.2506
+0.2473 0.2506
+0.2472 0.2505
+0.2472 0.2505
+0.2472 0.2505
+0.2471 0.2504
+0.2471 0.2504
+0.2471 0.2503
+0.2471 0.2503
+0.2470 0.2503
+0.2470 0.2502
+0.2470 0.2502
+0.2470 0.2502
+0.2469 0.2501
+0.2469 0.2501
+0.2469 0.2501
+0.2468 0.2500
+0.2468 0.2500
+0.2468 0.2500
+0.2468 0.2499
+0.2467 0.2499
+0.2467 0.2498
+0.2467 0.2498
+0.2467 0.2498
+0.2466 0.2497
+0.2466 0.2497
+0.2466 0.2497
+0.2466 0.2497
+0.2465 0.2496
+0.2465 0.2496
+0.2465 0.2496
+0.2465 0.2495
+0.2465 0.2495
+0.2464 0.2495
+0.2464 0.2494
+0.2464 0.2494
+0.2464 0.2494
+0.2463 0.2493
+0.2463 0.2493
+0.2463 0.2493
+0.2463 0.2492
+0.2462 0.2492
+0.2462 0.2492
+0.2462 0.2492
+0.2462 0.2491
+0.2462 0.2491
+0.2461 0.2491
+0.2461 0.2490
+0.2461 0.2490
+0.2461 0.2490
+0.2460 0.2490
+0.2460 0.2489
+0.2460 0.2489
diff --git a/src/blackbody.c b/src/blackbody.c
new file mode 100644
index 0000000..abe03a9
--- /dev/null
+++ b/src/blackbody.c
@@ -0,0 +1,130 @@
+/**
+ * Copyright © 2016  Mattias Andrée <maandree@member.fsf.org>
+ * 
+ * This program 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.
+ * 
+ * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#include "blackbody.h"
+#include "macros.h"
+#include <math.h>
+#include <errno.h>
+
+
+
+/**
+ * Convert from CIE xyY to [0, 1] sRGB.
+ * 
+ * @param  x  The 'x' component.
+ * @param  y  The 'y' component.
+ * @param  Y  The 'Y' component.
+ * @param  r  Output parameter for the “red” value.
+ *            (Seriously, sRGB red is orange, just look at it fullscreen.)
+ * @param  g  Output parameter for the green value.
+ * @param  b  Output parameter for the blue value.
+ */
+static void
+ciexyy_to_srgb(double x, double y, double Y, double *r, double *g, double *b)
+{
+#define SRGB(C)  (((C) <= 0.0031308) ? (12.92 * (C)) : ((1.0 + 0.055) * pow((C), 1.0 / 2.4) - 0.055))
+	double X, Z;
+
+	/* Convert CIE xyY to CIE XYZ. */
+	X = Y * (y == 0.0 ? 0.0 : (x / y));
+	Z = Y * (y == 0.0 ? 0.0 : ((1.0 - x - y) / y));
+
+	/* Convert CIE XYZ to [0, 1] linear RGB. (ciexyz_to_linear) */
+	*r = ( 3.240450 * X) + (-1.537140 * Y) + (-0.4985320 * Z);
+	*g = (-0.969266 * X) + ( 1.876010 * Y) + ( 0.0415561 * Z);
+	*b = (0.0556434 * X) + (-0.204026 * Y) + ( 1.0572300 * Z);
+
+	/* Convert [0, 1] linear RGB to [0, 1] sRGB. */
+	SRGB(*r), SRGB(*g), SRGB(*b);
+}
+
+
+/**
+ * Perform linear interpolation (considered very good)
+ * between the CIE xyY values for two colour temperatures
+ * and convert the result to sRGB. The two colours should
+ * be the closest below the desired colour temperature,
+ * and the closest above the desired colour temperature.
+ * 
+ * @param  x1    The 'x' component for the low colour.
+ * @param  y1    The 'y' component for the low colour.
+ * @param  x2    The 'x' component for the high colour.
+ * @param  y2    The 'y' component for the high colour.
+ * @param  temp  The desired colour temperature.
+ * @param  r     Output parameter for the “red” value.
+ * @param  g     Output parameter for the green value.
+ * @param  b     Output parameter for the blue value.
+ */
+static void
+interpolate(double x1, double y1, double x2, double y2, double temp, double *r, double *g, double *b)
+{
+	double weight = fmod(temp, (double)DELTA_TEMPERATURE) / (double)DELTA_TEMPERATURE;
+	double x = x1 * (1 - weight) + x2 * weight;
+	double y = y1 * (1 - weight) + y2 * weight;
+	ciexyy_to_srgb(x, y, 1.0, r, g, b);
+}
+
+
+/**
+ * Get the [0, 1] sRGB values of a colour temperature.
+ * 
+ * @param   fd    File descriptor for the colour table.
+ * @param   temp  The desired colour temperature.
+ * @param   r     Output parameter for the “red” value.
+ * @param   g     Output parameter for the green value.
+ * @param   b     Output parameter for the blue value.
+ * @return        0 on succeess, -1 on error.
+ * 
+ * @throws  0     The file did not have the expected size.
+ * @throws  EDOM  The selected temperature is below 1000 K.
+ */
+int
+get_colour(int fd, long int temp, double *r, double *g, double *b)
+{
+	double values[10]; /* low:x,y,r,g,b + high:x,y,r,g,b */
+	off_t offset;
+	double max;
+
+	/* We do not have any values for above 40 000 K, but
+	 * the differences will be unnoticeable, perhaps even
+	 * unencodeable. */
+	if (temp > HIGHEST_TEMPERATURE)  temp = HIGHEST_TEMPERATURE;
+	/* Things do not glow below 1000 K. Yes, fire is hot! */
+	if (temp < LOWEST_TEMPERATURE)   t ((errno = EDOM));
+
+	/* Read table. */
+	offset = ((off_t)temp - LOWEST_TEMPERATURE) / DELTA_TEMPERATURE;
+	offset *= (off_t)(sizeof(values) / 2);
+	errno = 0; xpread(fd, values, sizeof(values), offset);
+
+	/* Get colour. */
+	if (temp % DELTA_TEMPERATURE)
+		interpolate(values[0], values[1], values[6], values[7], (double)temp, r, g, b);
+	else
+		*r = values[2], *g = values[3], *b = values[4];
+
+	/* Adjust colours for use. */
+	max = fmax(fmax(fabs(*r), fabs(*g)), fabs(*b));
+	if (max != 0)  *r /= max, *g /= max, *b /= max;
+	*r = *r > 0.0 ? *r : 0.0;
+	*g = *g > 0.0 ? *g : 0.0;
+	*b = *b > 0.0 ? *b : 0.0;
+
+	return 0;
+fail:
+	return -1;
+}
+
diff --git a/src/blackbody.h b/src/blackbody.h
new file mode 100644
index 0000000..499e9b8
--- /dev/null
+++ b/src/blackbody.h
@@ -0,0 +1,51 @@
+/**
+ * Copyright © 2016  Mattias Andrée <maandree@member.fsf.org>
+ * 
+ * This program 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.
+ * 
+ * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+
+/**
+ * The highest colour temperature in the table.
+ */
+#define LIBRED_HIGHEST_TEMPERATURE  40000
+
+/**
+ * The lowest colour temperature in the table.
+ */
+#define LIBRED_LOWEST_TEMPERATURE  1000
+
+/**
+ * The temperature difference between the colours in the table.
+ */
+#define LIBRED_DELTA_TEMPERATURE  100
+
+
+
+/**
+ * Get the [0, 1] sRGB values of a colour temperature.
+ * 
+ * @param   fd    File descriptor for the colour table.
+ * @param   temp  The desired colour temperature.
+ * @param   r     Output parameter for the “red” value.
+ * @param   g     Output parameter for the green value.
+ * @param   b     Output parameter for the blue value.
+ * @return        0 on succeess, -1 on error.
+ * 
+ * @throws  0     The file did not have the expected size.
+ * @throws  EDOM  The selected temperature is below 1000 K.
+ */
+int libred_get_colour(int fd, long int temp, double *r, double *g, double *b);
+
diff --git a/src/parse_10deg.c b/src/parse_10deg.c
new file mode 100644
index 0000000..ef1c3d5
--- /dev/null
+++ b/src/parse_10deg.c
@@ -0,0 +1,96 @@
+/**
+ * Copyright © 2016  Mattias Andrée <maandree@member.fsf.org>
+ * 
+ * This program 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.
+ * 
+ * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#include "blackbody.h"
+#include "macros.h"
+#include <sys/stat.h>
+
+
+
+/**
+ * The number of measured temperatures.
+ */
+#define TEMPERATURES  ((LIBRED_HIGHEST_TEMPERATURE - LIBRED_LOWEST_TEMPERATURE) / LIBRED_DELTA_TEMPERATURE + 1)
+
+
+
+/**
+ * Convert from CIE xyY to [0, 1] sRGB.
+ * 
+ * This function is identical to that in 'blackbody.c'.
+ * 
+ * @param  x  The 'x' component.
+ * @param  y  The 'y' component.
+ * @param  Y  The 'Y' component.
+ * @param  r  Output parameter for the “red” value.
+ *            (Seriously, sRGB red is orange, just look at it fullscreen.)
+ * @param  g  Output parameter for the green value.
+ * @param  b  Output parameter for the blue value.
+ */
+static void ciexyy_to_srgb(double x, double y, double Y, double *r, double *g, double *b)
+{
+#define SRGB(C)  (((C) <= 0.0031308) ? (12.92 * (C)) : ((1.0 + 0.055) * pow((C), 1.0 / 2.4) - 0.055))
+  double X, Z;
+  
+  /* Convert CIE xyY to CIE XYZ. */
+  X = Y * (y == 0.0 ? 0.0 : (x / y));
+  Z = Y * (y == 0.0 ? 0.0 : ((1.0 - x - y) / y));
+  
+  /* Convert CIE XYZ to [0, 1] linear RGB. (ciexyz_to_linear) */
+  *r = ( 3.240450 * X) + (-1.537140 * Y) + (-0.4985320 * Z);
+  *g = (-0.969266 * X) + ( 1.876010 * Y) + ( 0.0415561 * Z);
+  *b = (0.0556434 * X) + (-0.204026 * Y) + ( 1.0572300 * Z);
+  
+  /* Convert [0, 1] linear RGB to [0, 1] sRGB. */
+  SRGB(*r), SRGB(*g), SRGB(*b);
+}
+
+
+/**
+ * Create the lookup table of temperatures.
+ * 
+ * Standard input should be the file '10deg',
+ * standard output should be the table file
+ * and must be a regular file.
+ * 
+ * @param   argc  Should be 1.
+ * @param   argv  Should only contain the name of the process.
+ * @return        0 on success, 1 on error.
+ */
+int main(int argc, char *argv[])
+{
+#define x  (xyrgb[0])
+#define y  (xyrgb[1])
+#define r  (xyrgb[2])
+#define g  (xyrgb[3])
+#define b  (xyrgb[4])
+  
+  double xyrgb[5];
+  struct stat attr;
+  long int temp = LIBRED_LOWEST_TEMPERATURE;
+  
+  for (; fscanf(stdin, "%lf %lf\n", &x, &y) == 2; temp += LIBRED_DELTA_TEMPERATURE)
+    {
+      (temp == 6500) ? (r = g = b = 1.0) : ciexyy_to_srgb(x, y, 1.0, &r, &g, &b);
+      xwrite(STDOUT_FILENO, xyrgb, sizeof(xyrgb));
+    }
+  xwrite(STDOUT_FILENO, xyrgb, sizeof(xyrgb)); /* sugar */
+  t (fstat(STDOUT_FILENO, &attr));
+  return ((size_t)(attr.st_size) != (TEMPERATURES + 1) * 5 * sizeof(double));
+ fail:
+  return perror(argc ? *argv : "parse_10deg"), 1;
+}
+
diff --git a/src/solar.c b/src/solar.c
new file mode 100644
index 0000000..8a57507
--- /dev/null
+++ b/src/solar.c
@@ -0,0 +1,312 @@
+/**
+ * Copyright © 2016  Mattias Andrée <maandree@member.fsf.org>
+ * 
+ * This program 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.
+ * 
+ * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#include <math.h>
+#include <time.h>
+#include <errno.h>
+
+
+
+#if !defined(CLOCK_REALTIME_COARSE)
+# define CLOCK_REALTIME_COARSE  CLOCK_REALTIME
+#endif
+
+
+
+/**
+ * Get current Julian Centuries time (100 Julian days since J2000.)
+ * 
+ * @return  The current Julian Centuries time.
+ * 
+ * @throws  0  On success.
+ * @throws     Any error specified for clock_gettime(3) on error.
+ */
+static double
+julian_centuries()
+{
+	struct timespec now;
+	double tm;
+	if (clock_gettime(CLOCK_REALTIME_COARSE, &now))  return 0.0;
+	tm = (double)(now.tv_nsec) / 1000000000.0 + (double)(now.tv_sec);
+	tm = (tm / 86400.0 + 2440587.5 - 2451545.0) / 36525.0;
+	return errno = 0, tm;
+}
+
+/**
+ * Convert a Julian Centuries timestamp to a Julian Day timestamp.
+ * 
+ * @param   tm  The time in Julian Centuries
+ * @return      The time in Julian Days
+ */
+static inline double
+julian_centuries_to_julian_day(double tm)
+{
+	return tm * 36525.0 + 2451545.0;
+}
+
+
+/**
+ * Convert an angle (or otherwise) from degrees to radians.
+ * 
+ * @param  deg  The angle in degrees.
+ * @param       The angle in radians.
+ */
+static inline double
+radians(double deg)
+{
+	return deg * (double)M_PI / 180.0;
+}
+
+/**
+ * Convert an angle (or otherwise) from radians to degrees.
+ * 
+ * @param  rad  The angle in radians.
+ * @param       The angle in degrees.
+ */
+static inline double
+degrees(double rad)
+{
+	return rad * 180.0 / (double)M_PI;
+}
+
+
+/**
+ * Calculates the Sun's elevation from the solar hour angle
+ * 
+ * @param   longitude    The longitude in degrees eastwards.
+ *                       from Greenwich, negative for westwards.
+ * @param   declination  The declination, in radians.
+ * @param   hour_angle   The solar hour angle, in radians.
+ * @return               The Sun's elevation, in radians.
+ */
+static inline double
+elevation_from_hour_angle(double latitude, double declination, double hour_angle)
+{
+	double rc = cos(radians(latitude));
+	rc *= cos(hour_angle) * cos(declination);
+	rc += sin(radians(latitude)) * sin(declination);
+	return asin(rc);
+}
+
+/**
+ * Calculates the Sun's geometric mean longitude.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The Sun's geometric mean longitude in radians.
+ */
+static inline double
+sun_geometric_mean_longitude(double tm)
+{
+	double rc = fmod(pow(0.0003032 * tm, 2.0) + 36000.76983 * tm + 280.46646, 360.0);
+#if defined(TIMETRAVELLER)
+	rc = rc < 0.0 ? (rc + 360.0) : rc;
+#endif
+	return radians(rc);
+}
+
+/**
+ * Calculates the Sun's geometric mean anomaly.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The Sun's geometric mean anomaly in radians.
+ */
+static inline double
+sun_geometric_mean_anomaly(double tm)
+{
+	return radians(pow(-0.0001537 * tm, 2.0) + 35999.05029 * tm + 357.52911);
+}
+
+/**
+ * Calculates the Earth's orbit eccentricity.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The Earth's orbit eccentricity.
+ */
+static inline double
+earth_orbit_eccentricity(double tm)
+{
+	return pow(-0.0000001267 * tm, 2.0) - 0.000042037 * tm + 0.016708634;
+}
+
+/**
+ * Calculates the Sun's equation of the centre, the difference
+ * between the true anomaly and the mean anomaly.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The Sun's equation of the centre, in radians.
+ */
+static inline double
+sun_equation_of_centre(double tm)
+{
+	double a = sun_geometric_mean_anomaly(tm), rc;
+	rc  = sin(1.0 * a) * (pow(-0.000014 * tm, 2.0) - 0.004817 * tm + 1.914602);
+	rc += sin(2.0 * a) * (-0.000101 * tm + 0.019993);
+	rc += sin(3.0 * a) * 0.000289;
+	return radians(rc);
+}
+
+/**
+ * Calculates the Sun's real longitudinal position.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The longitude, in radians.
+ */
+static inline double
+sun_real_longitude(double tm)
+{
+	return sun_geometric_mean_longitude(tm) + sun_equation_of_centre(tm);
+}
+
+/**
+ * Calculates the Sun's apparent longitudinal position.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The longitude, in radians.
+ */
+static inline double
+sun_apparent_longitude(double tm)
+{
+    double rc = degrees(sun_real_longitude(tm)) - 0.00569;
+    return radians(rc - 0.00478 * sin(radians(-1934.136 * tm + 125.04)));
+}
+
+/**
+ * Calculates the mean ecliptic obliquity of the Sun's
+ * apparent motion without variation correction.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The uncorrected mean obliquity, in radians.
+ */
+static double
+mean_ecliptic_obliquity(double tm)
+{
+	double rc = pow(0.001813 * tm, 3.0) - pow(0.00059 * tm, 2.0) - 46.815 * tm + 21.448;
+	return radians(23.0 + (26.0 + rc / 60.0) / 60.0);
+}
+
+/**
+ * Calculates the mean ecliptic obliquity of the Sun's
+ * parent motion with variation correction.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The mean obliquity, in radians.
+ */
+static double
+corrected_mean_ecliptic_obliquity(double tm)
+{
+	double rc = 0.00256 * cos(radians(-1934.136 * tm + 125.04));
+	return radians(rc + degrees(mean_ecliptic_obliquity(tm)));
+}
+
+/**
+ * Calculates the Sun's declination.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The Sun's declination, in radian.
+ */
+static inline double
+solar_declination(double tm)
+{
+	double rc = sin(corrected_mean_ecliptic_obliquity(tm));
+	return asin(rc * sin(sun_apparent_longitude(tm)));
+}
+
+/**
+ * Calculates the equation of time, the discrepancy
+ * between apparent and mean solar time.
+ * 
+ * @param   tm  The time in Julian Centuries.
+ * @return      The equation of time, in degrees.
+ */
+static inline double
+equation_of_time(double tm)
+{
+	double l = sun_geometric_mean_longitude(tm);
+	double e = earth_orbit_eccentricity(tm);
+	double m = sun_geometric_mean_anomaly(tm);
+	double y = pow(tan(corrected_mean_ecliptic_obliquity(tm) / 2.0), 2.0);
+	double rc = y * sin(2.0 * l);
+	rc += (4.0 * y * cos(2.0 * l) - 2.0) * e * sin(m);
+	rc -= pow(0.5 * y, 2.0) * sin(4.0 * l);
+	rc -= pow(1.25 * e, 2.0) * sin(2.0 * m);
+	return 4.0 * degrees(rc);
+}
+
+/**
+ * Calculates the Sun's elevation as apparent.
+ * from a geographical position.
+ * 
+ * @param   tm         The time in Julian Centuries.
+ * @param   latitude   The latitude in degrees northwards from 
+ *                     the equator, negative for southwards.
+ * @param   longitude  The longitude in degrees eastwards from
+ *                     Greenwich, negative for westwards.
+ * @return             The Sun's apparent elevation at the specified time as seen
+ *                     from the specified position, measured in radians.
+ */
+static inline double
+solar_elevation_from_time(double tm, double latitude, double longitude)
+{
+	double rc = julian_centuries_to_julian_day(tm);
+	rc = (rc - round(rc) - 0.5) * 1440;
+	rc = 720.0 - rc - equation_of_time(tm);
+	rc = radians(rc / 4.0 - longitude);
+	return elevation_from_hour_angle(latitude, solar_declination(tm), rc);
+}
+
+
+/**
+ * Calculates the Sun's elevation as apparent.
+ * from a geographical position.
+ * 
+ * @param   latitude   The latitude in degrees northwards from 
+ *                     the equator, negative for southwards.
+ * @param   longitude  The longitude in degrees eastwards from
+ *                     Greenwich, negative for westwards.
+ * @return             The Sun's apparent elevation as seen, right now,
+ *                     from the specified position, measured in degrees.
+ * 
+ * @throws  0  On success.
+ * @throws     Any error specified for clock_gettime(3) on error.
+ */
+double
+solar_elevation(double latitude, double longitude)
+{
+	double tm = julian_centuries();
+	return errno ? -1 : degrees(solar_elevation_from_time(tm, latitude, longitude));
+}
+
+
+/**
+ * Exit if time the is before year 0 in J2000.
+ * 
+ * @return  0 on success, -1 on error.
+ */
+#if defined(TIMETRAVELLER)
+int
+check_timetravel(void)
+{
+	struct timespec now;
+	if (clock_gettime(CLOCK_REALTIME, &now))  return -1;
+	if (now.tv_nsec < (time_t)946728000L)
+		fprintf(stderr, "We have detected that you are a time-traveller"
+		                "(or your clock is not configured correctly.)"
+		                "Please recompile with -DTIMETRAVELLER"
+	                        "(or correct your clock.)"), exit(1);
+	return 0;
+}
+#endif
+
diff --git a/src/solar.h b/src/solar.h
new file mode 100644
index 0000000..8e94cf9
--- /dev/null
+++ b/src/solar.h
@@ -0,0 +1,107 @@
+/**
+ * Copyright © 2016  Mattias Andrée <maandree@member.fsf.org>
+ * 
+ * This program 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.
+ * 
+ * This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+/**
+ * Approximate apparent size of the Sun in degrees.
+ */
+#define LIBRED_SOLAR_APPARENT_RADIUS  (32.0 / 60.0)
+
+
+/**
+ * The Sun's elevation at sunset and sunrise, measured in degrees.
+ */
+#define LIBRED_SOLAR_ELEVATION_SUNSET_SUNRISE   (0.0)
+
+/**
+ * The Sun's elevation at civil dusk and civil dawn, measured in degrees.
+ */
+#define LIBRED_SOLAR_ELEVATION_CIVIL_DUSK_DAWN   (-6.0)
+
+/**
+ * The Sun's elevation at nautical dusk and nautical dawn, measured in degrees.
+ */
+#define LIBRED_SOLAR_ELEVATION_NAUTICAL_DUSK_DAWN  (-12.0)
+
+/**
+ * The Sun's elevation at astronomical dusk and astronomical dawn, measured in degrees.
+ */
+#define LIBRED_SOLAR_ELEVATION_ASTRONOMICAL_DUSK_DAWN  (-18.0)
+
+
+/**
+ * Test whether it is twilight.
+ * 
+ * @param   ELEV:double  The current elevation.
+ * @return               1 if is twilight, 0 otherwise.
+ */
+#define LIBRED_IS_TWILIGHT(ELEV)  ((-18.0 <= (ELEV)) && ((ELEV) <= 0.0))
+
+/**
+ * Test whether it is civil twilight.
+ * 
+ * @param   ELEV:double  The current elevation.
+ * @return               1 if is civil twilight, 0 otherwise.
+ */
+#define LIBRED_IS_CIVIL_TWILIGHT(ELEV)  ((-6.0 <= (ELEV)) && ((ELEV) <= 0.0))
+
+/**
+ * Test whether it is nautical twilight.
+ * 
+ * @param   ELEV:double  The current elevation.
+ * @return               1 if is nautical twilight, 0 otherwise.
+ */
+#define LIBRED_IS_NAUTICAL_TWILIGHT(ELEV)  ((-12.0 <= (ELEV)) && ((ELEV) <= -6.0))
+
+/**
+ * Test whether it is astronomical twilight.
+ * 
+ * @param   ELEV:double  The current elevation.
+ * @return               1 if is astronomical twilight, 0 otherwise.
+ */
+#define LIBRED_IS_ASTRONOMICAL_TWILIGHT(ELEV)  ((-18.0 <= (ELEV)) && ((ELEV) <= -12.0))
+
+
+
+/**
+ * Calculates the Sun's elevation as apparent.
+ * from a geographical position.
+ * 
+ * @param   latitude   The latitude in degrees northwards from 
+ *                     the equator, negative for southwards.
+ * @param   longitude  The longitude in degrees eastwards from
+ *                     Greenwich, negative for westwards.
+ * @return             The Sun's apparent elevation as seen, right now,
+ *                     from the specified position, measured in degrees.
+ * 
+ * @throws  0  On success.
+ * @throws     Any error specified for clock_gettime(3) on error.
+ */
+double libred_solar_elevation(double latitude, double longitude);
+
+
+/**
+ * Exit if time the is before year 0 in J2000.
+ * 
+ * @return  0 on success, -1 on error.
+ */
+#if defined(TIMETRAVELLER)
+int libred_check_timetravel(void);
+#else
+# define libred_check_timetravel()  0
+#endif
+