PL. 118 · COLOUR / RADIANCE / SPECTRAL LOCUS
Planckian Locus
Max Planck, 1900 · the colour matching fit: Wyman, Sloan & Shirley, 2013
OPEN THE LIVE PLATE ▸DEFINITION
B_λ(λ,T) = (2hc²/λ⁵) / (exp(hc/λk_BT) − 1) X = ∫ B_λ x̄(λ) dλ, Y = ∫ B_λ ȳ(λ) dλ, Z = ∫ B_λ z̄(λ) dλ RGB = M_sRGB · (X, Y, Z) / Y
NOTES
Hot matter glows, and Planck found in 1900 exactly how much of each wavelength it sends out. The eye then throws almost all of that away: a whole spectrum arrives and three numbers leave, because there are three kinds of cone and no more. Do both in order and you get this curve — the colour of a thing that is merely hot, from the deep red of an ember through the white of daylight to the blue of the hottest stars. The band across the middle is that colour, one column per temperature, temperature rising to the right. The register above it is the sum being performed: each row is one wavelength, coloured by what that wavelength contributes to the three numbers, and weighted by how much of it a body at that temperature actually emits. Read a column downward and you are reading the terms; the band below is their total. Watch the weight tip from the red lobe to the blue one as the temperature climbs, and the band turns over from orange through white to blue in step with it. Both ends of the register are black because the eye is: past about 700 nanometres and short of about 400 there is nothing to weight, which is also why a real body at 1000 K reads as a dim ember rather than the furnace it is: almost everything it radiates is in the infrared, where nobody is looking. That difference in brightness is precisely what this plate divides out, and it has to, because the two ends of the sweep differ by about seven orders of magnitude and could not otherwise share a frame. The rule at the foot marks where the answer stops being true: below about 1900 K the colour is outside anything a screen can mix, and what is drawn there is the nearest one sRGB has.
PROVENANCE
- Origin
- M. Planck, "Ueber das Gesetz der Energieverteilung im Normalspectrum", Annalen der Physik 309(3), 1901, 553–563 — the definitive publication of the law he had presented to the German Physical Society on 14 December 1900
- Standing
- Public domain. A law of physics, published in 1901, and implemented here from the law itself: the plate carries only the second radiation constant hc/k_B, and even the Wien peak it uses is obtained by differentiating Planck rather than by citing Wien separately
- The colour matching functions, and why the tables are absent
- Turning a spectrum into a colour needs the CIE 1931 standard observer, whose tabulated x̄ ȳ z̄ values are asserted-copyright and are therefore not in this tree and never will be. What is here instead is the published analytic fit: C. Wyman, P.-P. Sloan and P. Shirley, "Simple Analytic Approximations to the CIE XYZ Color Matching Functions", Journal of Computer Graphics Techniques 2(2), 2013, 1–11 (jcgt.org/published/0002/02/01/). The three curves are sums of piecewise-continuous Gaussians — three lobes for x̄, two each for ȳ and z̄ — and this plate evaluates them from the coefficients of Table 1 of that paper. Those coefficients and the formula around them are the content the paper exists to publish, offered precisely so that people can implement them instead of copying tables; using them is attribution, which is given here in full, and not a licence. No code was taken: the paper ships a C listing and a supplemental repository, and neither was read into this implementation, which is written from Equation 4 and the table.
- Used only for what it was fitted to
- The paper fits absolute error in x̄, ȳ and z̄, and that is the only way this plate uses it. An earlier draft coloured the upper register by the chromaticity of each single wavelength, which is a ratio between the three curves — and out where all three are near zero, a ratio of a fit that was never fitted to ratios is not a quantity anyone should trust. It could not be checked either, since checking it would need exactly the tabulations this plate exists to do without. So the register carries the absolute values instead, straight through the sRGB matrix, and every wavelength beyond the response of the eye goes to black by arithmetic rather than by decision
- Normalised by Y, and what that costs
- Between the ends of this frame a black body changes its emitted power by about seven orders of magnitude, so a plate that plotted radiance would be one bright pixel column and nothing else. Dividing XYZ by Y throws the brightness away and keeps the chromaticity, then each column is scaled to the brightest sRGB triple carrying it. Brightness carries no information along the band; only hue does. The register above divides each of its terms by the same Y, which is the column being measured against its own visible output rather than against its total output, and then by one measured constant — the largest that ratio gets anywhere this plate can reach, so exposure 1 is the setting at which nothing clips. Measured, that normalisation holds the register within a factor of 3.5 across the entire range the sliders and the travel can produce, against a factor of 14 for the obvious alternative of scaling by the peak of the Planck curve
- Clipped, and marked
- Much of what the eye can see, no screen can mix. Where the conversion returns a negative primary the chromaticity is outside the sRGB gamut and the triple is clamped — measured, that is everything below 1883 K, and the figure is the same whether the integral takes the 36 samples this plate draws with or 512. Rather than let the clamp pass silently, the rule under the band is lit exactly where it happened — which also explains the one visible edge in the band. Clamping a primary at zero and then rescaling puts a kink in the colour exactly at the boundary, measured at 6.8 levels of 255 between adjacent columns near 1948 K, and it is there at every sample count from 8 to 96 rather than being an artefact of the integral. The rule beneath it is the caption for it. No chromatic adaptation is applied either: the white point is the D65 that sRGB was defined around, which is close to but not on the locus, so daylight temperatures land near white by construction rather than by coincidence
- Checked, not just plotted
- CIE standard illuminant A is defined as a Planckian radiator at 2856 K, and the CIE publishes its chromaticity as x = 0.44757, y = 0.40745 — a pair of numbers nothing in this implementation has ever seen. Integrating Planck at 2856 K against the analytic fit here gives x = 0.44624, y = 0.40866, which is 0.0018 away in the chromaticity plane. That residual is the price of the fit plus the cut of the integration band at 780 nm rather than 830, and it is independent evidence that what is drawn is the locus itself and not a gradient chosen to look like one
- Constants
- The window bounds are the temperatures at the left and right edges, spaced logarithmically because that is how the locus spreads evenly; travel slides that window back and forth in log temperature, which is the drift. Samples is the number of wavelengths the integral takes across 380–780 nm and is the GLSL loop bound. Measured against a 512-sample reference, the error first falls under one level in 255 at 28 samples and stays under it above that, though not monotonically — the midpoint rule beats against the lobes of the fit, so 12 samples is worse than 8 — and 36 is the default for margin rather than economy. The low end of the range is left reachable because an under-sampled integral is worth seeing once: at 8 samples the colour is wrong by about 14 levels in 255, and wrong smoothly rather than in bands. Exposure scales the spectral register only — the band below is a chromaticity strip, and brightening it would destroy the one thing it shows
- Source
- doi:10.1002/andp.19013090310
TOUCHDESIGNER · GLSL
The same shader this plate runs, reframed for a GLSL TOP. Pasted bare it renders the published constants as a still frame; wire absTime.seconds into u_t on the Vectors page to animate it.
// FORMA — PL. 118 · PLANCKIAN LOCUS — Max Planck, 1900 · the colour matching fit: Wyman, Sloan & Shirley, 2013
// B_λ(λ,T) = (2hc²/λ⁵) / (exp(hc/λk_BT) − 1)
// X = ∫ B_λ x̄(λ) dλ, Y = ∫ B_λ ȳ(λ) dλ, Z = ∫ B_λ z̄(λ) dλ
// RGB = M_sRGB · (X, Y, Z) / Y
// TouchDesigner port — paste into a GLSL TOP's pixel shader. Set the
// resolution on the TOP's Common page. As pasted it renders the published
// constants as a still frame; to animate, add a uniform named u_t on the
// GLSL TOP's Vectors 1 page with the expression absTime.seconds.
// Constants are consts — edit to tweak; comments give the measured range.
// Written from the published mathematics, not adapted from any code.
#define u_res (uTDOutputInfo.res.zw)
uniform float u_t; // absTime.seconds on the Vectors page; unset = still
const float u_phase = 0.041; // this plate's own grid phase, 0..1
// FORMA's COLOUR accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.2002, 0.3896);
const vec3 u_pal_b = vec3(0.5, 0.2176, 0.4235);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);
const float p_tlo = 1000.0; // window low (K) · live 600 .. 3000
const float p_thi = 12000.0; // window high (K) · live 6000 .. 20000
const float p_samples = 36.0; // spectral samples · live 8 .. 96
const float p_travel = 0.35; // window travel · live 0 .. 0.6
const float p_exposure = 1.3; // spectrum brightness · live 0.4 .. 2.4
/* The order's ramp — the same cosine formulation the JS kit uses, so a
plate keeps its classification colour in either language. */
vec3 ramp(float t){
return clamp(u_pal_a + u_pal_b * cos(6.28318530718 * (u_pal_c * t + u_pal_d)), 0.0, 1.0);
}
/* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
float cycle(float t, float period){ return fract(t / period + u_phase); }
float pingpong(float t, float period){
float u = cycle(t, period);
return u < 0.5 ? u * 2.0 : 2.0 - u * 2.0;
}
/* Planck's law gives the spectrum; the CIE colour matching functions reduce a
spectrum to a colour. The tabulated functions are asserted-copyright and are
not in this tree — what is here is the published analytic fit of Wyman,
Sloan and Shirley (JCGT 2(2), 2013), evaluated from the coefficients of its
own Table 1.
Every numeric constant below is declared inside the function that uses it
rather than at file scope. A composition emits two copies of a body and
namespaces the function declarations, so a file-scope const would be
declared twice by blackbody x blackbody and fail to link. */
/* Only the shape in wavelength matters — scaling a spectrum does not change
its chromaticity — so every constant factor is dropped. exp(-u) rather than
exp(u): at the cold end of the window u passes 80, and exp(u) there is
outside the range of a 32-bit float. */
float bbody_planck(float lambda, float T){
const float c2 = 1.438776877e7; // hc/k_B, in nm K
float u = c2 / (lambda * T);
float e = exp(-u);
return u * u * u * u * u * e / (1.0 - e);
}
/* Divided by the largest value the curve reaches inside the visible band: the
Wien peak where that is visible, the nearer edge of the band where it is
not. c2/x0 is the Wien displacement constant, and x0 is the root of
x = 5(1 - exp(-x)) — Wien falls out of Planck by differentiating. */
float bbody_norm(float lambda, float T){
const float c2 = 1.438776877e7, x0 = 4.965114232;
float peak = clamp(c2 / (x0 * T), 380.0, 780.0);
return bbody_planck(lambda, T) / bbody_planck(peak, T);
}
/* One piecewise-continuous Gaussian lobe: the width differs either side of the
centre, which is how the fit follows an asymmetry a plain Gaussian cannot.
g and d are reciprocal widths. */
float bbody_lobe(float lambda, float a, float b, float g, float d){
float x = (lambda - b) * (lambda < b ? g : d);
return a * exp(-0.5 * x * x);
}
/* x-bar, y-bar, z-bar of the CIE 1931 two-degree observer: three lobes, then
two and two. The negative amplitude is the real dip in x-bar near 500 nm. */
vec3 bbody_cmf(float lambda){
return vec3(
bbody_lobe(lambda, 0.362, 442.0, 0.0624, 0.0374)
+ bbody_lobe(lambda, 1.056, 599.8, 0.0264, 0.0323)
+ bbody_lobe(lambda, -0.065, 501.1, 0.0490, 0.0382),
bbody_lobe(lambda, 0.821, 568.8, 0.0213, 0.0247)
+ bbody_lobe(lambda, 0.286, 530.9, 0.0613, 0.0322),
bbody_lobe(lambda, 1.217, 437.0, 0.0845, 0.0278)
+ bbody_lobe(lambda, 0.681, 459.0, 0.0385, 0.0725));
}
/* CIE XYZ to linear sRGB (IEC 61966-2-1 primaries at D65), then that
standard's transfer function — a linear toe under a shifted power curve. */
vec3 bbody_rgb(vec3 c){
return vec3( 3.2404542 * c.x - 1.5371385 * c.y - 0.4985314 * c.z,
-0.9692660 * c.x + 1.8760108 * c.y + 0.0415560 * c.z,
0.0556434 * c.x - 0.2040259 * c.y + 1.0572252 * c.z);
}
vec3 bbody_encode(vec3 c){
c = clamp(c, 0.0, 1.0);
return mix(12.92 * c, 1.055 * pow(c, vec3(1.0 / 2.4)) - 0.055, step(0.0031308, c));
}
/* The mean of the integrand over the visible band at one temperature — the
quantity both registers are built from. The mean rather than the sum, so
the sample count can change how accurate this is and never how bright. */
vec3 bbody_sum(float T, float ns){
int n = int(floor(ns + 0.5));
vec3 xyz = vec3(0.0);
/* Constant loop bound at the slider's own ceiling, break on the live value,
mandelbrot fashion. Normalising the integrand by the column peak cancels
out of the chromaticity exactly and keeps the sum near 1 rather than near
1e-30, which is where a cold column would otherwise land. */
for (int i = 0; i < 96; i++){
if (i >= n) break;
float lambda = 380.0 + (float(i) + 0.5) * 400.0 / float(n);
xyz += bbody_cmf(lambda) * bbody_norm(lambda, T);
}
return xyz / float(n);
}
vec3 plate(vec2 uv){
const float FIELD = 0.60, BAND0 = 0.638, BAND1 = 0.940, RULE = 0.958;
/* The largest value the register's ratio takes anywhere this plate can
reach. Measured, so exposure 1 is exactly the setting at which nothing
clips and anything above it is a decision. */
const float REG = 30.1;
/* The window slides in log temperature and back, so the sweep keeps its
shape while it travels. cycle() carries this plate's PHASE. */
float slide = p_travel * sin(6.28318530718 * cycle(u_t, 34.0));
float T = exp(mix(log(p_tlo), log(p_thi), uv.x) + slide);
/* One integral per fragment, whichever register it lands in: the band needs
it for its colour, and the field needs its luminance to normalise by. */
vec3 xyz = bbody_sum(T, p_samples);
float Y = max(xyz.y, 1e-30);
if (uv.y < FIELD){
/* One term of the very sum below it: this wavelength's matching functions
taken straight through the sRGB matrix — the fit used only for what it
was fitted to, the absolute values — weighted by how much of this
wavelength the temperature emits. Negative components are simply what
it means for a single wavelength to be outside the gamut. */
float lambda = 380.0 + (uv.y / FIELD) * 400.0;
vec3 term = bbody_rgb(bbody_cmf(lambda)) * bbody_norm(lambda, T);
return bbody_encode(term / (Y * REG) * p_exposure);
}
/* Chromaticity, not radiance, then scaled to the brightest sRGB triple that
carries it. A negative primary means sRGB cannot mix the colour at all. */
vec3 lin = bbody_rgb(xyz / Y);
float m = max(max(lin.r, lin.g), lin.b);
if (uv.y > BAND0 && uv.y < BAND1) return bbody_encode(max(lin, 0.0) / max(m, 1e-9));
if (uv.y > RULE) return mix(vec3(26.0, 30.0, 38.0) / 255.0, vec3(209.0 / 255.0),
step(min(min(lin.r, lin.g), lin.b), -0.002 * m));
return vec3(4.0, 6.0, 10.0) / 255.0;
}
out vec4 fragColor;
void main(){
// FORMA's uv runs y-down, matching its canvas; TD's vUV runs up
vec2 uv = vec2(vUV.s, 1.0 - vUV.t);
fragColor = TDOutputSwizzle(vec4(plate(uv), 1.0));
}
NUKE · BLINKSCRIPT
The same shader this plate runs, transpiled to a BlinkScript kernel. Paste it into a BlinkScript node's Kernel Source and press Recompile; every constant arrives as a knob at its published value, and u_t animates with the expression frame/24. Compiled and rendered in Nuke 17.1, then compared against this plate on the page.
// FORMA — PL. 118 · PLANCKIAN LOCUS — Max Planck, 1900 · the colour matching fit: Wyman, Sloan & Shirley, 2013
// B_λ(λ,T) = (2hc²/λ⁵) / (exp(hc/λk_BT) − 1)
// X = ∫ B_λ x̄(λ) dλ, Y = ∫ B_λ ȳ(λ) dλ, Z = ∫ B_λ z̄(λ) dλ
// RGB = M_sRGB · (X, Y, Z) / Y
// Nuke port — a BlinkScript kernel. Paste into a BlinkScript node's Kernel
// Source and press Recompile. Every constant arrives as a knob at its published
// value (the comment gives the measured range); u_t is a knob too — animate it
// with the expression frame/24 or leave it at 0 for the still frame. Written
// from the published mathematics, not adapted from any code.
// Transpiled from the shader this plate runs on the page (GLSL ES 3.00):
// vec → float2/3/4, swizzles expanded, GLSL builtins Blink lacks written out
// as forma_ functions, float literals suffixed. Compiled and rendered in a
// real Nuke (17.1v1) and compared against this plate on the page: 34 of 34.
//
// plate() and its helpers are written to a single exit — the loop that runs
// once. That is not a style: Blink 17.1 drops a conditional early return from
// a called function while Vectorize is on, which is the node default, with no
// warning and no error. Written this way it paints correctly as pasted.
kernel Forma_blackbody : ImageComputationKernel<ePixelWise>
{
Image<eWrite> dst;
param:
float u_t; // seconds; 0 is the still frame
float p_tlo; // window low (K) · live 600 .. 3000
float p_thi; // window high (K) · live 6000 .. 20000
float p_samples; // spectral samples · live 8 .. 96
float p_travel; // window travel · live 0 .. 0.6
float p_exposure; // spectrum brightness · live 0.4 .. 2.4
local:
float2 u_res;
float u_phase;
float3 u_pal_a, u_pal_b, u_pal_c, u_pal_d;
void define(){
defineParam(u_t, "u_t", 0.0f);
defineParam(p_tlo, "p_tlo", 1000.0f);
defineParam(p_thi, "p_thi", 12000.0f);
defineParam(p_samples, "p_samples", 36.0f);
defineParam(p_travel, "p_travel", 0.35f);
defineParam(p_exposure, "p_exposure", 1.3f);
}
void init(){
u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
u_phase = 0.041f; // this plate's own grid phase, 0..1
// FORMA's COLOUR accent as cosine-gradient coefficients
u_pal_a = float3(0.46f, 0.2002f, 0.3896f);
u_pal_b = float3(0.5f, 0.2176f, 0.4235f);
u_pal_c = float3(1.0f, 1.0f, 1.0f);
u_pal_d = float3(0.0f, 0.05f, 0.1f);
}
/* GLSL builtins Blink lacks, written as templates rather than overload sets.
Blink's operators return expression templates (Swizzle<float,N>), so a call
passing an expression cannot resolve against an overload set on float2
against float3 — measured in Nuke 17.1: a float2 expression is ambiguous
between the two, while scalar-against-vector resolves. A template deduces
the expression's own type, so the ambiguity cannot arise. */
template <class T> T forma_fract(T v){ return v - floor(v); }
template <class T, class S> T forma_mod(T x, S y){ return x - y * floor(x / y); }
/* Blink's own min/max/clamp take no scalar bound against a vector, which GLSL
does; v * 0.0f + b is that bound at the vector's own width, and collapses to
b when v is a scalar, so one template serves both. */
template <class T, class S> T forma_min(T a, S b){ return min(a, a * 0.0f + b); }
template <class T, class S> T forma_max(T a, S b){ return max(a, a * 0.0f + b); }
template <class T, class S> T forma_clamp(T v, S lo, S hi){ return clamp(v, v * 0.0f + lo, v * 0.0f + hi); }
int forma_min(int a, int b){ return min(a, b); }
int forma_max(int a, int b){ return max(a, b); }
/* GLSL step(edge, x) is 1 where x >= edge; floor(sign(x - e) * 0.5 + 1) is
that exactly, equality included, out of builtins Blink does have. */
template <class T, class S> T forma_step(S e, T x){ return floor(sign(x - e) * 0.5f + 1.0f); }
template <class T, class S> T forma_smoothstep(S a, S b, T x){
T t = forma_clamp((x - a) / (b - a), 0.0f, 1.0f);
return t * t * (3.0f - 2.0f * t);
}
template <class T> float forma_distance(T a, T b){ return length(a - b); }
float forma_tanh(float x){ float e = exp(2.0f * x); return (e - 1.0f) / (e + 1.0f); }
float forma_radians(float d){ return d * 0.01745329252f; }
// the page's hash2 is exact uint32; Blink has int, so the shifts are made
// logical by masking and the read-back is lifted into 0 .. 2^32
/* A uint read back as a float. Blink has no unsigned type, so a value past
2^31 arrives as a negative int and float() of it is negative. Measured on
gabor, whose own generator then returned uniforms in [-0.5, 0.5) and drew
a different picture — it compiled, it rendered, and only comparing it with
/* The order's ramp — the same cosine formulation the JS kit uses, so a
plate keeps its classification colour in either language. */
float3 ramp(float t){
return forma_clamp(u_pal_a + u_pal_b * cos(6.28318530718f * (u_pal_c * t + u_pal_d)), 0.0f, 1.0f);
}
/* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
float cycle(float t, float period){ return forma_fract(t / period + u_phase); }
float pingpong(float t, float period){
float u = cycle(t, period);
return u < 0.5f ? u * 2.0f : 2.0f - u * 2.0f;
}
/* Planck's law gives the spectrum; the CIE colour matching functions reduce a
spectrum to a colour. The tabulated functions are asserted-copyright and are
not in this tree — what is here is the published analytic fit of Wyman,
Sloan and Shirley (JCGT 2(2), 2013), evaluated from the coefficients of its
own Table 1.f
Every numeric constant below is declared inside the function that uses it
rather than at file scope. A composition emits two copies of a body and
namespaces the function declarations, so a file-scope const would be
declared twice by blackbody x blackbody and fail to link. */
/* Only the shape in wavelength matters — scaling a spectrum does not change
its chromaticity — so every constant factor is dropped. exp(-u) rather than
exp(u): at the cold end of the window u passes 80, and exp(u) there is
outside the range of a 32-bit float. */
float bbody_planck(float lambda, float T){
const float c2 = 1.438776877e7f; // hc/k_B, in nm K
float u = c2 / (lambda * T);
float e = exp(-u);
return u * u * u * u * u * e / (1.0f - e);
}
/* Divided by the largest value the curve reaches inside the visible band: the
Wien peak where that is visible, the nearer edge of the band where it is
not. c2/x0 is the Wien displacement constant, and x0 is the root of
x = 5(1 - exp(-x)) — Wien falls out of Planck by differentiating. */
float bbody_norm(float lambda, float T){
const float c2 = 1.438776877e7f;
const float x0 = 4.965114232f;
float peak = forma_clamp(c2 / (x0 * T), 380.0f, 780.0f);
return bbody_planck(lambda, T) / bbody_planck(peak, T);
}
/* One piecewise-continuous Gaussian lobe: the width differs either side of the
centre, which is how the fit follows an asymmetry a plain Gaussian cannot.
g and d are reciprocal widths. */
float bbody_lobe(float lambda, float a, float b, float g, float d){
float x = (lambda - b) * (lambda < b ? g : d);
return a * exp(-0.5f * x * x);
}
/* x-bar, y-bar, z-bar of the CIE 1931 two-degree observer: three lobes, then
two and two. The negative amplitude is the real dip in x-bar near 500 nm. */
float3 bbody_cmf(float lambda){
return float3(
bbody_lobe(lambda, 0.362f, 442.0f, 0.0624f, 0.0374f)
+ bbody_lobe(lambda, 1.056f, 599.8f, 0.0264f, 0.0323f)
+ bbody_lobe(lambda, -0.065f, 501.1f, 0.0490f, 0.0382f),
bbody_lobe(lambda, 0.821f, 568.8f, 0.0213f, 0.0247f)
+ bbody_lobe(lambda, 0.286f, 530.9f, 0.0613f, 0.0322f),
bbody_lobe(lambda, 1.217f, 437.0f, 0.0845f, 0.0278f)
+ bbody_lobe(lambda, 0.681f, 459.0f, 0.0385f, 0.0725f));
}
/* CIE XYZ to linear sRGB (IEC 61966-2-1 primaries at D65), then that
standard's transfer function — a linear toe under a shifted power curve. */
float3 bbody_rgb(float3 c){
return float3( 3.2404542f * c.x - 1.5371385f * c.y - 0.4985314f * c.z,
-0.9692660f * c.x + 1.8760108f * c.y + 0.0415560f * c.z,
0.0556434f * c.x - 0.2040259f * c.y + 1.0572252f * c.z);
}
float3 bbody_encode(float3 c){
c = forma_clamp(c, 0.0f, 1.0f);
return lerp(12.92f * c, 1.055f * pow(c, float3(1.0f / 2.4f)) - 0.055f, forma_step(0.0031308f, c));
}
/* The mean of the integrand over the visible band at one temperature — the
quantity both registers are built from. The mean rather than the sum, so
the sample count can change how accurate this is and never how bright. */
float3 bbody_sum(float T, float ns){
int n = int(floor(ns + 0.5f));
float3 xyz = float3(0.0f);
/* Constant loop bound at the slider's own ceiling, break on the live value,
mandelbrot fashion. Normalising the integrand by the column peak cancels
out of the chromaticity exactly and keeps the sum near 1 rather than near
1e-30, which is where a cold column would otherwise land. */
for (int i = 0; i < 96; i++){
if (i >= n) break;
float lambda = 380.0f + (float(i) + 0.5f) * 400.0f / float(n);
xyz += bbody_cmf(lambda) * bbody_norm(lambda, T);
}
return xyz / float(n);
}
float3 plate(float2 uv){
float3 forma_r = float3(0.0f, 0.0f, 0.0f);
for (int forma_once = 0; forma_once < 1; forma_once++){
const float FIELD = 0.60f;
const float BAND0 = 0.638f;
const float BAND1 = 0.940f;
const float RULE = 0.958f;
/* The largest value the register's ratio takes anywhere this plate can
reach. Measured, so exposure 1 is exactly the setting at which nothing
clips and anything above it is a decision. */
const float REG = 30.1f;
/* The window slides in log temperature and back, so the sweep keeps its
shape while it travels. cycle() carries this plate's PHASE. */
float slide = p_travel * sin(6.28318530718f * cycle(u_t, 34.0f));
float T = exp(lerp(log(p_tlo), log(p_thi), uv.x) + slide);
/* One integral per fragment, whichever register it lands in: the band needs
it for its colour, and the field needs its luminance to normalise by. */
float3 xyz = bbody_sum(T, p_samples);
float Y = forma_max(xyz.y, 1e-30);
if (uv.y < FIELD){
/* One term of the very sum below it: this wavelength's matching functions
taken straight through the sRGB matrix — the fit used only for what it
was fitted to, the absolute values — weighted by how much of this
wavelength the temperature emits. Negative components are simply what
it means for a single wavelength to be outside the gamut. */
float lambda = 380.0f + (uv.y / FIELD) * 400.0f;
float3 term = bbody_rgb(bbody_cmf(lambda)) * bbody_norm(lambda, T);
{ forma_r = bbody_encode(term / (Y * REG) * p_exposure); break; }
}
/* Chromaticity, not radiance, then scaled to the brightest sRGB triple that
carries it. A negative primary means sRGB cannot mix the colour at all. */
float3 lin = bbody_rgb(xyz / Y);
float m = forma_max(forma_max(lin.x, lin.y), lin.z);
if (uv.y > BAND0 && uv.y < BAND1) { forma_r = bbody_encode(forma_max(lin, 0.0f) / forma_max(m, 1e-9)); break; }
if (uv.y > RULE) { forma_r = lerp(float3(26.0f, 30.0f, 38.0f) / 255.0f, float3(209.0f / 255.0f),
forma_step(forma_min(forma_min(lin.x, lin.y), lin.z), -0.002f * m)); break; }
{ forma_r = float3(4.0f, 6.0f, 10.0f) / 255.0f; break; }
}
return forma_r;
}
void process(int2 pos){
// FORMA's uv runs y-down like its canvas; Nuke's rows run up
float2 uv = float2((float(pos.x) + 0.5f) / u_res.x, 1.0f - (float(pos.y) + 0.5f) / u_res.y);
float3 c = plate(uv);
dst() = float4(c.x, c.y, c.z, 1.0f);
}
};