30:T460,Import the Zernike calculator using: ``` import {calculateZernike} from "@dkirkby/zernike" ``` then, for example, to calculate and print all results up to a max radial degree of 4: ``` const x = [-0.5, 0], y = [0.25, -0.25]; const {Z} = calculateZernike(4, x, y); for(let n = 0; n <= 4; ++n) { for(let l = -n; l <= n; l+=2) { for(let k = 0; k < x.length; ++k) { console.log(`Z[n=\${n},l=\${l}](x=\${x[k]},y=\${y[k]}) = \${Z(n,l,k)}`); } } } ``` If you need the polar coordinates calculated for your input Cartesian coordinates use: ``` const {Z, r, phi} = calculateZernike(4, x, y); ``` To obtain functions that evaluate the partial derivatives ${tex`\partial Z/\partial x`} and ${tex`\partial Z/\partial y`} use, for example: ``` const {Z, dZdx, dZdy} = calculateZernike(4, x, y); ``` You might find the convenience function `diskXYGrid` useful to initialize a rectangular grid of points on the unit disk at some resolution: ``` import {calculateZernike, diskXYGrid} from "@dkirkby/zernike" const [x, y] = diskXYGrid(10); // returns arrays of length ~ 10**2 (PI/4) const {Z} = getZernike(4, x, y); ```31:T14ab,function calculateZernike(maxRadialDegree, x, y, { debug=false, trace=false, profile=false, arrayType=Float64Array, outerValue=0}={}) { const start = profile ? performance.now() : 0; const {sqrt, atan2, floor, abs, cos, sin, pow} = Math; if(maxRadialDegree < 2) throw new Error("Expected maxRadialDegree >= 2."); if(x.length != y.length) throw new Error("Input arrays x and y must have the same length."); const nxy = x.length; // Calculate polar coordinates (r,phi) for each (x,y). const r = arrayType.from({length: nxy}, (_,i) => sqrt(x[i]**2 + y[i]**2)); const phi = arrayType.from({length: nxy}, (_,i) => atan2(y[i], x[i])); // Calculate C[m] = r**m cos(m phi) and S[m] = r**m sin(m phi) with recurrence. const N = maxRadialDegree + 1; const C = new arrayType(N * nxy), S = new arrayType(N * nxy); C[0] = 1; S[0] = 0; let X, Y, j; for(let k = 0; k < nxy; ++k) { if(r[k] > 1) continue; j = N * k; C[j] = 1, S[j] = 0; X = x[k], Y = y[k]; for(let n = 1; n < N; ++n) { j++; C[j] = X * C[j-1] - Y * S[j-1]; S[j] = Y * C[j-1] + X * S[j-1]; } } const getC = (l, k) => C[k * N + l]; const getS = (l, k) => S[k * N + l]; if(debug) { // Check the recurrence calculation against explicit const eps = 1e-6; for(let l = 0; l < N; ++l) { for(let k = 0; k < nxy; ++k) { const C2 = pow(r[k], l) * cos(l * phi[k]); const S2 = pow(r[k], l) * sin(l * phi[k]); if((abs(C2 - getC(l,k)) > eps) || (abs(S2 - getS(l,k)) > eps)) throw new Error(`C,S recurrence check failed for (l,k) = (${l},${k}).`); } } } // Calculate A[m,l,k] = A[m,l](r[k]) with recurrence (with l = |n - 2m|) const ncells = (i) => i % 2 == 0 ? (i / 2) * (i / 2 + 1) : (i + 1) * (i + 1) / 4; const M = ncells(N); const A = new arrayType(M * nxy).fill(NaN); const getIndex = (n, l, k) => k * M + ncells(n) + floor(l / 2); const setAdebug = (n, l, k, a) => { const idx = getIndex(n, l, k); if(trace) console.log(`setA(n=${n},l=${l},k=${k},idx=${idx}) = ${a}`); if(!Number.isNaN(A[idx])) throw new Error( `Recurrence overwrites (n,l,k) = (${n},${l},${k}) at index ${idx}.`); A[idx] = a; }; const setA = debug ? setAdebug : ((n, l, k, a) => { A[getIndex(n, l, k)] = a; }); const getAdebug = (n, l, k) => { const idx = getIndex(n, l, k) const a = A[idx]; if(trace) console.log(` getA(n=${n},l=${l},k=${k},idx=${idx}) = ${a}`); if(Number.isNaN(a)) throw new Error( `Recurrence uses un-initialized (n,l,k) = (${n},${l},${k}) at index ${idx}.`); return a; } const getA = debug ? getAdebug : ((n, l, k) => A[getIndex(n, l, k)]); for(let k = 0; k < nxy; ++k) { if(r[k] > 1) continue; const R2 = r[k] ** 2; setA(0, 0, k, 1); setA(1, 1, k, 1); for(let n = 2; n < N; ++n) { if(n % 2 == 0) { setA(n, 0, k, 2 * R2 * getA(n-1, 1, k) - getA(n-2, 0, k)); for(let l = 2; l < n - 1; l += 2) { setA(n, l, k, getA(n-1, l-1, k) + R2 * getA(n-1, l+1, k) - getA(n-2, l, k)); } setA(n, n, k, getA(n-1, n-1, k)); } else { for(let l = 1; l < n - 1; l += 2) { setA(n, l, k, getA(n-1, l-1, k) + R2 * getA(n-1, l+1, k) - getA(n-2, l, k)); } setA(n, n, k, getA(n-1, n-1, k)); } } } // Note that we exchange cos/sin relative to eqn (25) of https://doi.org/10.1117/12.2020389 // to match the convention in the wikipedia article. const Z = (n, l, k) => { if(r[k] > 1) return outerValue; if(l > 0) return sqrt(2 * (n + 1)) * getA(n,l,k) * getC(l,k); if(l < 0) return sqrt(2 * (n + 1)) * getA(n,-l,k) * getS(-l,k); return sqrt(n + 1) * getA(n,0,k) * getC(0,k); } // Implement Cartesian derivatives using eqns (23,24) from (but they use m for our l) // https://wp.optics.arizona.edu/visualopticslab/wp-content/uploads/sites/52/2016/08/Zernike-Notes-15Jan2016.pdf const dZdx = (n, l, k) => { if((r[k] > 1) || (n == 0)) return 0; let np, lp, norm, result = 0; const absl = abs(l), sign = l == 0 ? 1 : l / absl; for(np = absl+1; np <= n - 1; np += 2) { lp = sign * (absl + 1); result += sqrt(np + 1) * Z(np, lp, k); } if((l != 0) && (l != -1)) { norm = (l == 1) ? sqrt(2) : 1; for(np = absl-1; np <= n - 1; np += 2) { lp = sign * (absl - 1); result += norm * sqrt(np + 1) * Z(np, lp, k); } } norm = (l == 0) ? sqrt(2*(n + 1)) : sqrt(n + 1); return norm * result; } const dZdy = (n, l, k) => { if((r[k] > 1) || (n == 0)) return 0; let np, lp, norm, result = 0; const absl = abs(l), sign = l == 0 ? 1 : l / absl; for(np = absl+1; np <= n - 1; np += 2) { lp = sign * (absl + 1); result += sqrt(np + 1) * Z(np, -lp, k); } if((l != 0) && (l != 1)) { norm = (l == -1) ? sqrt(2) : 1; for(np = absl-1; np <= n - 1; np += 2) { lp = sign * (absl - 1); result -= norm * sqrt(np + 1) * Z(np, -lp, k); } } norm = (l == 0) ? sqrt(2*(n + 1)) : sqrt(n + 1); return norm * sign * result; } if(profile) { const stop = performance.now(); const elapsed = stop - start; return {Z, r, phi, elapsed}; } else return {Z, dZdx, dZdy, r, phi}; }32:T5b3,namedZernikes = { const {sin,cos,sqrt} = Math; return [ { n:0, l: 0, j: 1, Zj: (r,phi) => 1, name: "piston" }, { n:1, l:-1, j: 3, Zj: (r,phi) => 2*r*sin(phi), name: "vertical tilt" }, { n:1, l:+1, j: 2, Zj: (r,phi) => 2*r*cos(phi), name: "horizontal tilt" }, { n:2, l:-2, j: 5, Zj: (r,phi) => sqrt(6)*r*r*sin(2*phi), name: "oblique astigmatism" }, { n:2, l: 0, j: 4, Zj: (r,phi) => sqrt(3)*(2*r*r-1), name: "defocus" }, { n:2, l:+2, j: 6, Zj: (r,phi) => sqrt(6)*r*r*cos(2*phi), name: "vertical astigmatism" }, { n:3, l:-3, j: 9, Zj: (r,phi) => sqrt(8)*r*r*r*sin(3*phi), name: "vertical trefoil" }, { n:3, l:-1, j: 7, Zj: (r,phi) => sqrt(8)*(3*r*r*r-2*r)*sin(phi), name: "vertical coma" }, { n:3, l:+1, j: 8, Zj: (r,phi) => sqrt(8)*(3*r*r*r-2*r)*cos(phi), name: "horizontal coma" }, { n:3, l:+3, j:10, Zj: (r,phi) => sqrt(8)*r*r*r*cos(3*phi), name: "oblique trefoil" }, { n: 4, l:-4, j:15, Zj: (r,phi) => sqrt(10)*r*r*r*r*sin(4*phi), name: "oblique quadrafoil" }, { n: 4, l:-2, j:13, Zj: (r,phi) => sqrt(10)*(4*r*r*r*r-3*r*r)*sin(2*phi), name: "oblique secondary astigmatism" }, { n: 4, l: 0, j:11, Zj: (r,phi) => sqrt(5)*(6*r*r*r*r-6*r*r+1), name: "primary spherical" }, { n: 4, l:+2, j:12, Zj: (r,phi) => sqrt(10)*(4*r*r*r*r-3*r*r)*cos(2*phi), name: "vertical secondary astigmatism" }, { n: 4, l:+4, j:14, Zj: (r,phi) => sqrt(10)*r*r*r*r*cos(4*phi), name: "vertical quadrafoil" }, ]; }33:Tbb6,function plotZernike(maxRadialDegree, { plotSize=0.5*width, cellSize=1, numContours=21, }={}) { const ncell = Math.round(plotSize / cellSize); const xygrid = d3.range(ncell).map((_,i) => (i - (ncell - 1) / 2) / ((ncell - 1) / 2)); const [y, x] = d3.transpose(d3.cross(xygrid, xygrid)); const calculated = calculateZernike(maxRadialDegree, x, y); // Initialize transforms from data index space [0,ncell) x [0,ncell) to pixel space. const X = d3.scaleLinear().domain([0, ncell-1]).range([0, plotSize-1]); const Y = d3.scaleLinear().domain([0, ncell-1]).range([plotSize-1, 0]); const transform = ({type, value, coordinates}) => { return {type, value, coordinates: coordinates.map(rings => { return rings.map(points => { return points.map(([ix, iy]) => [ X(ix), Y(iy) ] ); }); })}; }; // Initialize the SVG output. const svg = d3.create("svg") .attr("viewBox", [0, 0, plotSize, plotSize]) .style("width", plotSize) .style("height", plotSize); // Construct the unit disk in pixel space (to re-use below). svg.append("circle") .attr("id", "circle") .attr("cx", plotSize / 2) .attr("cy", plotSize / 2) .attr("r", plotSize / 2 - 2); // Create a clipping path from the unit disk. svg.append("clipPath") .attr("id", "circleClip") .append("use") .attr("href", "#circle"); // Initialize the contours drawing. svg.append("g") .attr("id", "contours") .attr("fill", "none") .attr("stroke", "#777") .attr("stroke-opacity", 0.5) .attr("clip-path", "url(#circleClip)"); // Draw the disk outline. svg.append("use") .attr("href", "#circle") .attr("fill", "none") .attr("stroke-width", 2) .attr("stroke", "#000"); // Define a callback to compute the contours for specified Zernike coefficients. function update(nollIndices, coefs, {plotQuantity="Z"}={}) { if(nollIndices.length != coefs.length) throw new Error( `Coefficients array does not have expected length ${nollIndices.length}.`); const grid = Float64Array.from({length:x.length}, (_,k) => { let sum = 0.0; for(let i=0; i < coefs.length; ++i) { const {n,l} = nollDecode(nollIndices[i]); if(coefs[i] != 0) sum += coefs[i] * calculated[plotQuantity](n, l, k); } return sum; }); // See https://observablehq.com/@d3/contours const [lo, hi] = d3.extent(grid); const thresholds = d3.range(numContours).map(i => lo + (hi - lo) / (numContours - 1) * i); const contours = d3.contours() .size([ncell, ncell]).thresholds(thresholds)(grid) .map(transform); const color = d3.scaleSequential([hi, lo], d3.interpolateRdBu); // Draw the computed contours for this grid. svg.select("#contours") .selectAll("path") .data(contours) .join("path") .attr("fill", d => color(d.value)) .attr("d", d3.geoPath()); } update([], []); return Object.assign(svg.node(), { update }); }34:[] 2e:{"title":"Zernike Polynomials","cells":[{"id":0,"value":"# Zernike Polynomials\n\nA complete set of orthogonal polynomials for scalar fields on a disk, useful for analyzing optical systems. For the corresponding vector-field polynomials, see the [Zhao-Burge notebook](https://observablehq.com/@dkirkby/zhao-burge).","pinned":false,"hidden":false,"mode":"md"},{"id":783,"value":"| Version | Changes |\n| ------- | -------- |\n| v858 | Add Cartesian gradients dZdx, dZdy |\n| v780 | Initial release |","pinned":false,"hidden":false,"mode":"md"},{"id":694,"value":"Use the sliders below to mix different Zernike polynomials:","pinned":false,"hidden":false,"mode":"md"},{"id":601,"value":"viewof coefInputs = zernikeCoeffientInputs({nollIndices:[2,3,8,7,4], value:[-0.3,-0.5,0.5,0,0.4]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":848,"value":"Chose whether to plot the Zernike polynomial or one of its Cartesian derivatives,\n${tex`\\partial Z/\\partial x`} and ${tex`\\partial Z/\\partial y`}:","pinned":false,"hidden":false,"mode":"md"},{"id":791,"value":"viewof what = {\n const choices = [\"Z\", \"dZdx\", \"dZdy\"];\n return Inputs.radio(choices, {label: \"Quantity\", value: choices[0], format: d=>d + \"(n,l)\"});\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":469,"value":"plotZ = plotZernike(viewof coefInputs.maxRadialDegree)","pinned":false,"hidden":false,"mode":"ojs"},{"id":731,"value":"{\n plotZ.update(viewof coefInputs.nollIndices, coefInputs, {plotQuantity:what});\n return html`(the sliders and plot are linked by this cell)`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":770,"value":"Compare the plots below with the [corresponding figure](https://en.wikipedia.org/wiki/Zernike_polynomials#/media/File:Zernike_polynomials_with_read-blue_cmap.png) from the [wikipedia page](https://en.wikipedia.org/wiki/Zernike_polynomials):","pinned":false,"hidden":false,"mode":"md"},{"id":312,"value":"viewof n = Inputs.radio(d3.range(1, plotMaxRadialDegree+1), {label: \"n\", value: 3})","pinned":false,"hidden":false,"mode":"ojs"},{"id":316,"value":"viewof l = Inputs.radio(d3.range(-n, n+1, 2), {label: \"l\", value: n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":422,"value":"label = {\n const index = d3.index(namedZernikes, d => d.n, d => d.l);\n const {j} = nollIndex(n, l);\n let label = `n=${n},l=${l},j=${j}`;\n const name = index.get(n)?.get(l)?.name;\n if(name) label += `: ${name}`;\n plotSingle.update([j], [1.]);\n return html`${label}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":747,"value":"plotSingle = plotZernike(5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":397,"value":"### Usage Instructions","pinned":false,"hidden":false,"mode":"md"},{"id":399,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"id":394,"value":"### Implementation and Tests","pinned":false,"hidden":false,"mode":"md"},{"id":384,"value":"Evaluate arbitrary Zernike polynomials at arbitrary points on the unit disk using recurrence relations from [this paper](https://doi.org/10.1117/12.2020389) (but with sin/cos swapped). Returns 1 for any points outside the unit disk.","pinned":false,"hidden":false,"mode":"md"},{"id":6,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":828,"value":"Calculate the decompositions of ${tex`\\partial Z/\\partial x`} and ${tex`\\partial Z/\\partial y`} into Zernikes, for comparison with Tables 7 and 14 from [here](https://wp.optics.arizona.edu/visualopticslab/wp-content/uploads/sites/52/2016/08/Zernike-Notes-15Jan2016.pdf):","pinned":false,"hidden":false,"mode":"md"},{"id":807,"value":"function dZdxCoefs(n, l) {\n const {sqrt,abs} = Math;\n const coefs = [ ];\n if(n == 0) return coefs;\n const norm = (l == 0) ? 2*(n + 1) : n + 1;\n let np, lp, norm2;\n const absl = abs(l), sign = l == 0 ? 1 : l / absl;\n for(np = absl+1; np <= n - 1; np += 2) {\n lp = sign * (absl + 1);\n coefs.push({n,l,np,lp,csq:norm*(np + 1)});\n }\n if((l != 0) && (l != -1)) {\n norm2 = (l == 1) ? 2 : 1;\n for(np = absl-1; np <= n - 1; np += 2) {\n lp = sign * (absl - 1);\n coefs.push({n,l,np,lp,csq:norm*norm2*(np + 1)});\n }\n }\n return coefs;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":824,"value":"function dZdyCoefs(n, l) {\n const {sqrt,abs} = Math;\n const coefs = [ ];\n if(n == 0) return coefs;\n let norm = (l == 0) ? 2*(n + 1) : n + 1;\n let np, lp, norm2;\n const absl = abs(l), sign = l == 0 ? 1 : l / absl;\n norm *= sign;\n for(np = absl+1; np <= n - 1; np += 2) {\n lp = sign * (absl + 1);\n coefs.push({n,l,np,lp:-lp,csq:norm*(np + 1)});\n }\n if((l != 0) && (l != 1)) {\n norm2 = (l == -1) ? 2 : 1;\n for(np = absl-1; np <= n - 1; np += 2) {\n lp = sign * (absl - 1);\n coefs.push({n,l,np,lp:-lp,csq:-norm*norm2*(np + 1)});\n }\n }\n return coefs;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":263,"value":"Tabulate the named Zernike polynomials with radial degree up to 4:","pinned":false,"hidden":false,"mode":"md"},{"id":143,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":265,"value":"Define a utility function to create a rectangular grid covering the unit disk. Returns arrays [x,y] of length ~ ngrid**2 * PI/4.","pinned":false,"hidden":false,"mode":"md"},{"id":230,"value":"function diskXYGrid(ngrid, {clip=true}={}) {\n const xygrid = d3.range(ngrid).map((_,i) => (i - (ngrid-1)/2) / ((ngrid-1)/2));\n let grid = d3.cross(xygrid, xygrid);\n if(clip) grid = grid.filter(([x,y]) => x*x + y*y <= 1);\n return d3.transpose(grid);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":632,"value":"Convert between the Noll index `j` and the physical indices `n` (radial) and `l` (azimuthal):","pinned":false,"hidden":false,"mode":"md"},{"id":634,"value":"nollIndex = (n, l) => {\n const base = n * (n + 1) / 2 + Math.abs(l);\n const t = (n % 4) < 2;\n let j;\n if(t) j = base + ((l <= 0) ? 1 : 0);\n else j = base + ((l >= 0) ? 1 : 0);\n return {n, l, j};\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":647,"value":"nollDecode = (j) => {\n const {ceil, sqrt} = Math;\n const n = ceil(0.5 * (sqrt(8 * j + 1) - 3));\n const t = (n % 4) < 2;\n let l = j - n * (n + 1) / 2;\n if(t) {\n if(j % 2 == 1) l = -(l - 1);\n }\n else {\n if(j % 2 == 0) l -= 1;\n else {\n l = -l;\n if((l == -1) && (n % 2 == 0)) l = 0;\n }\n }\n return {n, l, j};\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":268,"value":"Compare the explicit named polynomials with the generic calculation:","pinned":false,"hidden":false,"mode":"md"},{"id":8,"value":"unitTestPassed = {\n const {abs} = Math;\n const atol = 1e-12;\n const [x, y] = diskXYGrid(3);\n const {Z,r,phi} = calculateZernike(4, x, y, {debug:true});\n namedZernikes.map(({n,l,j,Zj}) => {\n const {j:j2} = nollIndex(n, l);\n if(j != j2) throw new Error(`Invalid Noll index ${j2} for (n,l,j)=(${n},${l},${j})`);\n const {n:n2, l:l2} = nollDecode(j);\n if(n != n2 || l != l2) throw new Error(`Invalid Noll decode (${n},${l}) != (${n2},${l2}) for j=${j}`); \n const close = r.map((r,k) => {\n const Z1 = Zj(r, phi[k]);\n const Z2 = Z(n,l,k);\n const prefix = `(n,l,j,k)=(${n},${l},${j},${k}) (x,y)=(${x[k]},${y[k]}) (r,phi)=(${r},${phi[k]})`;\n if(!Number.isFinite(Z2)) throw new Error(`${prefix} is ${Z2}`);\n if(abs(Z1 - Z2) > atol) throw new Error(`${prefix} fails atol check: ${Z1} != ${Z2}`);\n });\n });\n return true;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":359,"value":"Benchmark different array types. The default is `Float64Array` but you can change this when you call `calculateZernike`.","pinned":false,"hidden":false,"mode":"md"},{"id":361,"value":"benchmark = {\n const [x, y] = diskXYGrid(150);\n const data = [Array, Float32Array, Float64Array].map(arrayType => {\n const {elapsed} = calculateZernike(15, x, y, {profile:true});\n return {\"Array Type\": arrayType.name, \"Elapsed Time (ms)\": elapsed};\n });\n return Inputs.table(data, {\n width: 400,\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":701,"value":"Utilities for setting coefficients and plotting arbitrary mixtures of Zernike polynomials:","pinned":false,"hidden":false,"mode":"md"},{"id":599,"value":"function zernikeCoeffientInputs({lo=-1, hi=+1, step=0.1, nollIndices=[2,3,4], value}={}) {\n value = value ?? nollIndices.map(_ => 0.5 * (lo + hi));\n const idx = d3.index(namedZernikes, d => d.j);\n const maxNollIndex = d3.max(nollIndices);\n const {n:maxRadialDegree} = nollDecode(maxNollIndex);\n const inputs = nollIndices.map((j, i) => {\n const {n,l} = nollDecode(j);\n let label = `(n=${n},l=${l},j=${j})`;\n if(idx.get(j)) label += ` ${idx.get(j).name}`;\n return Inputs.range([lo, hi], {label, value:value[i], step});\n });\n const form = Inputs.form(inputs);\n return Object.assign(form, {\n nollIndices,\n maxRadialDegree,\n });\n return form;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":467,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":453,"value":"plotData = {\n const ngrid = 51;\n const [x, y] = diskXYGrid(ngrid, {clip:false});\n const {Z} = calculateZernike(plotMaxRadialDegree, x, y, {outerValue:NaN});\n return {Z, xy:d3.transpose([x,y])};\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":123,"value":"plotMaxRadialDegree = 5","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}