From bcbca9e0650d8a1a2dfc092c7561faa08169a3c8 Mon Sep 17 00:00:00 2001
From: Kenny Levinsen <kl@kl.wtf>
Date: Wed, 11 Nov 2020 16:40:37 +0100
Subject: [PATCH] Use illuminant D for the majority of temperatures

D65, the natural whitepoint that wlsunset assumes, is defined on
illuminant D, which simulates daylight with atmospheric effects.
However, we used planckian locus for all values under 6500 K, which
meant that there was a significant jump in color - specifically, a
sudden reduction in green and blue - as we started reducing the color
temperature.

Instead, we purely use illuminant D down to 4000 K where it is well
defined. Below 4000 K, illuminant D starts unintentionally approaching
the planckian locus, before finally breaking completely at 2000 K. We
extend the boundary of illuminant D to 2500 K and perform a
sine-smoothed transition to planckian locus from 4000 K to 2500 K to
extend the range of wlsunset down to 1667 K, where planckian locus ends
and we finally give up.

The end-result is a smooth transition along all valid temperature
values, with no sudden jumps as we had before. However, we do end up
with slightly greener/bluer colors than earlier. We'll have to see how
that holds up.
---
 color_math.c | 37 +++++++++++++++++++++++++++++++++----
 1 file changed, 33 insertions(+), 4 deletions(-)

diff --git a/color_math.c b/color_math.c
index d7a8fbb..b850b96 100644
--- a/color_math.c
+++ b/color_math.c
@@ -71,9 +71,17 @@ enum sun_condition calc_sun(struct tm *tm, double latitude, struct sun *sun) {
 		condition(latitude, decl) : NORMAL;
 }
 
+/*
+ * Illuminant D, or daylight locus, is is a "standard illuminant" used to
+ * describe natural daylight. It is on this locus that D65, the whitepoint used
+ * by most monitors and assumed by wlsunset, is defined.
+ *
+ * This approximation is strictly speaking only well-defined between 4000K and
+ * 25000K, but we stretch it a bit further down for transition purposes.
+ */
 static int illuminant_d(int temp, double *x, double *y) {
 	// https://en.wikipedia.org/wiki/Standard_illuminant#Illuminant_series_D
-	if (temp >= 4000 && temp <= 7000) {
+	if (temp >= 2500 && temp <= 7000) {
 		*x = 0.244063 +
 			0.09911e3 / temp +
 			2.9678e6 / pow(temp, 2) -
@@ -91,7 +99,15 @@ static int illuminant_d(int temp, double *x, double *y) {
 	return 0;
 }
 
+/*
+ * Planckian locus, or black body locus, describes the color of a black body at
+ * a certain temperatures. This is not entirely equivalent to daylight due to
+ * atmospheric effects.
+ *
+ * This approximation is only valid from 1667K to 25000K.
+ */
 static int planckian_locus(int temp, double *x, double *y) {
+	// https://en.wikipedia.org/wiki/Planckian_locus#Approximation
 	if (temp >= 1667 && temp <= 4000) {
 		*x = -0.2661239e9 / pow(temp, 3) -
 			0.2343589e6 / pow(temp, 2) +
@@ -164,10 +180,23 @@ void calc_whitepoint(int temp, double *rw, double *gw, double *bw) {
 	}
 
 	double x = 1.0, y = 1.0;
-	if (temp > 1667 && temp <= 6500) {
-		planckian_locus(temp, &x, &y);
-	} else if (temp >= 6500 && temp <= 25000) {
+	if (temp >= 25000) {
+		illuminant_d(25000, &x, &y);
+	} else if (temp >= 4000) {
 		illuminant_d(temp, &x, &y);
+	} else if (temp >= 2500) {
+		double x1, y1, x2, y2;
+		illuminant_d(temp, &x1, &y1);
+		planckian_locus(temp, &x2, &y2);
+
+		double factor = (4000 - temp) / 1500;
+		double sinefactor = (cos(M_PI*factor) + 1.0) / 2.0;
+		x = x1 * sinefactor + x2 * (1.0 - sinefactor);
+		y = y1 * sinefactor + y2 * (1.0 - sinefactor);
+	} else if (temp >= 1667) {
+		planckian_locus(temp, &x, &y);
+	} else {
+		planckian_locus(1667, &x, &y);
 	}
 	double z = 1.0 - x - y;