PL. 119 · COLOUR / REFRACTION / SPECTRAL RAY FAN
Prism Dispersion
Ibn Sahl, 984 · Willebrord Snellius, c. 1621 · Isaac Newton, 1704 · Wilhelm Sellmeier, 1871
OPEN THE LIVE PLATE ▸DEFINITION
n²(λ) = 1 + Σᵢ Bᵢλ² / (λ² − Cᵢ) Sellmeier sin θ₁ = n(λ) sin θ₂ Snell, at both faces δ(λ) = θ₁ + θ₄ − A deviation through apex A
NOTES
White light meets a triangular prism and leaves it as a fan, because the refractive index of glass is not one number: it falls with wavelength, so violet is bent harder than red at the entry face and harder again at the exit face, and the two bends add. The index here is the Sellmeier formula carrying the six constants SCHOTT publishes for N-BK7, the commonest optical crown; each sampled wavelength is refracted in and out by the vector form of Snell, walked across the glass, and deposited along the ray it leaves on, in the colour that wavelength actually looks like. Two dials are the exhibit. The sample count is one: a spectrum is a continuum and this plate approximates it by a finite sum, so at four or six samples the fan is a few coloured lobes with gaps between them, and that banding is quadrature error rather than optics. Turn it up and watch the artefact leave. Watch also where it stops leaving — measured against a ninety-six sample reference, the picture is still moving at twenty and has stopped by about thirty, and that number is the answer to how finely a spectrum has to be cut before a screen cannot tell. The second dial is the exaggeration. A sixty degree crown prism spreads the whole visible band through two degrees, which is why Newton needed a darkened room and a wall a long way off, so at 1.0 — N-BK7 exactly as the data sheet publishes it — the spectrum here is a hairline. Above 1.0 the plate scales how far each index departs from the middle of the band while holding that middle fixed: a glass that does not exist, and the only way the effect fits on a card.
PROVENANCE
- Origin
- The law of refraction. Ibn Sahl set out the constant-ratio construction in his treatise on burning instruments, Baghdad, 984, which was identified and reconstructed by Roshdi Rashed, "A Pioneer in Anaclastics: Ibn Sahl on Burning Mirrors and Lenses", Isis 81(3), 1990, 464-491 (doi:10.1086/355456). Willebrord Snellius arrived at it again about 1621 and never published; it reached print in the Dioptrique of Rene Descartes, 1637, without the credit. Thomas Harriot had it earlier still, also unpublished. The name attached to the person who published least.
- The prism experiment
- Isaac Newton, Opticks, London, 1704, Book I — that the prism separates rays already present in white light rather than manufacturing colour, established by refracting one separated ray a second time and recombining the fan back to white. Pre-DOI; no resolvable identifier exists for it.
- Which dispersion formula, and why
- Two are in the literature and this plate runs the later one. Augustin-Louis Cauchy fitted the empirical n = A + B/λ² in his Memoire sur la dispersion de la lumiere, published at Prague inside the Nouveaux Exercices de Mathematiques, 1835-36; it is two constants and it is wrong near any absorption band, because it has no resonances in it. Wilhelm Sellmeier derived the resonant form from oscillators in "Zur Erklaerung der abnormen Farbenfolge im Spectrum einiger Substanzen", Annalen der Physik 219(6), 1871, 272-282 — the src link above — and that is the form optical glass is actually catalogued in, so it is the one a plate claiming to run real glass has to run. Cauchy is named here because he is the honest simpler exhibit and because the plate would be lying by omission if it let Sellmeier look like the only answer.
- The glass
- The six constants are the ones SCHOTT prints on the N-BK7 data sheet, revision 19 August 2010, under Constants of Dispersion Formula: B = 1.03961212, 0.231792344, 1.010469450 and C = 0.00600069867, 0.0200179144, 103.560653, with C in square micrometres. They are manufacturer measurements of one glass, not a derivation, and the plate says so rather than presenting them as physics. Checked rather than copied: evaluated here at the helium d line, 587.56 nm, the formula returns 1.516800 against the 1.51680 the same sheet prints as n_d; across the hydrogen F and C lines it returns 0.008054 against the printed 0.008054; and the Abbe number those imply is 64.17 against the printed 64.17. Three independent figures on the same sheet, all reproduced, which is what tells you six constants were transcribed without a slip.
- Wavelength into colour
- The CIE 1931 two-degree observer, evaluated through the multi-lobe Gaussian fit of Chris Wyman, Peter-Pike Sloan and Peter Shirley, "Simple Analytic Approximations to the CIE XYZ Color Matching Functions", Journal of Computer Graphics Techniques 2(2), 2013, jcgt.org/published/0002/02/01 — seven lobes, each a Gaussian whose width differs either side of its own peak. The fit is used because the CIE tabulations are asserted-copyright and never enter this tree; the seven amplitudes, centres and reciprocal widths are read from equation 4 and table 1 of that paper and evaluated here from those numbers, not from the code the paper lists beside them.
- What a screen cannot show
- A single wavelength lies outside the sRGB gamut almost everywhere, so converting XYZ to display primaries returns negative components, and they are clipped at zero. The ends of the fan therefore read less saturated than they are, and the deep violet reads dark because the observer itself is nearly blind there. The plate clips and says so rather than pretending a three-primary screen can show a spectral line.
- Standing
- Public domain. Every step is implemented from the published equations named above — Sellmeier from its resonant form, Snell in the vector form the tangential boundary condition gives directly, the observer from equation 4 of the fit. Nothing is taken from any optics or rendering codebase.
- Constants
- samples is locked from regenerate as a cost dial and stays live for the reader, because moving it is the point. apex, incidence and exaggeration are all jitterable, and all 320 random tuples drawn from the whole declared box render live — which is only true because the incidence dial names a position inside the transmitting window rather than an angle. Below a threshold incidence the beam never leaves the second face at all: the threshold is asin(n sin(A − asin(1/n))), measured here at 30.5 degrees for a 60 degree apex in N-BK7, 36.1 degrees at exaggeration 6, and 52.4 degrees at a 64 degree apex and exaggeration 12, with no threshold at all below about a 40.7 degree apex. A dial in degrees therefore has a dead lower half whose size depends on two other sliders: over 200,000 random tuples of a plain apex-and-incidence box, 26.7 per cent came out trapped and blank. The window is also computed from the most refractive sample, 380 nm, so that every wavelength transmits at every setting rather than the violet end vanishing first.
- Source
- doi:10.1002/andp.18712190612
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. 119 · PRISM DISPERSION — Ibn Sahl, 984 · Willebrord Snellius, c. 1621 · Isaac Newton, 1704 · Wilhelm Sellmeier, 1871
// n²(λ) = 1 + Σᵢ Bᵢλ² / (λ² − Cᵢ) Sellmeier
// sin θ₁ = n(λ) sin θ₂ Snell, at both faces
// δ(λ) = θ₁ + θ₄ − A deviation through apex A
// 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.5535; // 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_samples = 24.0; // spectral samples · live 4 .. 40
const float p_apex = 60.0; // prism apex angle (deg) · live 38 .. 64
const float p_incid = 0.22; // incidence within the transmitting window · live 0.1 .. 0.82
const float p_disp = 6.0; // dispersion exaggeration (1 = N-BK7) · live 1 .. 12
/* 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;
}
/* The same field, evaluated at the full backing resolution. Everything the
JS path computes once per frame is recomputed per pixel here, because a
fragment shader has nowhere to put it — the numbers are identical, so the
two paths are one picture sampled twice. */
/* Sellmeier with SCHOTT's published N-BK7 constants; lambda in micrometres,
because the C terms are in square micrometres. */
float disp_index(float um){
float l2 = um * um;
return sqrt(1.0 + 1.03961212 * l2 / (l2 - 0.00600069867)
+ 0.231792344 * l2 / (l2 - 0.0200179144)
+ 1.010469450 * l2 / (l2 - 103.560653));
}
/* p_disp scales the departure from the middle of the band and leaves the
middle fixed, so the fan opens without the whole beam moving. */
float disp_glass(float nm){
float mid = disp_index(0.540);
return mid + p_disp * (disp_index(nm * 0.001) - mid);
}
/* One lobe of the Wyman-Sloan-Shirley fit: a Gaussian whose width differs
either side of its own peak. The tabulated widths are reciprocals, so
they multiply. */
float disp_lobe(float nm, float amp, float mid, float wLo, float wHi){
float u = (nm - mid) * (nm < mid ? wLo : wHi);
return amp * exp(-0.5 * u * u);
}
/* CIE 1931 two-degree observer, then linear sRGB, clipped at zero — a
spectral line is out of gamut nearly everywhere and the negative
components are real. */
vec3 disp_light(float nm){
float X = disp_lobe(nm, 0.362, 442.0, 0.0624, 0.0374)
+ disp_lobe(nm, 1.056, 599.8, 0.0264, 0.0323)
+ disp_lobe(nm, -0.065, 501.1, 0.0490, 0.0382);
float Y = disp_lobe(nm, 0.821, 568.8, 0.0213, 0.0247)
+ disp_lobe(nm, 0.286, 530.9, 0.0613, 0.0322);
float Z = disp_lobe(nm, 1.217, 437.0, 0.0845, 0.0278)
+ disp_lobe(nm, 0.681, 459.0, 0.0385, 0.0725);
return max(vec3( 3.2404542 * X - 1.5371385 * Y - 0.4985314 * Z,
-0.9692660 * X + 1.8760108 * Y + 0.0415560 * Z,
0.0556434 * X - 0.2040259 * Y + 1.0572252 * Z), 0.0);
}
vec3 disp_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));
}
/* Snell in vector form. N is the unit normal on the side the ray arrives
from; the third component is 0 when the discriminant goes negative, which
is total internal reflection. */
vec3 disp_refract(vec2 I, vec2 N, float eta){
float ci = -dot(I, N);
float k = 1.0 - eta * eta * (1.0 - ci * ci);
if (k < 0.0) return vec3(0.0);
return vec3(eta * I + (eta * ci - sqrt(k)) * N, 1.0);
}
vec2 disp_turn(vec2 v, float c, float s){ return vec2(v.x * c - v.y * s, v.x * s + v.y * c); }
float disp_seg(vec2 q, vec2 a, vec2 b){
vec2 e = b - a, r = q - a;
float L2 = max(dot(e, e), 1e-12);
vec2 d = r - e * clamp(dot(r, e) / L2, 0.0, 1.0);
return dot(d, d);
}
float disp_glow(float d2, float s2){ float a = s2 / (s2 + d2); return a * a; }
/* One wavelength through the glass, in units of one prism side: refract in,
take the nearest of the three faces the internal ray is heading for,
refract out. xy is the exit point, zw the exit direction — zero when that
face reflected it instead. */
vec4 disp_trace(float nm, float hf, vec2 entry, vec2 inc){
float sh = sin(hf), ch = cos(hf);
float n = disp_glass(nm);
vec3 din = disp_refract(inc, vec2(-ch, -sh), 1.0 / n);
float tm = 1e9;
int face = -1;
vec2 hitN = vec2(0.0);
for (int e = 0; e < 3; e++){
vec2 fn = e == 0 ? vec2(-ch, -sh) : (e == 1 ? vec2(ch, -sh) : vec2(0.0, 1.0));
float fd = e == 2 ? ch : 0.0;
float den = dot(din.xy, fn);
if (den > 1e-9){
float tt = (fd - dot(entry, fn)) / den;
if (tt > 1e-5 && tt < tm){ tm = tt; face = e; hitN = fn; }
}
}
if (face < 0) return vec4(entry, 0.0, 0.0);
vec2 q = entry + tm * din.xy;
vec3 o = disp_refract(din.xy, -hitN, n);
return vec4(q, o.xy * o.z);
}
vec3 plate(vec2 uv){
float asp = u_res.y / u_res.x;
vec2 q = vec2(uv.x, uv.y * asp);
float hf = radians(p_apex) * 0.5;
float sh = sin(hf), ch = cos(hf);
vec2 entry = vec2(-0.45 * sh, 0.45 * ch);
/* The smallest incidence that still lets the most refracted sample out of
the second face. The dial names a position above it, which is what
keeps the whole declared box alive. */
float nv = disp_glass(380.0);
float crit = asin(1.0 / nv);
float apex = radians(p_apex);
float lo = apex > crit ? asin(min(1.0, nv * sin(apex - crit))) : 0.0;
float hi = max(lo + radians(4.0), radians(80.0));
float u = clamp(p_incid + 0.12 * (2.0 * pingpong(u_t, 34.0) - 1.0), 0.01, 0.99);
float th = lo + u * (hi - lo);
vec2 inc = vec2(cos(th - hf), -sin(th - hf));
/* The assembly turns so the mid-band ray always leaves on the same
bearing; what moves on screen is the prism. */
vec4 mid = disp_trace(540.0, hf, entry, inc);
vec2 md = dot(mid.zw, mid.zw) > 0.5 ? mid.zw : vec2(1.0, 0.0);
float phi = radians(18.0) - atan(md.y, md.x);
float cp = cos(phi), sp = sin(phi);
float side = 0.30 * min(1.0, asp);
vec2 org = vec2(0.30, 0.45 * asp) - side * disp_turn(vec2(0.0, 2.0 * ch / 3.0), cp, sp);
vec2 p1 = org + side * disp_turn(entry, cp, sp);
float SIG2 = 0.0055 * 0.0055;
/* The internal segments never leave the glass, so a pixel clear of the
prism cannot be near one. The JS path carries the same test with the
same twelve-half-width margin, where the kernel is under 5e-5 of peak. */
float circum = max(2.0 * ch / 3.0, sqrt(sh * sh + ch * ch / 9.0)) * side + 12.0 * 0.0055;
vec2 g = q - vec2(0.30, 0.45 * asp);
bool nearGlass = dot(g, g) < circum * circum;
int ns = max(1, int(p_samples + 0.5)); // the count is a divisor
vec3 acc = vec3(0.0), tot = vec3(0.0);
/* Constant bound at the slider's own ceiling, break on the live value. */
for (int k = 0; k < 40; k++){
if (k >= ns) break;
float nm = 380.0 + 320.0 * (float(k) + 0.5) / float(ns);
vec3 rgb = disp_light(nm);
tot += rgb;
vec4 r = disp_trace(nm, hf, entry, inc);
vec2 p2 = org + side * disp_turn(r.xy, cp, sp);
float w = nearGlass ? disp_glow(disp_seg(q, p1, p2), SIG2) : 0.0;
if (dot(r.zw, r.zw) > 0.5)
w += disp_glow(disp_seg(q, p2, p2 + disp_turn(r.zw, cp, sp) * 3.0), SIG2);
acc += rgb * w;
}
tot = max(tot, vec3(1e-4));
/* The beam going in is every sample at once, so it is laid down with the
sum the fan was divided by and comes out white by construction. */
vec2 back = p1 - disp_turn(inc, cp, sp) * 0.34;
acc += tot * disp_glow(disp_seg(q, back, p1), SIG2);
vec3 lit = disp_encode(3.2 * acc / tot);
/* Signed distance to the glass, read in the prism's own frame: for a
convex body the largest of the three plane distances is negative inside
and is the distance to the nearest face. */
vec2 loc = disp_turn(q - org, cp, -sp) / side;
float sd = max(max(dot(loc, vec2(-ch, -sh)), dot(loc, vec2(ch, -sh))), loc.y - ch) * side;
vec3 edge = ramp(0.95);
lit += edge * (0.55 * disp_glow(sd * sd, 0.0022 * 0.0022) + (sd < 0.0 ? 0.045 : 0.0));
return vec3(4.0, 6.0, 10.0) / 255.0 + lit;
}
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. 119 · PRISM DISPERSION — Ibn Sahl, 984 · Willebrord Snellius, c. 1621 · Isaac Newton, 1704 · Wilhelm Sellmeier, 1871
// n²(λ) = 1 + Σᵢ Bᵢλ² / (λ² − Cᵢ) Sellmeier
// sin θ₁ = n(λ) sin θ₂ Snell, at both faces
// δ(λ) = θ₁ + θ₄ − A deviation through apex A
// 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_dispersion : ImageComputationKernel<ePixelWise>
{
Image<eWrite> dst;
param:
float u_t; // seconds; 0 is the still frame
float p_samples; // spectral samples · live 4 .. 40
float p_apex; // prism apex angle (deg) · live 38 .. 64
float p_incid; // incidence within the transmitting window · live 0.1 .. 0.82
float p_disp; // dispersion exaggeration (1 = N-BK7) · live 1 .. 12
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_samples, "p_samples", 24.0f);
defineParam(p_apex, "p_apex", 60.0f);
defineParam(p_incid, "p_incid", 0.22f);
defineParam(p_disp, "p_disp", 6.0f);
}
void init(){
u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
u_phase = 0.5535f; // 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;
}
/* The same field, evaluated at the full backing resolution. Everything the
JS path computes once per frame is recomputed per pixel here, because a
fragment shader has nowhere to put it — the numbers are identical, so the
two paths are one picture sampled twice. */
/* Sellmeier with SCHOTT's published N-BK7 constants; lambda in micrometres,
because the C terms are in square micrometres. */
float disp_index(float um){
float l2 = um * um;
return sqrt(1.0f + 1.03961212f * l2 / (l2 - 0.00600069867f)
+ 0.231792344f * l2 / (l2 - 0.0200179144f)
+ 1.010469450f * l2 / (l2 - 103.560653f));
}
/* p_disp scales the departure from the middle of the band and leaves the
middle fixed, so the fan opens without the whole beam moving. */
float disp_glass(float nm){
float mid = disp_index(0.540f);
return mid + p_disp * (disp_index(nm * 0.001f) - mid);
}
/* One lobe of the Wyman-Sloan-Shirley fit: a Gaussian whose width differs
either side of its own peak. The tabulated widths are reciprocals, so
they multiply. */
float disp_lobe(float nm, float amp, float mid, float wLo, float wHi){
float u = (nm - mid) * (nm < mid ? wLo : wHi);
return amp * exp(-0.5f * u * u);
}
/* CIE 1931 two-degree observer, then linear sRGB, clipped at zero — a
spectral line is out of gamut nearly everywhere and the negative
components are real. */
float3 disp_light(float nm){
float X = disp_lobe(nm, 0.362f, 442.0f, 0.0624f, 0.0374f)
+ disp_lobe(nm, 1.056f, 599.8f, 0.0264f, 0.0323f)
+ disp_lobe(nm, -0.065f, 501.1f, 0.0490f, 0.0382f);
float Y = disp_lobe(nm, 0.821f, 568.8f, 0.0213f, 0.0247f)
+ disp_lobe(nm, 0.286f, 530.9f, 0.0613f, 0.0322f);
float Z = disp_lobe(nm, 1.217f, 437.0f, 0.0845f, 0.0278f)
+ disp_lobe(nm, 0.681f, 459.0f, 0.0385f, 0.0725f);
return forma_max(float3( 3.2404542f * X - 1.5371385f * Y - 0.4985314f * Z,
-0.9692660f * X + 1.8760108f * Y + 0.0415560f * Z,
0.0556434f * X - 0.2040259f * Y + 1.0572252f * Z), 0.0f);
}
float3 disp_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));
}
/* Snell in vector form. N is the unit normal on the side the ray arrives
from; the third component is 0 when the discriminant goes negative, which
is total internal reflection. */
float3 disp_refract(float2 I, float2 N, float eta){
float3 forma_r = float3(0.0f, 0.0f, 0.0f);
for (int forma_once = 0; forma_once < 1; forma_once++){
float ci = -dot(I, N);
float k = 1.0f - eta * eta * (1.0f - ci * ci);
if (k < 0.0f) { forma_r = float3(0.0f); break; }
{ forma_r = float3(eta * I + (eta * ci - sqrt(k)) * N, 1.0f); break; }
}
return forma_r;
}
float2 disp_turn(float2 v, float c, float s){ return float2(v.x * c - v.y * s, v.x * s + v.y * c); }
float disp_seg(float2 q, float2 a, float2 b){
float2 e = b - a;
float2 r = q - a;
float L2 = forma_max(dot(e, e), 1e-12);
float2 d = r - e * forma_clamp(dot(r, e) / L2, 0.0f, 1.0f);
return dot(d, d);
}
float disp_glow(float d2, float s2){ float a = s2 / (s2 + d2); return a * a; }
/* One wavelength through the glass, in units of one prism side: refract in,
take the nearest of the three faces the internal ray is heading for,
refract out. xy is the exit point, zw the exit direction — zero when that
face reflected it instead. */
float4 disp_trace(float nm, float hf, float2 entry, float2 inc){
float4 forma_r = float4(0.0f, 0.0f, 0.0f, 0.0f);
for (int forma_once = 0; forma_once < 1; forma_once++){
float sh = sin(hf);
float ch = cos(hf);
float n = disp_glass(nm);
float3 din = disp_refract(inc, float2(-ch, -sh), 1.0f / n);
float tm = 1e9;
int face = -1;
float2 hitN = float2(0.0f);
for (int e = 0; e < 3; e++){
float2 fn = e == 0 ? float2(-ch, -sh) : (e == 1 ? float2(ch, -sh) : float2(0.0f, 1.0f));
float fd = e == 2 ? ch : 0.0f;
float den = dot(float2(din.x, din.y), fn);
if (den > 1e-9){
float tt = (fd - dot(entry, fn)) / den;
if (tt > 1e-5 && tt < tm){ tm = tt; face = e; hitN = fn; }
}
}
if (face < 0) { forma_r = float4(entry, 0.0f, 0.0f); break; }
float2 q = entry + tm * float2(din.x, din.y);
float3 o = disp_refract(float2(din.x, din.y), -hitN, n);
{ forma_r = float4(q, float2(o.x, o.y) * o.z); break; }
}
return forma_r;
}
float3 plate(float2 uv){
float asp = u_res.y / u_res.x;
float2 q = float2(uv.x, uv.y * asp);
float hf = forma_radians(p_apex) * 0.5f;
float sh = sin(hf);
float ch = cos(hf);
float2 entry = float2(-0.45f * sh, 0.45f * ch);
/* The smallest incidence that still lets the most refracted sample out of
the second face. The dial names a position above it, which is what
keeps the whole declared box alive. */
float nv = disp_glass(380.0f);
float crit = asin(1.0f / nv);
float apex = forma_radians(p_apex);
float lo = apex > crit ? asin(forma_min(1.0f, nv * sin(apex - crit))) : 0.0f;
float hi = forma_max(lo + forma_radians(4.0f), forma_radians(80.0f));
float u = forma_clamp(p_incid + 0.12f * (2.0f * pingpong(u_t, 34.0f) - 1.0f), 0.01f, 0.99f);
float th = lo + u * (hi - lo);
float2 inc = float2(cos(th - hf), -sin(th - hf));
/* The assembly turns so the mid-band ray always leaves on the same
bearing; what moves on screen is the prism. */
float4 mid = disp_trace(540.0f, hf, entry, inc);
float2 md = dot(float2(mid.z, mid.w), float2(mid.z, mid.w)) > 0.5f ? float2(mid.z, mid.w) : float2(1.0f, 0.0f);
float phi = forma_radians(18.0f) - atan2(md.y, md.x);
float cp = cos(phi);
float sp = sin(phi);
float side = 0.30f * forma_min(1.0f, asp);
float2 org = float2(0.30f, 0.45f * asp) - side * disp_turn(float2(0.0f, 2.0f * ch / 3.0f), cp, sp);
float2 p1 = org + side * disp_turn(entry, cp, sp);
float SIG2 = 0.0055f * 0.0055f;
/* The internal segments never leave the glass, so a pixel clear of the
prism cannot be near one. The JS path carries the same test with the
same twelve-half-width margin, where the kernel is under 5e-5 of peak. */
float circum = forma_max(2.0f * ch / 3.0f, sqrt(sh * sh + ch * ch / 9.0f)) * side + 12.0f * 0.0055f;
float2 g = q - float2(0.30f, 0.45f * asp);
bool nearGlass = dot(g, g) < circum * circum;
int ns = forma_max(1, int(p_samples + 0.5f)); // the count is a divisor
float3 acc = float3(0.0f);
float3 tot = float3(0.0f);
/* Constant bound at the slider's own ceiling, break on the live value. */
for (int k = 0; k < 40; k++){
if (k >= ns) break;
float nm = 380.0f + 320.0f * (float(k) + 0.5f) / float(ns);
float3 rgb = disp_light(nm);
tot += rgb;
float4 r = disp_trace(nm, hf, entry, inc);
float2 p2 = org + side * disp_turn(float2(r.x, r.y), cp, sp);
float w = nearGlass ? disp_glow(disp_seg(q, p1, p2), SIG2) : 0.0f;
if (dot(float2(r.z, r.w), float2(r.z, r.w)) > 0.5f)
w += disp_glow(disp_seg(q, p2, p2 + disp_turn(float2(r.z, r.w), cp, sp) * 3.0f), SIG2);
acc += rgb * w;
}
tot = forma_max(tot, float3(1e-4));
/* The beam going in is every sample at once, so it is laid down with the
sum the fan was divided by and comes out white by construction. */
float2 back = p1 - disp_turn(inc, cp, sp) * 0.34f;
acc += tot * disp_glow(disp_seg(q, back, p1), SIG2);
float3 lit = disp_encode(3.2f * acc / tot);
/* Signed distance to the glass, read in the prism's own frame: for a
convex body the largest of the three plane distances is negative inside
and is the distance to the nearest face. */
float2 loc = disp_turn(q - org, cp, -sp) / side;
float sd = forma_max(forma_max(dot(loc, float2(-ch, -sh)), dot(loc, float2(ch, -sh))), loc.y - ch) * side;
float3 edge = ramp(0.95f);
lit += edge * (0.55f * disp_glow(sd * sd, 0.0022f * 0.0022f) + (sd < 0.0f ? 0.045f : 0.0f));
return float3(4.0f, 6.0f, 10.0f) / 255.0f + lit;
}
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);
}
};