I still have a ways to go with more explanatory diagrams but I added an amusing picture to the top of

… so it’s not just formulas anymore. Warning: moving the points around and watching the little rainbow lines move in response is kind of addicting, and might distract you from whatever work you were otherwise trying to do.

@john-clarke

Try this self-contained version:

```
// Given an array of knots and an array of values,
// return a function which evaluates a natural cubic spline
// extrapolated by linear functions on either end.
NCS = function NCS(knots, values) {
const n = knots.length - 1;
const t = new Float64Array(n+3); t.set(knots, 1);
const y = new Float64Array(n+3); y.set(values, 1);
const s = new Float64Array(n+2); // horizontal scales
const m = new Float64Array(n+3); // slope at each knot
const d = new Float64Array(n+3); // diagonal matrix
// Natural cubic spline algorithm
for (let i = 1; i < n+1; i++) {
s[i] = 1 / (t[i+1] - t[i]);
m[i] += (m[i+1] = 3*s[i]*s[i]*(y[i+1] - y[i])); }
d[1] = 0.5 / s[1];
m[1] = d[1] * m[1];
for (let i = 2; i <= n+1; i++) {
d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);
m[i] = d[i]*(m[i] - s[i-1]*m[i-1]); }
for (let i = n; i >= 1; i--) {
m[i] -= d[i]*s[i]*m[i+1]; }
// Linear extrapolation
t[0] = t[1] - 1; t[n+2] = t[n+1] + 1;
y[0] = y[1] - m[1]; y[n+2] = y[n+1] + m[n+1];
s[0] = s[n+1] = 1;
m[0] = m[1]; m[n+2] = m[n+1];
// Set up Chebyshev coefficients
const coeffs = new Float64Array(4*(n+2));
for (let i = 0; i < n+2; i++) {
const h = t[i+1] - t[i];
const y0 = y[i], y1 = y[i+1], m0 = h*m[i], m1 = h*m[i+1], j = 4*i;
coeffs[j] = 0.5*(y0 + y1 + 0.125*(m0 - m1));
coeffs[j+1] = 0.03125*(18*(y1 - y0) - m0 - m1);
coeffs[j+2] = 0.0625*(m1 - m0);
coeffs[j+3] = 0.03125*(2*(y0 - y1) + m0 + m1); }
return Object.assign(
function eval_spline(u) {
// Binary search
let nn = t.length - 1, low = 0, half;
while ((half = (nn/2)|0) > 0) {
low += half * (t[low + half] <= u);
nn -= half; }
const i = low, j = 4*i;
u = (2*u - t[i] - t[i+1]) * s[i]; // scale to [–1, 1]
// Clenshaw's method.
let b1 = coeffs[j+3], b2 = coeffs[j+2] + 2 * u * b1;
b1 = coeffs[j+1] + 2 * u * b2 - b1;
return coeffs[j] + u * b1 - b2; },
{ knots: t, values: y, dvalues: m, scales: s, coeffs: coeffs });
}
```

Example usage:

```
{
const
t = [1, 2.5, 3, 4, 5.5, 6, 7, 7.2, 9],
y = [1, 1.5, 1, 2, 1, 1.5, 1.3, 1.1, 1],
spline = NCS(t, y),
x = [...Array(500).keys()].map(i=> i/50),
div = DOM.element('div');
Plotly.newPlot(div, [
{ x: x, y: x.map(spline), mode: 'lines', type: 'scatter'},
{ x: t, y: y, mode: 'markers', type: 'scatter'}
], {width: width});
return div;
}
```