30:Ta10,snowflake_map = { const height = width * .88; const scale = width / 2 - 20; const ctx = DOM.context2d(width, height); const piemaps = [snowflake3, snowflake2, snowflake1, gon3, snowflake1, snowflake2] const graticule = d3.geoGraticule() .extentMinor([[-180, -85], [180, 85]]) .stepMinor([30, 30])(); const graticule2 = d3.geoGraticule() .extentMajor([[-180, -85 -1e-5], [180, 85 + 1e-5]]) .extentMinor([[-180, -80], [180, 80]]) .stepMajor([10, 5]) .stepMinor([5, 5])(); let index = 0, time = 0; animate(ctx.canvas, invalidation, () => { if (+new Date - time < 3000) return; // leave each step up for 3 seconds time = +new Date; const piemap = piemaps[index]; const projection = d3 .geoProjection((l, p) => piemap(...d3.geoStereographicRaw(l, p))) .scale(scale) .translate([width / 2, height / 2]) .rotate([-75,-13.92, 11.3]) .precision(1e-1) .clipAngle(90-1e-12); const path = d3.geoPath(projection, ctx); let V = piemap.vertices.map( (x, i) => (i % 2) ? (x-piemap.center[1]) * (-scale) + height/2 : (x-piemap.center[0]) * scale + width/2); ctx.fillStyle = "#ffffff"; ctx.fillRect(0, 0, width, height); ctx.fillStyle = JCh_to_hex([0.21, 0.013, -2]); ctx.beginPath(); ctx.moveTo(V[0], V[1]); for (let i = 0, n = V.length; i < n; i+=2) { ctx.lineTo(V[i], V[i+1]); } ctx.closePath(); ctx.fill(); ctx.strokeStyle = JCh_to_hex([0.17, 0.046, -2]); ctx.lineWidth = 0.1 + 6.5e-4 * width; ctx.beginPath(); path(graticule2); ctx.stroke(); ctx.strokeStyle = JCh_to_hex([0.14, 0.082, -2]); ctx.lineWidth = 0.15 + 7e-4 * width; ctx.beginPath(); path(graticule); ctx.stroke(); ctx.fillStyle = JCha_to_rgba([0.053, 0.07, -1.9], .85); ctx.beginPath(); path(land50); ctx.fill(); ctx.strokeStyle = JCha_to_rgba([0.07, 0.06, -2], .8); ctx.miterLimit = 50; ctx.lineWidth = 0.3 + 1.2e-3 * width; ctx.beginPath(); ctx.moveTo(V[0], V[1]); for (let i = 0, n = V.length; i < n; i+=2) { ctx.lineTo(V[i], V[i+1]); } ctx.closePath(); ctx.stroke(); index = (index + 1) % piemaps.length; }); return md`

${ctx.canvas}

Oblique aspect of a conformal map of a hemisphere in a [snowflake][].
Thanks to [Philippe Rivière](../@fil) for help with map projection code.

[snowflake]: https://en.wikipedia.org/wiki/Koch_snowflake * * * ` }31:T4ca,Intro = md` ## Introduction A [Schwarz–Christoffel map][SC] is a [conformal map][] (that is, it locally preserves shapes) of the [upper half of the complex plane][halfplane] or the complex [unit disk][] to an arbitrary polygon. That this is always possible for any shape is the [Riemann mapping theorem][], stated by [Bernhard Riemann][] [in 1851][rthesis]. In the late 1860s, [Elwin Christoffel][] ([1867][]) and [Hermann Schwarz][] ([1869][]) independently came up with a constructive formula for the case of polygons. [SC]: https://en.wikipedia.org/wiki/Schwarz–Christoffel_mapping [conformal map]: https://en.wikipedia.org/wiki/Conformal_map [halfplane]: https://en.wikipedia.org/wiki/Upper_half-plane [unit disk]: https://en.wikipedia.org/wiki/Unit_disk [Riemann mapping theorem]: https://en.wikipedia.org/wiki/Riemann_mapping_theorem [Bernhard Riemann]: https://en.wikipedia.org/wiki/Bernhard_Riemann [rthesis]: http://www.emis.de/classics/Riemann/Grund.pdf [Elwin Christoffel]: https://en.wikipedia.org/wiki/Elwin_Bruno_Christoffel [1867]: https://doi.org/10.1007/BF02419161 [Hermann Schwarz]: https://en.wikipedia.org/wiki/Hermann_Schwarz [1869]: https://doi.org/10.1515/crll.1869.70.105 ` 32:Tef9,scmap_grid_plot = { const g = polygon; const scale = 640/4; const {PI} = Math; const bezplotwindow = [-2, 2, -2, 2]; const spokes = (n) => { const h = Math.sqrt(2) - 1; return Array.from( linspace(0, 2*Math.PI, n, false), (a) => { const x = Math.cos(a), y = Math.sin(a); return circle_interp([0,0], [h*x, h*y], [x, y]); }); } const circles = (n) => Array.from( linspace(0, 1, n, false), (r) => { return t => [r*Math.cos(2*t*PI), r*Math.sin(2*t*PI)]; }).slice(1); const fillcolor = JCh_to_hex([0.215, 0.005, 0.3]), strokecolor2 = JCh_to_hex([0.17, 0.05, 0.3]), strokecolor1 = JCh_to_hex([0.14, 0.08, 0.3]), strokecolor0 = JCh_to_hex([0.10, 0.12, 0.3]); const math = (markup, x, y) => { const rendered = MathJax.tex2svg(markup).querySelector("svg"); return svg` ${rendered}`} let circleplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [[0,0], ...chunk(g.prevertices, 2)] .map(([cx, cy]) => svg``) const transform = `scale(${scale} ${scale}) translate(.9, 1) scale(0.85, -0.85)`; circleplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} ${math(`z_k`, 41, 273)} `; } let polyplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [g.center, ...chunk(g.vertices, 2)] .map(([cx, cy]) => svg``) const outline = polypts_to_svgpath(g.vertices, 6); const transform = `scale(${scale} ${scale}) translate(3, 1.06) scale(0.85, -0.85)`; polyplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} ${math(`w_k`, 313, 292)} `; } return html`

`}33:Tf87,scmap_angles_plot = { const g = polygon; const scale = 640/1.8; const {PI} = Math; const bezplotwindow = [-2, 2, -2, 2]; const spokes = (n) => { const h = Math.sqrt(2) - 1; return Array.from( linspace(0, 2*Math.PI, n, false), (a) => { const x = Math.cos(a), y = Math.sin(a); return circle_interp([0,0], [h*x, h*y], [x, y]); }); } const circles = (n) => Array.from( linspace(0, 1, n, false), (r) => { return t => [r*Math.cos(2*t*PI), r*Math.sin(2*t*PI)]; }).slice(1); const fillcolor = JCh_to_hex([0.217, 0.0025, 0.3]), strokecolor3 = JCh_to_hex([0.19, 0.02, 0.3]), strokecolor2 = JCh_to_hex([0.18, 0.04, 0.3]), strokecolor1 = JCh_to_hex([0.13, 0.09, 0.3]), strokecolor0 = JCh_to_hex([0.10, 0.12, 0.3]), angleecolor = JCh_to_hex([0.10, 0.12, -2.2]); let polyplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [[0,0], ...chunk(g.vertices, 2)] .map(([cx, cy]) => svg``) const outline = polypts_to_svgpath(g.vertices, 6); const transform = `scale(${scale} ${scale}) translate(1.15, 0.05) scale(0.85, -0.85)`; const math = (markup, x, y) => { const rendered = MathJax.tex2svg(markup).querySelector("svg"); return svg` ${rendered}`} const V = g.vertices; polyplot = svg` ${spokespath(24*6*6)} ${circlespath(6*6*6)} ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${dots} ${arc_arrow([V[18], V[19]], [V[20], V[21]], [V[16], V[17]], false, .11, .007)} ${arc_arrow([V[18], V[19]], [2*V[18]-V[16], 2*V[19]-V[17]], [V[20], V[21]], true, .11, .007)} ${arc_arrow([V[20], V[21]], [V[22], V[23]], [V[18], V[19]], false, .085, .007)} ${arc_arrow([V[20], V[21]], [2*V[20]-V[18], 2*V[21]-V[19]], [V[22], V[23]], true, .12, .007)} ${math(`w_{k-1}`, 96, 26)} ${math(`w_k`, 62, 277)} ${math(`w_{k+1}`, 419, 159)} ${math(`w_{k+2}`, 476, 284)} ${math(`\\alpha_k\\pi`, 116, 243)} ${math(`\\beta_k\\pi`, 120, 306)} ${math(`\\alpha_{k+1}\\pi`, 428, 96)} ${math(`\\beta_{k+1}\\pi`, 490, 152)} `; } return html`

`}34:T979, md` Toby Driscoll’s fantastic [SC Toolbox][] ([user’s guide][sct guide], [github][sct github]) is a [Matlab][] library for constructing and evaluating such maps. The trickiest part – which we won’t discuss here – is a non-linear optimization problem to find which points on the real line or unit circle (the “prevertices” ${$`\{z_k\}`}) should map onto the vertices of the polygon. After the prevertices have been found, evaluating the map is a [numerical integration][] problem, and the map can be inverted using [Newton’s method][] or similar. While [conceptually relatively simple][driscoll square], this evaluation can still be computationally expensive if the polygon has a large number of vertices (and finding inverses in particular is slow, because it requires several forward evaluations), and accuracy sometimes suffers from rounding errors very close to the corners. In this notebook I’ll propose a method for evaluating maps from the disk, using prevertices found by SC Toolbox, that is nearly as accurate as its built-in evaluation but two orders of magnitude faster. I am building on the work of [Abi Gopal][] & Nick Trefethen (2019) [“Representation of Conformal Maps by Rational Functions”][conformal rat], who globally approximate a Schwarz–Christoffel map (among other types of conformal maps) using a single rational function. Their method is much faster than numerical integration, but since Schwarz–Christoffel maps in particular involve a power function singularity at each vertex, the degree of the approximation must grow higher than I’d prefer to fit the data accurately near any sharp corners, and accuracy also suffers near sharp corners from rounding errors in the construction of the approximation, based in part on the precision limits of floating point numbers away from the origin. [SC Toolbox]: http://www.math.udel.edu/~driscoll/SC/ [sct guide]: http://www.math.udel.edu/~driscoll/SC/guide.pdf [sct github]: https://github.com/tobydriscoll/sc-toolbox [Matlab]: https://en.wikipedia.org/wiki/MATLAB [Newton’s method]: https://en.wikipedia.org/wiki/Newton%27s_method [numerical integration]: https://en.wikipedia.org/wiki/Numerical_integration [driscoll square]: http://www.chebfun.org/examples/complex/ConformalSquare.html [Abi Gopal]: https://www.maths.ox.ac.uk/people/abinand.gopal/ [conformal rat]: https://arxiv.org/abs/1804.08127 `35:T214b,scpie_steps_and_slices_plot = { const g = polygon; const k = 9; const exponent = g.exponents[k] const scale = 640/4; const {PI, sqrt, atan2, sin, exp, log, atan, cos} = Math; const bezplotwindow = [-4, 4, -2, 4]; // elementwise midpoint between two arrays of points on the unit circle const midpoint = (v0, v1) => { const sum = v0.map(([x, y], i) => [v1[i][0] + x, v1[i][1] + y]); return sum.map(([x, y]) => [x / Math.hypot(x, y), y / Math.hypot(x, y)]); } let wedges0, wedges1, wedgehues, hplmap, complex_power; { const v0 = Array.from(chunk(g.prevertices, 2)); const v1 = rotate(v0.slice(), 1); const piewidths = Array.from(chunk(g.piewidths, 2)); const m = midpoint(v0, v1); wedgehues = m.map(([x, y]) => atan2(y, x)); const m0 = midpoint(v0, m); const m1 = midpoint(m, v1); const h = ([x, y]) => [x * (sqrt(2) - 1), y * (sqrt(2) - 1)]; wedges0 = v0.map((_,i) => [ circle_interp(v0[i], h(v0[i]), [0,0]), circle_interp([0,0], h(m[i]), m[i]), circle_interp(m[i], m0[i], v0[i]) ]); wedges1 = v1.map((_,i) => [ circle_interp(m[i], h(m[i]), [0,0]), circle_interp([0,0], h(v1[i]), v1[i]), circle_interp(v1[i], m1[i], m[i]) ]); hplmap = ([x, y]) => { const d0x = x - v0[k][0], d0y = y - v0[k][1], d1x = x - v1[k][0], d1y = y - v1[k][1], hx = piewidths[k][0], hy = piewidths[k][1]; const // q = denominator of moebius map qx = d0x - (hx * d1x + hy * d1y), qy = d0y - (hx * d1y - hy * d1x); x = (d0x * qx + d0y * qy) * (1 / (qx * qx + qy * qy)); y = (d0y * qx - d0x * qy) * (1 / (qx * qx + qy * qy)); x = Math.max(-10, Math.min(x, 10)); y = Math.max(-10, Math.min(y, 10)); return [x, y]; }; complex_power = ([x, y]) => { const // complex power r = exp(0.5 * log(x*x + y*y) * exponent), a = atan2(y, x) * exponent; return [1.75*r * cos(a), 1.75*r * sin(a)]; } } // The functions for chroma here are a bit arbitrary, but are trying // to work around the unavailability of colorful teal in sRGB. // The more appropriate remedy would be to write some code to bring // JzCzhz colors into sRGB gamut along constant lightness/hue rays. const wedge0colors = wedgehues.map((h) => JCh_to_hex([0.125, 0.062 + .015 * sin(h) + .015 * sin(h) * sin(h), h])) const wedge1colors = wedgehues.map((h) => JCh_to_hex([0.17, 0.06 + .011 * sin(h) + .011 * sin(h) * sin(h), h])) const strokecolor0 = JCh_to_hex([0.10, 0.0, 0.3]); let circleplot; { const wedges0path = wedges0 .map((fs) => fs.map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const transform = `scale(${scale} ${scale}) translate(.9, 1) scale(0.85, -0.85)`; circleplot = svg` ${wedges0path} ${wedges1path} `; } let hplplot; { const clip = DOM.uid("hplclip"); const wedges0path = wedges0 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const transform = `scale(${scale} ${scale}) translate(3.0, 1.8) scale(0.85, -0.85)`; hplplot = svg` ${wedges0path}Z ${wedges1path}Z `; } let wedgeplot; { const bezplotwindow = [-1, 3, -1, 3]; const clip = DOM.uid("wedgeclip"); const line = ([ax, ay], [bx, by]) => { if (Math.hypot(ax, ay) > 10) return []; return [ ax, ay, 2*ax/3 + bx/3, 2*ay/3 + by/3, ax/3 + 2*bx/3, ay/3 + 2*by/3, bx, by]; }; const wedges0path = wedges0 .map((fs) => ([ line(complex_power(hplmap(fs[0](1e-10))), complex_power(hplmap(fs[0](0.01)))), bezeval((t) => complex_power(hplmap(fs[0](t))), [0.01, 1], bezplotwindow, 1e-4), bezeval((t) => complex_power(hplmap(fs[1](t))), [0, 0.99], bezplotwindow, 1e-4), line(complex_power(hplmap(fs[1](.99))), complex_power(hplmap(fs[1](1)))), ]).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => ([ line(complex_power(hplmap(fs[0](0))), complex_power(hplmap(fs[0](0.01)))), bezeval((t) => complex_power(hplmap(fs[0](t))), [0.01, 1], bezplotwindow, 1e-4), bezeval((t) => complex_power(hplmap(fs[1](t))), [0, 0.99], bezplotwindow, 1e-4), line(complex_power(hplmap(fs[1](.99))), complex_power(hplmap(fs[1](1-1e-10)))), ]).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const transform = `scale(${scale} ${scale}) translate(2.25, 3.95) scale(0.85, -0.85)`; wedgeplot = svg` ${wedges0path}Z ${wedges1path}Z `; } let polyplot; { const wedges0path = wedges0 .map((fs) => fs.map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const outline = polypts_to_svgpath(g.vertices, 6); const transform = `scale(${scale} ${scale}) translate(0.95, 3.2) scale(0.85, -0.85)`; polyplot = svg` ${wedges0path} ${wedges1path} `; } const math = (markup, x, y) => { const rendered = MathJax.tex2svg(markup).querySelector("svg"); return svg` ${rendered}`} const r = String.raw; return html`

`}36:T463, md` 1. A [Möbius map][] ${$`M_k^+`} or ${$`M_k^-`} from the unit disk to the upper or lower half-plane, where the nearest prevertex ${$`z_k`} is mapped to ${$`0`}, and the adjacent prevertex is mapped to ${$`1`}. 2. The power function ${$`z \mapsto z^{\alpha_k}`} to cancel the singularity for this corner. 3. A [barycentric][] rational approximation for a map from the wedge of the complex plane to a polygon centered at the vertex and scaled so that ${$`0 \mapsto 0`} and ${$`1 \mapsto 1`}. 4. A [similarity][] map of the form ${$`\left.z \mapsto w_k(1-z) + w_{k+1}z\right.`} to put the output into its final place in the desired output polygon, scaled and rotated appropriately to fit the edge between points ${$`w_k`} and ${$`w_{k+1}`}. Because floating point numbers have high precision near the origin, the rational approximation in step (3) can be fast and accurate, even very close to the corner. [Möbius map]: https://en.wikipedia.org/wiki/Möbius_transformation [barycentric]: https://people.maths.ox.ac.uk/trefethen/barycentric.pdf [similarity]: https://en.wikipedia.org/wiki/Similarity_(geometry) `37:T197a,moebius_plot = { const g = polygon; const k = 9; const scale = 640/4; const {PI, sqrt, atan2, sin, exp, log, atan, cos} = Math; const bezplotwindow = [-20, 20, -4, 4]; // elementwise midpoint between two arrays of points on the unit circle const midpoint = (v0, v1) => { const sum = v0.map(([x, y], i) => [v1[i][0] + x, v1[i][1] + y]); return sum.map(([x, y]) => [x / Math.hypot(x, y), y / Math.hypot(x, y)]); } let wedges0, wedges1, wedgehues, hplmap, hplmap2; { const v0 = Array.from(chunk(g.prevertices, 2)); const v1 = rotate(v0.slice(), 1); const piewidths = Array.from(chunk(g.piewidths, 2)); const m = midpoint(v0, v1); wedgehues = m.map(([x, y]) => atan2(y, x)); const m0 = midpoint(v0, m); const m1 = midpoint(m, v1); const h = ([x, y]) => [x * (sqrt(2) - 1), y * (sqrt(2) - 1)]; wedges0 = v0.map((_,i) => [ circle_interp(v0[i], h(v0[i]), [0,0]), circle_interp([0,0], h(m[i]), m[i]), circle_interp(m[i], m0[i], v0[i]) ]); wedges1 = v1.map((_,i) => [ circle_interp(m[i], h(m[i]), [0,0]), circle_interp([0,0], h(v1[i]), v1[i]), circle_interp(v1[i], m1[i], m[i]) ]); hplmap = ([x, y]) => { const d0x = x - v0[k][0], d0y = y - v0[k][1], d1x = x - v1[k][0], d1y = y - v1[k][1], hx = piewidths[k][0], hy = piewidths[k][1]; const // q = denominator of moebius map qx = d0x - (hx * d1x + hy * d1y), qy = d0y - (hx * d1y - hy * d1x); x = (d0x * qx + d0y * qy) * (1 / (qx * qx + qy * qy)); y = (d0y * qx - d0x * qy) * (1 / (qx * qx + qy * qy)); x = Math.max(-10, Math.min(x, 10)); y = Math.max(-10, Math.min(y, 10)); return [x, y]; }; hplmap2 = ([x, y]) => { [x, y] = hplmap([x, y]); return [1 - x, -y]; } } // The functions for chroma here are a bit arbitrary, but are trying // to work around the unavailability of colorful teal in sRGB. // The more appropriate remedy would be to write some code to bring // JzCzhz colors into sRGB gamut along constant lightness/hue rays. const wedge0colors = wedgehues.map((h) => JCh_to_hex([0.125, 0.062 + .015 * sin(h) + .015 * sin(h) * sin(h), h])) const wedge1colors = wedgehues.map((h) => JCh_to_hex([0.17, 0.06 + .011 * sin(h) + .011 * sin(h) * sin(h), h])) const strokecolor0 = JCh_to_hex([0.10, 0.0, 0.3]); let circleplot; { const wedges0path = wedges0 .map((fs) => fs.map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const transform = `scale(${scale} ${scale}) translate(.9, 1.12) scale(0.85, -0.85)`; circleplot = svg` ${wedges0path} ${wedges1path} `; } let hplplot; { const clip = DOM.uid("hplclip"); const wedges0path = wedges0 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const transform = `scale(${scale} ${scale}) translate(3.2, 1.0) scale(0.425, -0.425)`; hplplot = svg` ${wedges0path}Z ${wedges1path}Z `; } let hplplot2; { const clip = DOM.uid("hplclip"); const wedges0path = wedges0 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap2(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap2(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const transform = `scale(${scale} ${scale}) translate(3.2, 1.2) scale(0.425, -0.425)`; hplplot2 = svg` ${wedges0path}Z ${wedges1path}Z `; } const math = (markup, x, y) => { const rendered = MathJax.tex2svg(markup).querySelector("svg"); return svg` ${rendered}`} const r = String.raw; return html`

`}38:T5a9,md` These maps (which I recommend evaluating directly as written below instead of converting to a standard form) are: ${responsive_katex(invalidation, $$`\begin{aligned} M_k^+(z) = \tilde{z} &= \dfrac{z - z_k}{(z - z_k) - \bar{h}_k(z - z_{k+1})} \\[1.5em] M_{k+1}^-(z) = 1 - \tilde{z} &= \dfrac{z - z_{k+1}}{(z - z_{k+1}) - h_k(z - z_k)} \end{aligned}`)} We can compute ${$`m_k`} and ${$`h_k`} as: ${responsive_katex(invalidation, $$`\begin{aligned} m_k &= \pm \dfrac{z_k + z_{k+1}}{\left| z_k + z_{k+1} \right|} = i\dfrac{z_k - z_{k+1}}{\left| z_k - z_{k+1} \right|} \\[1.5em] h_k &= \pm \dfrac{(z_k + z_{k+1}) \bar{z}_k}{\left| (z_k + z_{k+1}) \bar{z}_k \right|} = i\dfrac{(z_k - z_{k+1})\bar{z}_k}{\left| (z_k - z_{k+1})\bar{z}_k \right|} \end{aligned}`)} To avoid loss of precision and guarantee the correct sign, if the angle between ${$`z_k`} and ${$`z_{k+1}`} is less than 1/4 turn, use the additive form. If the angle is between 1/4 turn and 3/4 turn, use the subtractive form. If the angle is more than 3/4 turn, use the additive form but flip the sign. The inverses of ${$`M_k^+`} and ${$`M_{k+1}^-`} can be computed, respectively: ${responsive_katex(invalidation, $$`\begin{aligned} \tilde{z} &\mapsto z = z_k + \dfrac{(z_{k+1} - z_k)\tilde{z}}{\tilde{z} + h_k(1 - \tilde{z})} \\[1.5em] 1 - \tilde{z} &\mapsto z = z_{k+1} + \dfrac{(z_k - z_{k+1})(1 - \tilde{z})}{(1 - \tilde{z}) + \bar{h}_k\tilde{z}} \end{aligned}`)} ` 39:T83b,SC_formula = md` ## Schwarz–Christoffel formula The typical formula for the Schwarz–Christoffel map from the upper half plane to a polygon looks like (up to similarity): ${$$` \int^{\tilde{z}} \prod_j\left(\zeta - \tilde{z}_j\right)^{-\beta_j} \mathrm{d}\zeta `} Since we are starting from the origin (where our prevertex ${$`z_k`} has been mapped) I propose to use a (scaled and rotated) form which retains slightly more precision near the origin, multiplying the above by ${$`\prod_{j\neq k} (-\tilde{z}_j)^{\beta_j}`}, a complex constant: ${responsive_katex(invalidation, $$`\begin{aligned} & \int^{\tilde{z}} { \zeta^{-\beta_k} \prod_{j \neq k} {\left(1 - \dfrac{\zeta}{\tilde{z}_j}\right)^{\!-\beta_j}\! \mathrm{d}\zeta}} \\[1.5em] & \quad = \int^{\tilde{z}} {\exp \Biggl( -\beta_k \log{\zeta} - \sum_{j \neq k} { \beta_j \operatorname{log1p} \biggl(\dfrac{-\zeta}{\tilde{z}_j}\biggr) } \Biggr) \mathrm{d}\zeta} \end{aligned}`)} Where ${$`\operatorname{log1p}(z) = \log(1 + z)`}; see [Wikipedia][log1p]. One convenient feature of this version is that as ${$`\zeta \to 0`}, the inside of the integral approaches ${$`\zeta^{-\beta_k}`}. Dividing each term in the product by the prevertex is same idea as the common form of the Schwarz–Christoffel formula used for the disk: ${responsive_katex(invalidation, $$`\begin{aligned} & \int^z { \prod_j{\left(1 - \frac{\zeta}{z_j}\right)^{\!-\beta_j}} \!\mathrm{d}\zeta } \\ & \quad = \int^z {\exp \Biggl( - \sum_j { \beta_j \operatorname{log1p} \biggl(\dfrac{-\zeta}{z_j}\biggr) } \Biggr) \mathrm{d}\zeta} \end{aligned}`)} Except that we need a special case to avoid dividing by ${$`\tilde{z}_k = 0`}. If we parametrize a contour in the complex plane by a function ${$`\zeta(t)`} with its domain as the interval ${$`[0, 1]`}, then we can be more explicit: ${responsive_katex(invalidation, $$` \int_0^1 \exp\Biggl( -\beta_k \log{\zeta(t)} - \sum_{j \neq k} {\beta_j \operatorname{log1p}\biggl(\dfrac{-\zeta(t)}{\tilde{z}_j}\biggr)} \Biggr) \zeta'(t)\, \mathrm{d}t `)} [log1p]: https://en.wikipedia.org/wiki/Natural_logarithm#lnp1 ` 3a:Tf56,md` Since the inside of this integral has a singularity of type ${$`z \mapsto z^{-\beta_k}`} at the origin, we could use [Gauss–Jacobi quadrature][] (named after [Carl Gauss][] ([1814][gauss 1814]) and [Carl Jacobi][] ([1826][jacobi 1826])) to evaluate it numerically at a particular point. See Nick Trefethen (1980) [“Numerical Computation of the Schwarz–Christoffel Transformation”][lnt 1980] and [Nick Hale][] & [Alex Townsend][] (2013) [“Fast and Accurate Computation of Gauss–Legendre and Gauss–Jacobi Quadrature Nodes and Weights”][hale townsend]. We can save doing redundant work by separating out the singularity, and approximating the smooth part as a polynomial (relying on [Chebfun’s][chebfun] ability to pick whatever degree polynomial is necessary for an accurate approximation). We can use Nick Hale & al.’s [\`singfun/cumsum\`][cumsum] function to get an indefinite integral over our interval which we can cheaply evaluate at whatever intermediate points we like. The math behind \`singfun/cumsum\` is described in an [unpublished technical note][cumsum ref] by Nick Hale & [Sheehan Olver][], “Numerical Computation of Indefinite Integrals for Functions with Poles or Algebraic Singularities”, and the original prototype came from Mark Richardson’s 2009 thesis, [*Approximating Divergent Functions in the Chebfun System*][richardson 2009]. In [\`singfun\`][singfun], we would represent the function to be integrated as the product of the singularity at the interval’s left endpoint ${$`(2t)^{-\beta_k}`} and a polynomial approximation of: ${responsive_katex(invalidation, $$` \exp\Biggl( -\beta_k \log\biggl(\dfrac{\zeta(t)}{2t}\biggr) - \sum_{j \neq k} {\beta_j \operatorname{log1p}\biggl(\dfrac{-\zeta(t)}{\tilde{z}_j}\biggr)} \Biggr) \zeta'(t) `)} Which for typical Schwarz–Christoffel maps and the contours we care about should be smooth and relatively compact to represent. The ${$`2`} here shows up because \`singfun\` internally uses the interval ${$`[-1, 1]`}, with left and right singularities of the form ${$`(1 + x)^{\lambda}`} and ${$`(1 - x)^{\rho}`}, but we are instead using the interval ${$`[0, 1]`}, which is half as wide. After taking an indefinite integral with \`cumsum\`, which will result in a product of a polynomial approximation and ${$`(2t)^{\alpha_k}`}, we should separately evaluate the polynomial and the power function and multiply the results together, because if we use the built-in \`feval\` in \`singfun\`, it will first transplant our domain to ${$`[-1,1]`}, which will lose precision for values very close to ${$`0`}. In Matlab code: > \`\`\`matlab left = chebfun(@(t) 2*t, [0 1]); df = left .^ (-beta) .* ... chebfun(@(t) exp(... stuff ...).*dzeta(t), [0 1]); fsmooth = left .^ (beta - 1) .* cumsum(df); f = @(t) fsmooth(t) .* (2*t) .^ (1 - beta); \`\`\` [Gauss–Jacobi quadrature]: https://en.wikipedia.org/wiki/Gauss–Jacobi_quadrature [Carl Gauss]: https://en.wikipedia.org/wiki/Carl_Friedrich_Gauss [gauss 1814]: https://gallica.bnf.fr/ark:/12148/bpt6k2412190.langEN [Carl Jacobi]: https://en.wikipedia.org/wiki/Carl_Gustav_Jacob_Jacobi [jacobi 1826]: https://doi.org/10.1515/crll.1826.1.301 [lnt 1980]: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1980_1.pdf [Nick Hale]: http://appliedmaths.sun.ac.za/~nhale/ [Alex Townsend]: http://pi.math.cornell.edu/~ajt/ [hale townsend]: http://www.chebfun.org/publications/HaleTownsend2013a.pdf [chebfun]: http://www.chebfun.org [cumsum]: https://github.com/chebfun/chebfun/blob/master/%40singfun/cumsum.m [cumsum ref]: https://github.com/chebfun/chebfun/blob/master/%40singfun/cumsum.m#L29-L30 [Sheehan Olver]: http://wwwf.imperial.ac.uk/~solver/ [richardson 2009]: http://solo.bodleian.ox.ac.uk/primo-explore/fulldisplay?docid=oxfaleph017190265&vid=SOLO [singfun]: https://github.com/chebfun/chebfun/blob/master/%40singfun/singfun.m `3b:T52a,Contours = md` ## Sector boundary contours To find points along the boundary of our sector of the unit disk (the “pie slice”) for constructing a boundary correspondence to use in building an AAA approximation, we can apply our Möbius map above to the spoke between ${$`z_k`} and the origin, the circular arc from ${$`z_k`} to ${$`m_k`}, and the spoke from ${$`m_k`} to the origina. After transformation, the circular arc becomes the interval ${$`\left[0, \tfrac{1}{2}\right]`}, parametrized by ${$`\zeta(t) = \tfrac{1}{2}t`}. The two transformed spokes can be [written as parametric functions][circle interp] of the interval ${$`\left[0, 1\right]`} (via another Möbius map): [circle interp]: https://observablehq.com/@jrus/circle-arc-interpolation ${responsive_katex(invalidation, $$`\begin{aligned} \zeta_{z_k}(t) &= \dfrac{\tilde{s}_k \tilde{o}t} {(\tilde{o} - \tilde{s}_k)(1-t) + \tilde{s}_kt} \\[2em] \zeta_{m_k}(t) &= \tfrac{1}{2} + i \dfrac{\operatorname{Im}(\tilde{r}_k)\operatorname{Im}(\tilde{o})t} {\operatorname{Im}(\tilde{o} - \tilde{r}_k)(1-t) + \operatorname{Im}(\tilde{r}_k)t} \end{aligned}`)} Where ${$`\tilde{s}_k`} and ${$`\tilde{r}_k`} are the Möbius-mapped midpoints of the spokes, and ${$`\tilde{o}`} is the Möbius-mapped origin. `3c:T1354,contours_plot = { const g = polygon; const k = 9; const exponent = g.exponents[k] const scale = 640/4; const {PI, sqrt, atan2, sin, exp, log, atan, cos} = Math; const bezplotwindow = [-4, 4, -2, 4]; // elementwise midpoint between two arrays of points on the unit circle const midpoint = (v0, v1) => { const sum = v0.map(([x, y], i) => [v1[i][0] + x, v1[i][1] + y]); return sum.map(([x, y]) => [x / Math.hypot(x, y), y / Math.hypot(x, y)]); } let wedges0, wedges1, wedgehues, hplmap, complex_power; { const v0 = Array.from(chunk(g.prevertices, 2)); const v1 = rotate(v0.slice(), 1); const piewidths = Array.from(chunk(g.piewidths, 2)); const m = midpoint(v0, v1); wedgehues = m.map(([x, y]) => atan2(y, x)); const m0 = midpoint(v0, m); const m1 = midpoint(m, v1); const h = ([x, y]) => [x * (sqrt(2) - 1), y * (sqrt(2) - 1)]; wedges0 = v0.map((_,i) => [ circle_interp(v0[i], h(v0[i]), [0,0]), circle_interp([0,0], h(m[i]), m[i]), circle_interp(m[i], m0[i], v0[i]) ]); wedges1 = v1.map((_,i) => [ circle_interp(m[i], h(m[i]), [0,0]), circle_interp([0,0], h(v1[i]), v1[i]), circle_interp(v1[i], m1[i], m[i]) ]); hplmap = ([x, y]) => { const d0x = x - v0[k][0], d0y = y - v0[k][1], d1x = x - v1[k][0], d1y = y - v1[k][1], hx = piewidths[k][0], hy = piewidths[k][1]; const // q = denominator of moebius map qx = d0x - (hx * d1x + hy * d1y), qy = d0y - (hx * d1y - hy * d1x); x = (d0x * qx + d0y * qy) * (1 / (qx * qx + qy * qy)); y = (d0y * qx - d0x * qy) * (1 / (qx * qx + qy * qy)); x = Math.max(-10, Math.min(x, 10)); y = Math.max(-10, Math.min(y, 10)); return [x, y]; }; } const wedge0colors = wedgehues.map((h) => JCh_to_hex([0.175, 0.030 + .008 * sin(h) + .008 * sin(h) * sin(h), h])) const wedge1colors = wedgehues.map((h) => JCh_to_hex([0.195, 0.02 + .005 * sin(h) + .005 * sin(h) * sin(h), h])) const strokecolor0 = JCh_to_hex([0.10, 0.0, 0.3]); const outlinecolor = JCh_to_hex([0.05, 0.09, wedgehues[k]]); let hplplot; { const clip = DOM.uid("hplclip"); const wedges0path = wedges0 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const outline = wedges0 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``)[k]; const wedges1path = wedges1 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const dots = wedges0[9].slice(0,2).map((f) => [hplmap(f(0)), hplmap(f(.5)), hplmap(f(1))]).flat(1) .map(([cx, cy]) => svg``) const transform = `scale(${scale} ${scale}) translate(1.5, 1.7) scale(1, -1)`; hplplot = svg` ${wedges0path}Z ${wedges1path}Z ${outline} ${dots} `; } const math = (markup, x, y, scale=1) => { const rendered = MathJax.tex2svg(markup).querySelector("svg"); return svg` ${rendered}`} const r = String.raw; return html`

`}3d:T992,md` The chosen midpoints could be anywhere along the spoke, e.g. halfway instead of ${$`\sqrt{2} - 1`}, the midpoint when the disk represents a stereographically projected hemisphere. I haven’t investigated whether choosing a different midpoint makes any practical difference. Note that ${$`\tfrac{1}{2}`}, ${$`\tilde{o}`}, and ${$`\tilde{r}_k`} all lie along the line with real part ${$`\tfrac{1}{2}`}, so the differences between them are purely imaginary. If we like we could write ${$`i\operatorname{Im}(\tilde{o}) = \tilde{o} - \tfrac{1}{2}`}, ${$`i\operatorname{Im}(\tilde{r}_k) = \tilde{r}_k - \tfrac{1}{2}`}, and ${$`i\operatorname{Im}(\tilde{o} - \tilde{r}_k) = \tilde{o} - \tilde{r}_k`}. These functions have derivatives: ${responsive_katex(invalidation, $$`\begin{aligned} \zeta_{z_k}'(t) &= \dfrac{\tilde{s}_k \tilde{o}(\tilde{o} - \tilde{s}_k)} {\left((\tilde{o} - \tilde{s}_k)(1-t) + \tilde{s}_kt\right)^2} \\[2em] \zeta_{m_k}'(t) &= i\dfrac{\operatorname{Im}(\tilde{r}_k)\operatorname{Im}(\tilde{o})\operatorname{Im}(\tilde{o} - \tilde{r}_k)} {\left(\operatorname{Im}(\tilde{o} - \tilde{r}_k)(1-t) + \operatorname{Im}(\tilde{r}_k)t\right)^2} \end{aligned} `)} The Schwarz–Christoffel map of the line segment ${$`\zeta_{m_k}`} from ${$`\tfrac{1}{2}`} to ${$`\tilde{o}`} should be representable as a simple polynomial, as there are no singularities nearby. The Schwarz–Christoffel maps of both the interval ${$`\left[0, \tfrac{1}{2}\right]`} and the circular arc ${$`\zeta_{z_{k}}`} from ${$`0`} to ${$`\tilde{o}`} can be represented as the product of a polynomial times a function with ${$`z \mapsto z^{\alpha_k}`} singularity. If we integrate from the origin to the point ${$`\tilde{o}`} along the path consisting of the interval ${$`\left[0, \tfrac{1}{2}\right]`} followed by the arc ${$`\zeta_{m_k}(t)`}, we sometimes get a (slightly) different result than integrating directly along the arc ${$`\zeta_{z_k}(t)`}. Both of these endpoints can also fail to exactly match the “conformal center” according to the data in the SC Toolbox diskmap, due to some combination of rounding errors and inaccuracy in the provided prevertices. This leads approximations from different slices to slightly disagree along the seams and at their common center point. I haven’t found a good way to fix these discrepancies that doesn’t degrade accuracy of the subsequent approximations. `3e:T79a,Scale = { const narrow = width < 550; const r = String.raw; return md` ## Scaling the output To fit our approximation into its final position, we can examine our two integrals for the segments sharing the same output edge, ${responsive_katex(invalidation, $$`\begin{aligned} ${narrow ? `&` : ``} \int^{\tilde{z}} { \zeta^{-\beta_k} \prod_{j \neq k} {\left(1 -\dfrac{\zeta}{\tilde{z}_j}\right)^{\!-\beta_j}\! \mathrm{d}\zeta}} ${narrow ? r`\\[1em] &\qquad\qquad = ` : `& =`} \prod_{j\neq k} \left(-\tilde{z}_j\right)^{\beta_j}\int^{\tilde{z}} \prod_j\left(\zeta - \tilde{z}_j\right)^{-\beta_j} \mathrm{d}\zeta \\[2em] ${narrow ? `&` : ``} \int^{1 - \tilde{z}} { \zeta^{-\beta_{k+1}}\!\! \prod_{j \neq k+1}\! {\left(1 -\dfrac{\zeta}{1 - \tilde{z}_j}\right)^{\!-\beta_j}\! \mathrm{d}\zeta}} ${narrow ? r`\\[1em] &\qquad\qquad = ` : `& =`} -\!\! \prod_{j\neq {k+1}}\! (1 - \tilde{z}_j)^{\beta_j}\int^{\tilde{z}} \prod_j\left(\zeta - \tilde{z}_j\right)^{-\beta_j} \mathrm{d}\zeta \end{aligned}`)} Remember that ${$`M_{k+1}^-(z) = 1 - \tilde{z} = 1 - M_k^+(z)`}. If we want to normalize the scales of these two integrals, we can multiply the second one (or divide the first one) by: ${responsive_katex(invalidation, $$`\begin{aligned} & \prod_{\substack{j \neq k+1\\[.2em] j \neq k}} \left(\dfrac{1 - \tilde{z}_j}{-\tilde{z}_j}\right)^{-\beta_j} = \prod_{\substack{j \neq k+1\\[.2em] j \neq k}} \left(1 - \dfrac{1}{\tilde{z}_j}\right)^{-\beta_j} \\[2em] & \qquad = \exp \Biggl( -\! \sum_{\substack{j \neq k+1\\[.2em] j \neq k}} { \beta_j \operatorname{log1p} \biggl(\dfrac{-1}{\tilde{z}_j}\biggr) } \Biggr) \end{aligned}`)} Since the two adjacent prevertices have been mapped to ${$`\tilde{z}_k = 0`} and ${$`\tilde{z}_{k+1} = 1`}, for every other value of ${$`j`} we must have ${$`\tilde{z}_j < 0`} or ${$`1 < \tilde{z}_j`}, which implies that ${$`1 > 1 / \tilde{z}_j`}, so this product is always a positive real number. `} 3f:T435,md` Notice that: ${$$`\operatorname{log1p}\biggl(\dfrac{-1}{\tilde{z}_j}\biggr) = -\operatorname{log1p}\biggl(\dfrac{-1}{1 - \tilde{z}_j}\biggr)`} So in practice in my code, I use either ${$`M_{k}^+(z_j) = \tilde{z}_j`} or ${$`M_{k+1}^-(z_j) = 1 - \tilde{z}_j`} based on whether ${$`z_j`} is closer to ${$`z_k`} or ${$`z_{k+1}`}, respectively, to avoid some loss of precision in cases where the two pie slices from ${$`k`} to ${$`k+1`} are wider than the nearby slices: ${responsive_katex(invalidation, $$` \exp \Biggl( - \sum_{\substack{j \\[.15em] \tilde{z}_j < 0}} { \beta_j \operatorname{log1p} \biggl(\dfrac{-1}{\tilde{z}_j}\biggr) } + \sum_{\substack{j \\[.15em] \tilde{z}_j > 1}} { \beta_j \operatorname{log1p} \biggl(\dfrac{-1}{1 - \tilde{z}_j}\biggr) } \Biggr) `)} After normalization, the magnitudes of each integral along the interval ${$`\left[0, \tfrac{1}{2}\right]`} in its respective domain should sum to the length of the whole edge of the output polygon. This then tells us how to scale each piece of the approximation to fit into its final shape. `40:T6c0,Approx = md` ## Constructing a rational approximation Once we can approximate the boundaries of our segment by complex functions of the unit interval, we can form an approximation over the whole sector by finding a large number of points on the boundary both in the “input” half-plane and in the “output” polygon, and applying our power function to these “input” values to eliminate the singularity at the origin. We can feed these into Chebfun’s [\`aaa\`][aaa] function, treated as a black box, to obtain forward and inverse rational approximations. My approach to choosing points on the boundary has been to try to roughly space them along each curved edge like the [Chebyshev points][chebpts] *after* the power function has been applied. So for boundary edges adjacent to the prevertex, I take the Chebyshev points for the interval ${$`[0, 1]`}, apply a power function ${$`t \mapsto t^{1/\alpha_k}`}, and then for the AAA “input” values apply ${$`t \mapsto \zeta(t)`} followed by ${$`z \mapsto z^{\alpha_k}`}, whereas for the “output” values I apply the result of \`cumsum\` from above. For the edge that is not adjacent to a prevertex, there’s no need to skew the input points towards one end of the interval. I haven’t experimented much with varying the point distribution or the number of points used here; it’s possible equally spaced points would work just as well, or that substantially more or fewer points would make some material difference. (The illustration below uses the same point distribution, but for clarity reduces the number of points.) [aaa]: https://github.com/chebfun/chebfun/blob/master/aaa.m [chebpts]: https://www.chebfun.org/ATAP/atap-first6chapters.pdf#page=10 ` 41:T22b9,boundary_points_plot = { const g = polygon; const k = 9; const exponent = g.exponents[k] let corner_to_origin; { const vert = Array.from(chunk(g.vertices, 2)); // [sidex, sidey] = inverse of displacement from w_k to w_k+1 let sidex = vert[k+1][0] - vert[k][0], sidey = vert[k+1][1] - vert[k][1]; const sider2 = sidex*sidex + sidey*sidey; sidex /= sider2; sidey /= -sider2; corner_to_origin = ([x, y]) => { x -= vert[k][0]; y -= vert[k][1]; return [x*sidex - y*sidey, x*sidey + y*sidex]; }; } const origin_vertices = Array.from(chunk(g.vertices, 2)) .map(corner_to_origin).flat(1); const scale = 640/4; const {PI, sqrt, atan2, sin, exp, log, atan, cos} = Math; const bezplotwindow = [-4, 4, -2, 4]; // elementwise midpoint between two arrays of points on the unit circle const midpoint = (v0, v1) => { const sum = v0.map(([x, y], i) => [v1[i][0] + x, v1[i][1] + y]); return sum.map(([x, y]) => [x / Math.hypot(x, y), y / Math.hypot(x, y)]); } let wedges0, wedges1, wedgehues, hplmap, complex_power; { const v0 = Array.from(chunk(g.prevertices, 2)); const v1 = rotate(v0.slice(), 1); const piewidths = Array.from(chunk(g.piewidths, 2)); const m = midpoint(v0, v1); wedgehues = m.map(([x, y]) => atan2(y, x)); const m0 = midpoint(v0, m); const m1 = midpoint(m, v1); const h = ([x, y]) => [x * (sqrt(2) - 1), y * (sqrt(2) - 1)]; wedges0 = v0.map((_,i) => [ circle_interp(v0[i], h(v0[i]), [0,0]), circle_interp([0,0], h(m[i]), m[i]), circle_interp(m[i], m0[i], v0[i]) ]); wedges1 = v1.map((_,i) => [ circle_interp(m[i], h(m[i]), [0,0]), circle_interp([0,0], h(v1[i]), v1[i]), circle_interp(v1[i], m1[i], m[i]) ]); hplmap = ([x, y]) => { const d0x = x - v0[k][0], d0y = y - v0[k][1], d1x = x - v1[k][0], d1y = y - v1[k][1], hx = piewidths[k][0], hy = piewidths[k][1]; const // q = denominator of moebius map qx = d0x - (hx * d1x + hy * d1y), qy = d0y - (hx * d1y - hy * d1x); x = (d0x * qx + d0y * qy) * (1 / (qx * qx + qy * qy)); y = (d0y * qx - d0x * qy) * (1 / (qx * qx + qy * qy)); x = Math.max(-10, Math.min(x, 10)); y = Math.max(-10, Math.min(y, 10)); return [x, y]; }; complex_power = ([x, y]) => { const // complex power r = exp(0.5 * log(x*x + y*y) * exponent), a = atan2(y, x) * exponent; return [r * cos(a), r * sin(a)]; } } const wedge0colors = wedgehues.map((h) => JCh_to_hex([0.175, 0.030 + .008 * sin(h) + .008 * sin(h) * sin(h), h])) const wedge1colors = wedgehues.map((h) => JCh_to_hex([0.195, 0.02 + .005 * sin(h) + .005 * sin(h) * sin(h), h])) const strokecolor0 = JCh_to_hex([0.10, 0.0, 0.3]); const dotscolor = JCh_to_hex([0.05, 0.09, wedgehues[k]]); const dots_exp = 1 / g.exponents[k]; const pts = [ chebpts(33).map((x) => 0.5 * (x + 1)).map((x)=> Math.pow(x, dots_exp)), chebpts(33).map((x) => 0.5 * (x + 1)), chebpts(33).map((x) => 0.5 * (x + 1)).map((x)=> 1 - Math.pow(x, dots_exp)) ]; let hplplot; { const clip = DOM.uid("hplclip"); const wedges0path = wedges0 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.slice(0,2).map((f) => bezeval((t) => hplmap(f(t)), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const dots = wedges0[k] .map((f, i) => pts[i].map((t) => hplmap(f(t)))) .flat(1) .map(([cx, cy]) => svg``); const transform = `scale(${scale} ${scale}) translate(0.15, 2.5) scale(1.7, -1.7)`; hplplot = svg` ${wedges0path}Z ${wedges1path}Z ${dots} `; } let wedgeplot; { const bezplotwindow = [-1, 3, -1, 3]; const clip = DOM.uid("wedgeclip"); const line = ([ax, ay], [bx, by]) => { if (Math.hypot(ax, ay) > 10) return []; return [ ax, ay, 2*ax/3 + bx/3, 2*ay/3 + by/3, ax/3 + 2*bx/3, ay/3 + 2*by/3, bx, by]; }; const wedges0path = wedges0 .map((fs) => ([ line(complex_power(hplmap(fs[0](1e-10))), complex_power(hplmap(fs[0](0.01)))), bezeval((t) => complex_power(hplmap(fs[0](t))), [0.01, 1], bezplotwindow, 1e-4), bezeval((t) => complex_power(hplmap(fs[1](t))), [0, 0.99], bezplotwindow, 1e-4), line(complex_power(hplmap(fs[1](.99))), complex_power(hplmap(fs[1](1)))), ]).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => ([ line(complex_power(hplmap(fs[0](0))), complex_power(hplmap(fs[0](0.01)))), bezeval((t) => complex_power(hplmap(fs[0](t))), [0.01, 1], bezplotwindow, 1e-4), bezeval((t) => complex_power(hplmap(fs[1](t))), [0, 0.99], bezplotwindow, 1e-4), line(complex_power(hplmap(fs[1](.99))), complex_power(hplmap(fs[1](1-1e-10)))), ]).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const dots = wedges0[k] .map((f, i) => pts[i].map((t) => complex_power(hplmap(f(t))))) .flat(1) .map(([cx, cy]) => svg``); const transform = `scale(${scale} ${scale}) translate(2.25, 1.15) scale(1.7, -1.7)`; wedgeplot = svg` ${wedges0path}Z ${wedges1path}Z ${dots} `; } let polyplot; { const clip = DOM.uid("wedgeclip"); const wedges0path = wedges0 .map((fs) => fs.map((f) => bezeval((t) => corner_to_origin(g(...f(t))), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const wedges1path = wedges1 .map((fs) => fs.map((f) => bezeval((t) => corner_to_origin(g(...f(t))), [0, 1], bezplotwindow, 1e-4)).flat(1)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d, i) => svg``); const dots = wedges0[k] .map((f, i) => pts[i].map((t) => corner_to_origin(g(...f(t))))) .flat(1) .map(([cx, cy]) => svg``); const outline = polypts_to_svgpath(origin_vertices, 6); // const transform = `scale(${scale} ${scale}) translate(2.95, 2.2) scale(0.85, -0.85)`; const transform = `scale(${scale} ${scale}) translate(2.25, 2.4) scale(1.7, -1.7)`; polyplot = svg` ${wedges0path} ${wedges1path} ${dots} `; } const math = (markup, x, y) => { const rendered = MathJax.tex2svg(markup).querySelector("svg"); return svg` ${rendered}`} const r = String.raw; return html`

`}42:T4a3,LNT_lecture = { const ytwidth = Math.min(width, 640) - 40; const ytheight = ytwidth * 9/16; return md` For details about the approximation method, see [Yuji Nakatsukasa][], [Olivier Sète][], & Nick Trefethen (2018) [“The AAA Algorithm for Rational Approximation”][aaa paper]. Or for context watch Nick Trefethen’s July 2020 lecture, “Rational Functions”:

${youtube('S1upJPMIFfg', ytwidth, ytheight)}

While we’re at it, let’s clip the slide at 56 minutes:

We might summarize this notebook as a switch from using a more general method which “aims for root-exponential convergence” in favor of a more specific method where we “break the domain into pieces and ‘nail the singularities’ for each piece”. [Yuji Nakatsukasa]: https://people.maths.ox.ac.uk/nakatsukasa [Olivier Sète]: https://www.math.tu-berlin.de/fachgebiete_ag_modnumdiff/fg_numerische_lineare_algebra/v_menue/mitarbeiter/dr_olivier_sete/publications [aaa paper]: https://arxiv.org/abs/1612.00337 `} 43:T7e9,Evaluation = md` ## Evaluation in slices Once we have saved the data for a particular diskmap, we can evaluate any list of points in the domain by first finding which of them fit in each slice, and then applying the 4 steps described above for that slice to the relevant points. I cribbed the idea and some of the structure of the Matlab code for this from [\`chebfun/feval\`][feval]. See [Ricardo Pachón][], [Rodrigo Platte][], & Nick Trefethen (2010) [“Piecewise-Smooth Chebfuns”][piecewise] for context. (In Javascript, it instead makes more sense to iterate one point at a time.) To find whether a complex-valued input lies between a pair of adjacent prevertices, I check that it lies to the left of the first prevertex and to the right of the next: namely ${$`\operatorname{Im}(\bar{z}_{k}z) \ge 0`} and ${$`\operatorname{Im}(\bar{z}_{k+1}z) \lt 0`}. It would probably work just as well to compute the [phase][] of each point and directly compare them. To choose between the two adjacent slices closest to ${$`z_k`} and ${$`z_{k+1}`}, I compare the displacements from the input to the two prevertices, ${$`|z - z_k|`} and ${$`|z - z_{k+1}|`}. To evaluate a rational function in barycentric form, given nodes ${tex`\{\xi_\kappa\}`}, function values ${tex`\{\varphi_\kappa\}`}, and barycentric weights ${tex`\{\lambda_\kappa\}`} (using arbitrary Greek letters to avoid confusion with the notation elsewhere in this notebook), we use the formula: ${responsive_katex(invalidation, $$` z \mapsto \left. \displaystyle{\sum_{\kappa}{\frac{\lambda_\kappa\varphi_\kappa}{z - \xi_\kappa}}} \middle/ \displaystyle{\sum_{\kappa}{\frac{\lambda_\kappa}{z - \xi_\kappa}}} \right. `)} [piecewise]: https://people.maths.ox.ac.uk/trefethen/publication/PDF/2010_134.pdf [feval]: https://github.com/chebfun/chebfun/blob/master/%40chebfun/feval.m#L132-L152 [Ricardo Pachón]: https://ricardopachon.com [Rodrigo Platte]: https://math.la.asu.edu/~platte/ [phase]: https://en.wikipedia.org/wiki/Argument_(complex_analysis) `44:Tebc,lshape_plot = { const g = lshape; const scale = 640/4; const {PI} = Math; const bezplotwindow = [-2, 2, -2, 2]; const spokes = (n) => { const h = Math.sqrt(2) - 1; return Array.from( linspace(0, 2*Math.PI, n, false), (a) => { const x = Math.cos(a), y = Math.sin(a); return circle_interp([0,0], [h*x, h*y], [x, y]); }); } const circles = (n) => Array.from( linspace(0, 1, n, false), (r) => { return t => [r*Math.cos(2*t*PI), r*Math.sin(2*t*PI)]; }).slice(1); const fillcolor = JCh_to_hex([0.215, 0.005, -2.7]), strokecolor2 = JCh_to_hex([0.17, 0.05, -2.7]), strokecolor1 = JCh_to_hex([0.14, 0.08, -2.7]), strokecolor0 = JCh_to_hex([0.10, 0.12, -2.7]); const math = (markup, x, y) => { const rendered = MathJax.tex2svg(markup).querySelector("svg"); return svg` ${rendered}`} let circleplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [[0,0], ...chunk(g.prevertices, 2)] .map(([cx, cy]) => svg``) const transform = `scale(${scale} ${scale}) translate(.9, 1) scale(0.85, -0.85)`; circleplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} `; } let polyplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [g.center, ...chunk(g.vertices, 2)] .map(([cx, cy]) => svg``) const outline = polypts_to_svgpath(g.vertices, 6); const transform = `scale(${scale} ${scale}) translate(2.2, 1.85) scale(0.85, -0.85)`; polyplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} `; } return html`

`}45:Td75,md` For example, for an [L-shaped region][] SC Toolbox constructs a map with an “apparent accuracy” of ${$`4.06 \times 10^{-13}`}. The \`scpie\` approximation involves 12 rational-function slices of respective degree [16, 15, 16, 23, 23, 16, 15, 16, 16, 11, 11, 16], with the approximations of degree 23 near the inside corner of the L, and the approximations of degree 11 near the outside corner of the L. > ${md` \`\`\`matlab Lpoly = [2; 2+1i; 1+1i; 1+2i; 2i; 0]; Lshape = center(diskmap(polygon(Lpoly), ... scmapopt('Tolerance', 1e-15)), 0.7+0.7i); [Lpie, Lpie_data] = scpie(Lshape); % zm = 1 million quasirandom points in the unit disk plastic = @(j) 2 .* ( ... mod(0.7548776662466927 .* j, 1) - 0.5 + ... (mod(0.5698402909980532 .* j, 1) - 0.5) .* 1i); zm = plastic(transpose(1:2e6)); zm = zm(abs(zm) <= 1); zm = zm(1:1e6); tic; Lzm = Lshape(zm); toc tic; LPzm = Lpie(zm); toc accuracy(Lshape) rms(Lzm - LPzm) norm(Lzm - LPzm, Inf) \`\`\` `} Evaluating the SCToolbox diskmap \`Lshape\` at 1 million [quasirandom points][] in the disk takes 60 seconds on my computer. Evaluating the SCPie approximation \`Lpie\` at the same 1 million points takes less than 0.4 seconds (in Javascript in my web browser, less than 0.25 seconds). The [RMS error][] between the two outputs is ${$`1.69 \times 10^{-13}`}, the median error is ${$`1.18 \times 10^{-13}`}, and the maximum error is ${$`7.70\times 10^{-12}`}. Among the million points only two of them have error worse than ${$`10^{-12}`}, and both are extremely close to the inside corner of the L; for these points the SCPie approximation should be more accurate than the SCToolbox approximation, if we assume the prevertices are accurate. After that, errors are highest in skinny segments where the conformal center as computed from the prevertices slightly differed from the conformal center provided by SC Toolbox. We can get some idea how we are doing by examining the errors between two different maps constructed by SC Toolbox itself: > ${md` \`\`\`matlab LshapeB = center(diskmap(polygon(Lpoly), ... scmapopt('Tolerance', 1e-14)), 0.7+0.7i); LBzm = LshapeB(zm); accuracy(LshapeB) rms(LBzm - Lzm) norm(Lzm - LPzm, Inf) \`\`\` `} This diskmap where we asked for 14 digits instead of 15 yields a map with “apparent accuracy” ${$`4.81 \times 10^{-13}`}. On our set of 1 million quasirandom points, the RMS error between the two maps is ${$`9.89 \times 10^{-13}`}, the median error is ${$`4.36 \times 10^{-13}`}, and the maximum error is ${$`4.49 \times 10^{-11}`}. This is not exactly a fair comparison, because using different prevertices necessarily results in big discrepancies for points nearby. But it shows that our approximation is doing pretty well, and highlights just how important the prevertex data are: for the square root singularities here, if we try to evaluate the prevertices from one map using the other, we get an error on the order of ${$`10^{-6}`}, even though the prevertices themselves only differ by about ${$`10^{-12}`}: > \`\`\`matlab pre1 = prevertex(Lshape); pre2 = prevertex(Lshape2); max(abs(Lshape(pre1) - Lshape2(pre1))) max(abs(pre1 - pre2)) \`\`\` [L-shaped region]: https://www.mathworks.com/content/dam/mathworks/tag-team/Objects/t/72943_92021v00Cleve_L_Shaped_Membrane_Nov_2003.pdf [quasirandom points]: https://observablehq.com/@jrus/plastic-sequence [RMS error]: https://en.wikipedia.org/wiki/Root-mean-square_deviation ` 46:T8d6,Future = md` ## Next steps The obvious immediate thing to do is figure out how to get inverse evaluation working. The inverse for each pie slice is no problem: just run the steps backwards. The tricky bit is (efficiently and reliably) figuring out which pie slice a given point in the polygon belongs to. Perhaps this should be done with a lower-fidelity global approximation. It would also be nice to figure out how to handle SC Toolbox maps from the disk to an unbounded “polygon” with one or more vertices at infinity. Beyond that, it should be straightforward to extend this method to handle SC maps from the half plane to a polygon. Maps from other domains such as a two-infinite-ended strip or a rectangle are less obvious, but I would expect some piecewise rational approximation to work well for them too. The same general approach should also work for “polygons” bounded by circular arcs (or geodesics on the [sphere][] or [hyperbolic plane][]), though it might take a different idea than the power function used in this notebook to handle a corner bounded by arcs tangent or nearly tangent to each-other. Since SC Toolbox does not support these shapes another tool would need to be used to compute the boundary correspondence. See Louis Howell (1993) [“Numerical Conformal Mapping of Circular Arc Polygons”][howell 1993]. With more effort, I expect the approach to quadrature of mapping the unit circle to the real line followed by approximating the SC derivative as a product of a polynomial and a power function should also be effective in the nonlinear optimization problem of determining relative positions of the prevertices. Among other advantages, those integrals can be computed entirely with real arithmetic. I speculate that we could compute SC maps from the disk at close to full precision (possibly by judiciously using double-double arithmetic in a few places), even when the polygon is elongated or prevertices are extremely crowded. But I don’t have the immediate time/inclination to tackle this research project. [sphere]: https://en.wikipedia.org/wiki/Stereographic_projection [hyperbolic plane]: https://en.wikipedia.org/wiki/Poincaré_disk_model [howell 1993]: https://core.ac.uk/download/pdf/82038727.pdf `47:Te1c,drum_grid_plot = { const g = drum; const scale = 640/4; const {PI} = Math; const bezplotwindow = [-5, 5, -5, 5]; const spokes = (n) => { const h = Math.sqrt(2) - 1; return Array.from( linspace(0, 2*Math.PI, n, false), (a) => { const x = Math.cos(a), y = Math.sin(a); return circle_interp([0,0], [h*x, h*y], [x, y]); }); } const circles = (n) => Array.from( linspace(0, 1, n, false), (r) => { return t => [r*Math.cos(2*t*PI), r*Math.sin(2*t*PI)]; }).slice(1); const fillcolor = JCh_to_hex([0.215, 0.005, -0.6]), strokecolor2 = JCh_to_hex([0.17, 0.05, -0.6]), strokecolor1 = JCh_to_hex([0.14, 0.08, -0.6]), strokecolor0 = JCh_to_hex([0.10, 0.10, -0.6]); let circleplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [[0,0], ...chunk(g.prevertices, 2)] .map(([cx, cy]) => svg``) const transform = `scale(${scale} ${scale}) translate(.8, 1.2) scale(0.75, -0.75)`; circleplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} `; } let polyplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [g.center, ...chunk(g.vertices, 2)] .map(([cx, cy]) => svg``) const outline = polypts_to_svgpath(g.vertices, 6); const transform = `scale(${scale} ${scale}) translate(1.7, 2.3) scale(0.75, -0.75)`; polyplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} `; } return html`

`}48:Te0c,gon7_grid_plot = { const g = gon7; const scale = 640/4; const {PI} = Math; const bezplotwindow = [-2, 2, -2, 2]; const spokes = (n) => { const h = Math.sqrt(2) - 1; return Array.from( linspace(0, 2*Math.PI, n, false), (a) => { const x = Math.cos(a), y = Math.sin(a); return circle_interp([0,0], [h*x, h*y], [x, y]); }); } const circles = (n) => Array.from( linspace(0, 1, n, false), (r) => { return t => [r*Math.cos(2*t*PI), r*Math.sin(2*t*PI)]; }).slice(1); const fillcolor = JCh_to_hex([0.215, 0.005, -2.3]), strokecolor2 = JCh_to_hex([0.17, 0.05, -2.3]), strokecolor1 = JCh_to_hex([0.14, 0.08, -2.3]), strokecolor0 = JCh_to_hex([0.10, 0.12, -2.3]); let circleplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval(f, [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [[0,0], ...chunk(g.prevertices, 2)] .map(([cx, cy]) => svg``) const transform = `scale(${scale} ${scale}) translate(.9, 1) scale(0.85, -0.85)`; circleplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} `; } let polyplot; { const spokespath = (n) => spokes(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const circlespath = (n) => circles(n) .map((f) => bezeval((t) => g(...f(t)), [0, 1], bezplotwindow, 1e-4)) .map((b) => bezpts_to_svgpath(b, 6)) .map((d) => svg``); const dots = [g.center, ...chunk(g.vertices, 2)] .map(([cx, cy]) => svg``) const outline = polypts_to_svgpath(g.vertices, 6); const transform = `scale(${scale} ${scale}) translate(3, 1) scale(0.85, -0.85)`; polyplot = svg` ${spokespath(24*6)} ${circlespath(6*6)} ${spokespath(24)} ${circlespath(6)} ${spokespath(1)} ${dots} `; } return html`

`}49:T683,Matlab = md` * * * ## Matlab code ### Example usage: \`\`\` dmap = @(vert, cen, tol) ... center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen); Lshape = dmap([2; 2+1i; 1+1i; 1+2i; 2i; 0], 0.7+0.7i, 1e-14) [Lpie, Lpie_data] = scpie(Lshape) % see http://www.math.udel.edu/~driscoll/research/drums.html Drum = dmap([1+1i; 2+1i; 2; 3+1i; 3+2i; 1+2i; 1+3i; 2i], 1.5+1.5i, 1e-12) [Dpie, Dpie_data] = scpie(Drum) Rectangle = dmap([2; 2+1i; -2+1i; -2], 0.5i, 1e-13) [Rpie, Rpie_data] = scpie(Rectangle, 1e-12) % zp = 10,000 quasirandom points in the unit disk plastic = @(j) 2 .* ( ... mod(0.7548776662466927 .* j, 1) - 0.5 + ... (mod(0.5698402909980532 .* j, 1) - 0.5) .* 1i); zp = plastic((1:2e4).'); zp = zp(abs(zp) <= 1); zp = zp(1:1e4); max_degree = @(data) max(cellfun(@(a) length(a), data.nodes(:))); fstring = '%.1f times faster; rms error %.2e; max error %.2e; max degree %u'; tic; Lz = Lshape(zp); t1 = toc; tic; LPz = Lpie(zp); t2 = toc; sprintf(fstring, t1/t2, rms(Lz - LPz), norm(Lz - LPz, Inf), max_degree(Lpie_data)) tic; Dz = Drum(zp); t1 = toc; tic; DPz = Dpie(zp); t2 = toc; sprintf(fstring, t1/t2, rms(Dz - DPz), norm(Dz - DPz, Inf), max_degree(Dpie_data)) tic; Rz = Rectangle(zp); t1 = toc; tic; RPz = Rpie(zp); t2 = toc; sprintf(fstring, t1/t2, rms(Rz - RPz), norm(Rz - RPz, Inf), max_degree(Rpie_data)) % This plot is nothing new; it uses SC Toolbox. % It is just here to illustrate the distribution of points. figure; hold on; daspect([1 1 1]); axis([-.2 3.2 -.2 3.2]); Dpts = [Dpie_data.prevertices; Dpie_data.prevertices .* Dpie_data.piewidths]; plot(Drum, sqrt((1:15)./16), angle(Dpts)); plot(DPz, 'r.'); \`\`\` ` 4a:T777,md` The above prints out e.g. > \`\`\` 130.2 times faster; rms error 1.69e-13; max error 6.31e-13; max degree 22 127.3 times faster; rms error 1.32e-12; max error 2.96e-11; max degree 25 294.3 times faster; rms error 2.78e-11; max error 3.28e-11; max degree 17 \`\`\` Compare with the global rational approximations from Gopal & Trefethen (2009) for SC Toolbox maps of the same polygons: > \`\`\` 42.3 times faster; rms error 1.4e-08; max error 2.6e-03; degree 58 32.3 times faster; rms error 1.7e-08; max error 4.6e-02; degree 67 87.9 times faster; rms error 1.4e-08; max error 3.0e-02; degree 31 \`\`\` This comparison is just to give a rough idea; it isn’t a fair uniform benchmark. Gopal & Trefethen evaluated a different, smaller set of points; their comparative speed might be a bit faster with the same point set. They were running on a different computer and probably a different version of Matlab. Also, they set SC Toolbox to use lower precision internally since they were only aiming for 8 digits. On the other hand, if we used this piecewise pie-slice method aiming for accuracy of only a few digits, we could get by with a lower degree approximation for each slice, and make evaluation twice as fast. Here’s a plot showing the piece boundaries on the drum shape and the set of quasirandom points we used above (the [plastic sequence][], see [Martin Roberts][] (2018), [“The Unreasonable Effectiveness of Quasirandom Sequences”][roberts 2018]). The radial contours show segment boundaries, and the angular contours are square-root spaced to reach better into the corners:

[plastic sequence]: https://observablehq.com/@jrus/plastic-sequence [Martin Roberts]: http://extremelearning.com.au/about/ [roberts 2018]: http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/ * * * ` 4b:Tbc7,md` Something a bit slower, the second iteration of the [Koch snowflake][]. Cf. Lehel Banjai & Nick Trefethen (2003) [“A Multipole Method for Schwarz–Christoffel Mapping of Polygons with Thousands of Sides”][multipole sc]. On my desktop computer constructing the SC Toolbox diskmap takes 3.5 seconds and constructing the 96-piece rational approximation takes an additional 10 seconds: [Koch snowflake]: https://en.wikipedia.org/wiki/Koch_snowflake [multipole sc]: https://ora.ox.ac.uk/objects/uuid:e0fbd9e3-fd76-4447-b839-0df31a5309e4/download_file?file_format=pdf&safe_filename=NA-02-14.pdf > \`\`\` side = @(p) circshift(p, -1, 1) - p; sf_iter = @(p) reshape(([p,p,p,p] + bsxfun(@times, side(p), ... [0, 1/3, 0.5 - 0.5i/sqrt(3), 2/3])).', [], 1); triangle = [1; -0.5 + 0.5i*sqrt(3); -0.5 - 0.5i*sqrt(3)]; snowflake2 = dmap(circshift(sf_iter(sf_iter(triangle)), -1, 1), 0, 1e-11) [sfpie2, sfpie2_data] = scpie(snowflake2) tic; SF2z = snowflake2(zp); t1 = toc; tic; SFP2w = sfpie2(zp); t2 = toc; sprintf(fstring, t1/t2, rms(SF2z - SFP2w), norm(SF2z - SFP2w, Inf), max_degree(sfpie2_data)) \`\`\` Which prints out e.g.: > \`\`\` 88.7 times faster; rms error 1.35e-11; max error 1.87e-10; max degree 23 \`\`\` If we add one more iteration and drop the requested accuracy to 8 digits, constructing the diskmap takes 48 seconds and throws up an ominous “Warning: Severe crowding”, and then constructing the 288-piece rational approximation takes an additional 43 seconds. Even in this case, my computer can evaluate the map 1 million times in about 2 seconds. (My Javascript implementation used for rendering the hemisphere map at the top of this notebook ends up about 4 times faster than my Matlab implementation, evaluating about 2 million points per second when running on 1 CPU core of my computer.) > \`\`\` snowflake3 = dmap(circshift(sf_iter(sf_iter(sf_iter(triangle))), -1, 1), 0, 1e-8) [sfpie3, sfpie3_data] = scpie(snowflake3) tic; SF3z = snowflake3(zp); t1 = toc; tic; SFP3z = sfpie3(zp); t2 = toc; sprintf(fstring, t1/t2, rms(SF3z - SFP3z), norm(SF3z - SFP3z, Inf), max_degree(sfpie3_data)) figure; hold on; daspect([1 1 1]); axis([-1.15 1.15 -1 1]); S3pts = [sfpie3_data.prevertices; sfpie3_data.prevertices .* sfpie3_data.piewidths]; plot(snowflake3, sqrt((1:15)./16), angle(S3pts)); plot(SFP3z, 'r.'); \`\`\` Which prints out e.g.: > \`\`\` 80.5 times faster; rms error 5.96e-10; max error 7.66e-09; max degree 23 \`\`\` So we are still getting the requested 8 digits here accurately (but this is as about as accurate as SC Toolbox can handle for this polygon). You can see how the individual pie slices start getting pretty skinny, and at this point the overhead of testing every point against every pie slice starts to get a bit expensive; we could probably further improve evaluation speed of this type of map if we put some work in optimizing it.

* * * `4c:T267d,Matlab_scpie = matlab_file('scpie.m', String.raw` function [scpie_diskmap_eval, scpie_diskmap_data] = ... scpie(scdiskmap, aaa_tolerance, npts_per_side) if nargin < 2 aaa_tolerance = accuracy(scdiskmap); end if nargin < 3 npts_per_side = 2 ^ 9; end warning('off', 'CHEBFUN:aaa:Froissart'); % rational quadratic circle-arc interpolant, a Moebius % map f such that f(0) = a, f(1/2) = m, f(1) = b rerp = @(a, m, b) ... (@(t) ((b-m).*a.*(1-t) + (m-a).*b.*t) ... ./ ((b-m).*(1-t) + (m-a).*t)); % derivative of the above function rerpdiff = @(a, m, b) ... (@(t) ((b-m).*(m-a).*(b-a)) ... ./ ((b-m).*(1-t) + (m-a).*t) .^ 2); % same as t -> rerp(0, m, b)(t) / t, but well defined for t = 0 rerp0_quotient = @(m, b) ... (@(t) (m.*b) ... ./ ((b-m).*(1-t) + m.*t)); unit_spoke = rerp(0, sqrt(2) - 1, 1); % a Moebius map from the disk to the upper half plane assuming % a and b are points on the disk and assuming h is a unit complex % number representing half the angle between them, which sends % a -> 0, b -> 1, and 0 to a point on the line Re(z) = 1/2 disk2hpl = @(a, b, h) ... (@(z) (z - a) ./ ((z - a) - conj(h) .* (z - b))); % z -> -1 / disk2hpl(a, b, h)(z) disk2hpl_flip = @(a, b, h) ... (@(z) (conj(h) .* (z - b) - (z - a)) ./ (z - a)); % a chebfun we can raise to a power to represent a power-function % singularity to turn a smooth chebfun into a singfun leftcheb = chebfun(@(x) 2 .* x, [0 1]); left = @(exponent) leftcheb .^ exponent; % Chebyshev points on the interval [0, 1] with endpoints chopped off interval = chebpts(npts_per_side+1, [0 1]); interval = interval(2:end-1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% pre0 = prevertex(scdiskmap); pre1 = circshift(pre0, -1, 1); n = length(pre0); vert0 = vertex(polygon(scdiskmap)); vert1 = circshift(vert0, -1, 1); side = vert1 - vert0; turn_angles0 = 1 - angle(polygon(scdiskmap)); turn_angles1 = circshift(turn_angles0, -1, 1); % 'mid' = midpoint between prevertices % 'piewidths' = mid ./ pre mid = (pre0 + pre1) .* sign(imag(pre1 .* conj(pre0))); far_apart = (mid .* conj(mid)) < 2; mid(far_apart) = 1i .* (pre0(far_apart) - pre1(far_apart)); piewidths = mid .* conj(pre0); mid = mid ./ abs(mid); % normalize piewidths = piewidths ./ abs(piewidths); data = struct; data.center = center(scdiskmap); data.vertices = vert0; data.prevertices = pre0; data.exponents = 1 - turn_angles0; data.piewidths = piewidths; data.nodes = cell(n, 2); data.values = cell(n, 2); data.weights = cell(n, 2); data.errors = cell(n, 2); for k=1:n turn0 = turn_angles0(k); turn1 = turn_angles1(k); pow0 = 1 - turn0; pow1 = 1 - turn1; pow0inv = 1/pow0; pow1inv = 1/pow1; % Moebius maps of the disk to the upper/lower half plane % which put the vertices for this interval at points 0, 1. moe0 = disk2hpl(pre0(k), pre1(k), piewidths(k)); moe1 = disk2hpl(pre1(k), pre0(k), conj(piewidths(k))); moe0flip = disk2hpl_flip(pre0(k), pre1(k), piewidths(k)); moe1flip = disk2hpl_flip(pre1(k), pre0(k), conj(piewidths(k))); % Moebius mapped prevertices, inverted and negated hpre0flip = real(moe0flip(pre0)); hpre1flip = real(moe1flip(pre1)); % sort prevertices roughly from furthest to closest [~, idx0] = sort(abs(turn_angles0 .* hpre0flip)); [~, idx1] = sort(abs(turn_angles1 .* hpre1flip)); idx0 = idx0(idx0 ~= k); idx1 = idx1(idx1 ~= k); hturn0 = turn_angles0(idx0).'; hturn1 = turn_angles1(idx1).'; hpre0flip = hpre0flip(idx0).'; hpre1flip = hpre1flip(idx1).'; % relative scale between adjacent sectors rs = zeros(size(hpre0flip)); n0 = hpre0flip > 0; n1 = hpre1flip > 0; rs(n0) = rs(n0) - hturn0(n0) .* log1p(hpre0flip(n0)); rs(n1) = rs(n1) + hturn1(n1) .* log1p(hpre1flip(n1)); relscale = exp(sum(rs)); % construct endpoints and midpoints to use in defining % parametric functions for boundaries. sk0 = moe0((sqrt(2)-1) * pre0(k)); sk1 = moe1((sqrt(2)-1) * pre1(k)); rk0 = 0.5 + 1i * imag(moe0((sqrt(2)-1) * mid(k))); rk1 = 0.5 + 1i * imag(moe1((sqrt(2)-1) * mid(k))); o0 = 0.5 + 1i * imag(moe0(0)); o1 = 0.5 + 1i * imag(moe1(0)); % parametric functions of the interval [0, 1] onto the % projected pie-slice boundaries in the half plane spoke0 = rerp(0, sk0, o0); spoke1 = rerp(0, sk1, o1); dspoke0 = rerpdiff(0, sk0, o0); dspoke1 = rerpdiff(0, sk1, o1); spoke0q = rerp0_quotient(sk0, o0); % spoke0 / t spoke1q = rerp0_quotient(sk1, o1); % spoke1 / t spoke0m = rerp(0.5, rk0, o0); spoke1m = rerp(0.5, rk1, o1); dspoke0m = rerpdiff(0.5, rk0, o0); dspoke1m = rerpdiff(0.5, rk1, o1); arc = rerp(0, 0.25, 0.5); darc = rerpdiff(0, 0.25, 0.5); arcq = rerp0_quotient(0.25, 0.5); % inside of schwarz-christoffel integrals d_sc_spoke0_smooth = @(t) dspoke0(t) .* ... exp(-turn0 * log(0.5 * spoke0q(t)) - sum(bsxfun(@times, ... hturn0, log1p(bsxfun(@times, spoke0(t), hpre0flip))), 2)); d_sc_spoke1_smooth = @(t) dspoke1(t) .* ... exp(-turn1 * log(0.5 * spoke1q(t)) - sum(bsxfun(@times, ... hturn1, log1p(bsxfun(@times, spoke1(t), hpre1flip))), 2)); d_sc_arc0_smooth = @(t) darc(t) .* ... exp(-turn0 * log(0.5 * arcq(t)) - sum(bsxfun(@times, ... hturn0, log1p(bsxfun(@times, arc(t), hpre0flip))), 2)); d_sc_arc1_smooth = @(t) darc(t) .* ... exp(-turn1 * log(0.5 * arcq(t)) - sum(bsxfun(@times, ... hturn1, log1p(bsxfun(@times, arc(t), hpre1flip))), 2)); d_sc_spoke0m = @(t) dspoke0m(t) .* ... exp(-turn0 * log(spoke0m(t)) - sum(bsxfun(@times, ... hturn0, log1p(bsxfun(@times, spoke0m(t), hpre0flip))), 2)); d_sc_spoke1m = @(t) dspoke1m(t) .* ... exp(-turn1 * log(spoke1m(t)) - sum(bsxfun(@times, ... hturn1, log1p(bsxfun(@times, spoke1m(t), hpre1flip))), 2)); % indefinite integrals sc_spoke0_smooth = left(-pow0) .* ... cumsum(left(-turn0) .* chebfun(d_sc_spoke0_smooth, [0 1])); sc_spoke1_smooth = left(-pow1) .* ... cumsum(left(-turn1) .* chebfun(d_sc_spoke1_smooth, [0 1])); sc_spoke0 = @(t) (2*t) .^ pow0 .* sc_spoke0_smooth(t); sc_spoke1 = @(t) (2*t) .^ pow1 .* sc_spoke1_smooth(t); sc_arc0_smooth = left(-pow0) .* ... cumsum(left(-turn0) .* chebfun(d_sc_arc0_smooth, [0 1])); sc_arc1_smooth = left(-pow1) .* ... cumsum(left(-turn1) .* chebfun(d_sc_arc1_smooth, [0 1])); sc_arc0 = @(t) (2*t) .^ pow0 .* sc_arc0_smooth(t); sc_arc1 = @(t) (2*t) .^ pow1 .* sc_arc1_smooth(t); sc_spoke0m_cheb = cumsum(chebfun(d_sc_spoke0m, [0 1])); sc_spoke1m_cheb = cumsum(chebfun(d_sc_spoke1m, [0 1])); sc_spoke0m = @(t) sc_arc0(1) + sc_spoke0m_cheb(t); sc_spoke1m = @(t) sc_arc1(1) + sc_spoke1m_cheb(t); scale0 = sc_arc0(1) + sc_arc1(1) .* relscale; scale1 = sc_arc1(1) + sc_arc0(1) ./ relscale; inpts0 = [0; 0.5; o0; ... arc(interval .^ pow0inv); ... spoke0m(interval); ... spoke0(interval .^ pow0inv)] .^ pow0; inpts1 = [0; 0.5; o1; ... arc(interval .^ pow1inv); ... spoke1m(interval); ... spoke1(interval .^ pow1inv)] .^ pow1; outpts0 = [0; sc_arc0(1); sc_spoke0(1); ... sc_arc0(interval .^ pow0inv); ... sc_spoke0m(interval); ... sc_spoke0(interval .^ pow0inv)] ./ scale0; outpts1 = [0; sc_arc1(1); sc_spoke1(1); ... sc_arc1(interval .^ pow1inv); ... sc_spoke1m(interval); ... sc_spoke1(interval .^ pow1inv)] ./ scale1; [r0, ~, ~, ~, rnodes0, rvalues0, rweights0, err0] = ... aaa(outpts0, inpts0, 'tol', aaa_tolerance); [r1, ~, ~, ~, rnodes1, rvalues1, rweights1, err1] = ... aaa(outpts1, inpts1, 'tol', aaa_tolerance); data.errors{k, 1} = err0 * abs(side(k)); data.errors{k, 2} = err1 * abs(side(k)); % TODO: figure out inverse evaluation % [s0, ~, ~, ~, snodes0, svalues0, sweights0] = ... % aaa(inpts0, outpts0, 'tol', aaa_tolerance); % [s1, ~, ~, ~, snodes1, svalues1, sweights1] = ... % aaa(inpts1, outpts1, 'tol', aaa_tolerance); data.nodes{k, 1} = rnodes0; data.values{k, 1} = rvalues0; data.weights{k, 1} = rweights0; data.nodes{k, 2} = rnodes1; data.values{k, 2} = rvalues1; data.weights{k, 2} = rweights1; end scpie_diskmap_data = data; scpie_diskmap_eval = @(z) scpie_eval(scpie_diskmap_data, z); end function [out] = scpie_eval(scpie_diskmap_data, z) % equivalent to imag(conj(a) .* b) wedge = @(a, b) real(a) .* imag(b) - imag(a) .* real(b); data = scpie_diskmap_data; % initialize values to 'center' out = ones(size(z)) .* data.center; vertices0 = data.vertices; vertices1 = circshift(vertices0, -1, 1); n = length(vertices0); prevertices0 = data.prevertices; prevertices1 = circshift(prevertices0, -1, 1); exponents0 = data.exponents; exponents1 = circshift(exponents0, -1, 1); % Loop over the sectors corresponding to each polygon side. for k = 1:n z0 = prevertices0(k); z1 = prevertices1(k); if real(data.piewidths(k)) > 0 I = (wedge(z0, z) >= 0) & (wedge(z1, z) < 0); else % rare edge case where the sector is more than half the disk I = (wedge(z0, z) >= 0) | (wedge(z1, z) < 0); end if any(I(:)) zI = z(I); outI = zeros(size(zI)); d0 = zI - z0; d1 = zI - z1; h = data.piewidths(k); % first half sector I0 = abs(d0) <= abs(d1); zI0 = d0(I0) ./ (d0(I0) - conj(h) .* d1(I0)); % moebius map zI0 = zI0 .^ exponents0(k); zI0 = reval(zI0, data.nodes{k, 1}, data.values{k, 1}, data.weights{k, 1}); zI0 = vertices0(k) .* (1 - zI0) + vertices1(k) .* zI0; outI(I0) = zI0; % second half sector I1 = ~I0; zI1 = d1(I1) ./ (d1(I1) - h .* d0(I1)); % moebius map zI1 = zI1 .^ exponents1(k); zI1 = reval(zI1, data.nodes{k, 2}, data.values{k, 2}, data.weights{k, 2}); zI1 = vertices1(k) .* (1 - zI1) + vertices0(k) .* zI1; outI(I1) = zI1; out(I) = outI; end end % infinite or NaN inputs can't be mapped out(isnan(z) | isinf(z)) = NaN; end `)4d:T790,Matlab_scpie_json = matlab_file('scpie_json.m', String.raw` function [] = scpie_json(scpie_diskmap_data, filename) fileID = fopen(filename,'w'); printz1 = @(z) sprintf(' % .16e, % .16e,\n', [real(z), imag(z)].'); printz2 = @(z) sprintf(' % .16e, % .16e,\n', [real(z), imag(z)].'); data = scpie_diskmap_data; output = cell(0); output{end+1} = sprintf('{\n'); txt = printz1(data.center); txt = [sprintf('"center": [\n'), txt(1:end-2), sprintf('\n],\n')]; output{end+1} = txt; txt = printz1(data.vertices); txt = [sprintf('"vertices": [\n'), txt(1:end-2), sprintf('\n],\n')]; output{end+1} = txt; txt = printz1(data.prevertices); txt = [sprintf('"prevertices": [\n'), txt(1:end-2), sprintf('\n],\n')]; output{end+1} = txt; txt = sprintf(' % .16e,\n', data.exponents); txt = [sprintf('"exponents": [\n'), txt(1:end-2), sprintf('\n],\n')]; output{end+1} = txt; txt = printz1(data.piewidths); txt = [sprintf('"piewidths": [\n'), txt(1:end-2), sprintf('\n],\n')]; output{end+1} = txt; txt = cellfun(printz2, data.nodes, 'UniformOutput', false).'; txt = cellfun(@(t) t(1:end-2), txt, 'UniformOutput', false); txt = strjoin(txt(:).', '\n ], [\n'); txt = [sprintf('"nodes": [\n [\n'), txt, sprintf('\n ]\n],\n')]; output{end+1} = txt; txt = cellfun(printz2, data.values, 'UniformOutput', false).'; txt = cellfun(@(t) t(1:end-2), txt, 'UniformOutput', false); txt = strjoin(txt(:).', '\n ], [\n'); txt = [sprintf('"values": [\n [\n'), txt, sprintf('\n ]\n],\n')]; output{end+1} = txt; txt = cellfun(printz2, data.weights, 'UniformOutput', false).'; txt = cellfun(@(t) t(1:end-2), txt, 'UniformOutput', false); txt = strjoin(txt(:).', '\n ], [\n'); txt = [sprintf('"weights": [\n [\n'), txt, sprintf('\n ]\n]\n')]; output{end+1} = txt; output{end+1} = sprintf('}'); cellfun(@(s) fprintf(fileID, '%s', s), output); fclose(fileID); end `)4e:Te8b,scpie_eval = function (scpie_diskmap_data) { const {atan2, cos, exp, log, sin} = Math, rotate1 = (array) => [...array.slice(1), array[0]], // returns [A[1], A[2], ... A[n-1], A[0]] even = (array) => array.filter((x, i) => !(i % 2)), odd = (array) => array.filter((x, i) => i % 2), evenodd = (array) => [even(array), odd(array)], displacement = (array0, array1) => array0.map((x, i) => array1[i] - x), data = scpie_diskmap_data, center = data.center, n = data.vertices.length >> 1, [nodes0, nodes1] = evenodd(data.nodes), [values0, values1] = evenodd(data.values), [weights0, weights1] = evenodd(data.weights), [vertices0x, vertices0y] = evenodd(data.vertices), vertices1x = rotate1(vertices0x), vertices1y = rotate1(vertices0y), sidesx = displacement(vertices0x, vertices1x), sidesy = displacement(vertices0y, vertices1y), [prevertices0x, prevertices0y] = evenodd(data.prevertices), prevertices1x = rotate1(prevertices0x), prevertices1y = rotate1(prevertices0y), [piewidthsx, piewidthsy] = evenodd(data.piewidths), exponents0 = data.exponents, exponents1 = rotate1(exponents0); const scpie_eval = function scpie_eval (x, y) { if (!isFinite(x + y)) return NaN; // figure out which sector the point is in let k; for (k = 0;; k++) { if (k === n) return center.slice(); const z_left_z0 = (prevertices0x[k] * y - prevertices0y[k] * x) >= 0; const z_right_z1 = (prevertices1x[k] * y - prevertices1y[k] * x) < 0; if (z_left_z0 & z_right_z1) break; // special case when the sector between k & k+1 is wider than half the disk: if (piewidthsx[k] <= 0) { if (z_left_z0 | z_right_z1) break; } } const d0x = x - prevertices0x[k], d0y = y - prevertices0y[k], d1x = x - prevertices1x[k], d1y = y - prevertices1y[k], hx = piewidthsx[k], hy = piewidthsy[k]; // half-sector closest to prevertex k if ((d0x*d0x + d0y*d0y) <= (d1x*d1x + d1y*d1y)) { const // q = denominator of moebius map qx = d0x - (hx * d1x + hy * d1y), qy = d0y - (hx * d1y - hy * d1x); x = (d0x * qx + d0y * qy) * (1 / (qx * qx + qy * qy)); y = (d0y * qx - d0x * qy) * (1 / (qx * qx + qy * qy)); const // complex power r = exp(0.5 * log(x*x + y*y) * exponents0[k]), a = atan2(y, x) * exponents0[k]; x = r * cos(a); y = r * sin(a); // rational approximation let out = reval(nodes0[k], values0[k], weights0[k], x, y); x = out[0]; y = out[1]; // transform to final position; reuse out array out[0] = vertices0x[k] + (sidesx[k] * x - sidesy[k] * y); out[1] = vertices0y[k] + (sidesy[k] * x + sidesx[k] * y); return out; } // half-sector closest to prevertex k+1 else { const // q = denominator of moebius map qx = d1x - (hx * d0x - hy * d0y), qy = d1y - (hx * d0y + hy * d0x); x = (d1x * qx + d1y * qy) * (1 / (qx * qx + qy * qy)); y = (d1y * qx - d1x * qy) * (1 / (qx * qx + qy * qy)); const // complex power r = exp(0.5 * log(x*x + y*y) * exponents1[k]), a = atan2(y, x) * exponents1[k]; x = r * cos(a); y = r * sin(a); // rational approximation let out = reval(nodes1[k], values1[k], weights1[k], x, y); x = out[0]; y = out[1]; // transform to final position; reuse out array out[0] = vertices1x[k] - (sidesx[k] * x - sidesy[k] * y); out[1] = vertices1y[k] - (sidesy[k] * x + sidesx[k] * y); return out; } } Object.assign(scpie_eval, scpie_diskmap_data); return scpie_eval; }4f:T436,reval = function reval(nodes, values, weights, x, y) { const n = nodes.length; // running total for numerator & denominator let px = 0, py = 0, qx = 0, qy = 0; for (let j = 0; j < n; j += 2) { // q += w_j / (z - z_j) // p += f_j * w_j / (z - z_j) const xj = x - nodes[j], yj = nodes[j+1] - y, wxj = weights[j], wyj = weights[j+1], fxj = values[j], fyj = values[j+1], dj = 1 / (xj*xj + yj*yj), qxj = dj * (wxj*xj - wyj*yj), qyj = dj * (wxj*yj + wyj*xj); px += fxj * qxj - fyj * qyj; py += fxj * qyj + fyj * qxj; qx += qxj, qy += qyj; } let d = 1 / (qx * qx + qy * qy), fx = (px * qx + py * qy) * d, fy = (py * qx - px * qy) * d; // Edge case where x + iy is one of the nodes; directly use value if (fx + fy !== fx + fy) { // true if x or y is NaN for (let j = 0; j < n; j += 2) { const xj = x - nodes[j], yj = nodes[j+1] - y; if ((xj === 0) & (yj === 0)) { fx = values[j]; fy = values[j+1]; break; } } } return [fx, fy]; }; 50:T691,arc_arrow = (center, endpoint0, endpoint1, internal, radius, linewidth) => { var r0x = endpoint0[0] - center[0], r0y = endpoint0[1] - center[1], r0q = radius / Math.hypot(r0x, r0y), r1x = endpoint1[0] - center[0], r1y = endpoint1[1] - center[1], r1q = radius / Math.hypot(r1x, r1y), clockwise = Math.sign(r0x * r1y - r0y * r1x), clockwise_arrow = +!((clockwise < 0) & internal), arrow_sign = 1 - 2*clockwise_arrow; r0x *= r0q, r0y *= r0q, r1x *= r1q, r1y *= r1q; var startx = center[0] + r0x, starty = center[1] + r0y, endx = center[0] + r1x, endy = center[1] + r1y, // rotate r1 a bit backward so it doesn’t collide with the arrowhead rotc = Math.cos(2*linewidth/radius), rots = arrow_sign * Math.sin(2*linewidth/radius), arcendx = center[0] + r1x * rotc - r1y * rots, arcendy = center[1] + r1y * rotc + r1x * rots; var arrowback0x = endx - 3*linewidth/radius * r1x - arrow_sign * 4 * linewidth/radius * r1y, arrowback0y = endy - 3*linewidth/radius * r1y + arrow_sign * 4 * linewidth/radius * r1x, arrowback1x = endx + 3*linewidth/radius * r1x - arrow_sign * 4 * linewidth/radius * r1y, arrowback1y = endy + 3*linewidth/radius * r1y + arrow_sign * 4 * linewidth/radius * r1x; var arcpath = `M${startx} ${starty} A${radius} ${radius} 0 ${+((clockwise < 0) & !internal)} ${clockwise_arrow} ${arcendx} ${arcendy}` var arrowheadpath = `M${endx} ${endy} L${arrowback0x} ${arrowback0y} L${arrowback1x} ${arrowback1y}Z`; return svg` ` }51:[{"id":"1bc3a9a37deb92ef23a182cb82e0320da15a46de0a300b4fee792f038635841d189fc0973f181133a89de9ead975f8e6d620b87cffedf913c98470d68990422c","href":"https://static.observableusercontent.com/files/1bc3a9a37deb92ef23a182cb82e0320da15a46de0a300b4fee792f038635841d189fc0973f181133a89de9ead975f8e6d620b87cffedf913c98470d68990422c","name":"lshape@1.json","type":"application/json","lastModified":1597735061165.72,"size":34385},{"id":"e037ae26488b58643fb2fd7b907e10b298a666f653aee07c522bb792261e09e47daf2b170102b53d08f8e33cb5dade6994c15046591f61fe2bbdcb5e4de35a7f","href":"https://static.observableusercontent.com/files/e037ae26488b58643fb2fd7b907e10b298a666f653aee07c522bb792261e09e47daf2b170102b53d08f8e33cb5dade6994c15046591f61fe2bbdcb5e4de35a7f","name":"rectangle@1.json","type":"application/json","lastModified":1597735060994.195,"size":20987},{"id":"c53731f0b540c8278abd13b8a7f4c08fc8e01d0ff739ef67102813adc2081c18dbed663a1f4b8e2efd3d0b21c731e15659303d32386b497fbcd1d14ee2212311","href":"https://static.observableusercontent.com/files/c53731f0b540c8278abd13b8a7f4c08fc8e01d0ff739ef67102813adc2081c18dbed663a1f4b8e2efd3d0b21c731e15659303d32386b497fbcd1d14ee2212311","name":"drum@1.json","type":"application/json","lastModified":1597735059729.615,"size":49241},{"id":"b5c19db66c9711dbc5d94b34f78face6cd7fe03cb191b52831b27148fd837ff8d5877bff75dd100d8bcdb336260e8eb65080c9229f48e6b77c0a9937afe6ba0d","href":"https://static.observableusercontent.com/files/b5c19db66c9711dbc5d94b34f78face6cd7fe03cb191b52831b27148fd837ff8d5877bff75dd100d8bcdb336260e8eb65080c9229f48e6b77c0a9937afe6ba0d","name":"gon5.json","type":"application/json","lastModified":1597977558983.673,"size":35300},{"id":"9ec77bda4551301403dbf29042d5b1d717af470c90812be0f55b80afc9545a61fae5e7d5f837b6cdd5d67e4efa4ef7cac92ca2300819808fbc56300aa655e541","href":"https://static.observableusercontent.com/files/9ec77bda4551301403dbf29042d5b1d717af470c90812be0f55b80afc9545a61fae5e7d5f837b6cdd5d67e4efa4ef7cac92ca2300819808fbc56300aa655e541","name":"gon6.json","type":"application/json","lastModified":1597977559374.283,"size":42647},{"id":"bd112728a388a420f379df54f18925b63ce2f3786efbd7e409d69d1242569460388bb02b011823ad87629d45a883cb57036d4acca3716ff073ecda14d547b8ea","href":"https://static.observableusercontent.com/files/bd112728a388a420f379df54f18925b63ce2f3786efbd7e409d69d1242569460388bb02b011823ad87629d45a883cb57036d4acca3716ff073ecda14d547b8ea","name":"gon4.json","type":"application/json","lastModified":1597982966175.828,"size":18719},{"id":"65c3abf453b4268e2e4ace49af6ac2128827d7fc131a24b1e6f5def5f4df7e47bd21e66dd47f08368d6ca536cc60de8cc996c2748906134a04f7f45c93a1551f","href":"https://static.observableusercontent.com/files/65c3abf453b4268e2e4ace49af6ac2128827d7fc131a24b1e6f5def5f4df7e47bd21e66dd47f08368d6ca536cc60de8cc996c2748906134a04f7f45c93a1551f","name":"gon7.json","type":"application/json","lastModified":1597981640173.239,"size":47078},{"id":"55b84a317dfa3f70d0d9bd1bfd1e4debdcf8e70b2db81cf9c5c7088bc1ed5ed4df5d0d80e583056f2c464a8911a8d73438bbdef52f77065874847714f5c12a50","href":"https://static.observableusercontent.com/files/55b84a317dfa3f70d0d9bd1bfd1e4debdcf8e70b2db81cf9c5c7088bc1ed5ed4df5d0d80e583056f2c464a8911a8d73438bbdef52f77065874847714f5c12a50","name":"gon3.json","type":"application/json","lastModified":1597982999034.558,"size":13478},{"id":"f9aa6eac6e445c38ffa36de35ba2a5643eb3537e21413dee4e61854b9331f1761ac877bed248e0855a20c2f57e37e0f33a53778b750be79330c5eb8e825773ea","href":"https://static.observableusercontent.com/files/f9aa6eac6e445c38ffa36de35ba2a5643eb3537e21413dee4e61854b9331f1761ac877bed248e0855a20c2f57e37e0f33a53778b750be79330c5eb8e825773ea","name":"gon8@1.json","type":"application/json","lastModified":1597983227073.769,"size":51185},{"id":"3d1d45c33ffc35c7c1b890376677b2832106f719c56182029bdb30f720f521cb33a23eff5882430768f3668d09bdee462a9846b38fdfeaef657b2140369085cc","href":"https://static.observableusercontent.com/files/3d1d45c33ffc35c7c1b890376677b2832106f719c56182029bdb30f720f521cb33a23eff5882430768f3668d09bdee462a9846b38fdfeaef657b2140369085cc","name":"snowflake1.json","type":"application/json","lastModified":1598386085427.481,"size":80087},{"id":"a0b448303d615734d3229bdb9ef6d9a628d6c49f95b3b4fad0272957fad55b11114fc25bfce19c2d95e4f6f4af3ef97210733a94391126fd87ffa7c07e275a00","href":"https://static.observableusercontent.com/files/a0b448303d615734d3229bdb9ef6d9a628d6c49f95b3b4fad0272957fad55b11114fc25bfce19c2d95e4f6f4af3ef97210733a94391126fd87ffa7c07e275a00","name":"snowflake2@2.json","type":"application/json","lastModified":1598383312036.193,"size":337937},{"id":"0892373c5732d9a1b0893fc560ae743544fa215444652cecdb6c5c0bd4ebff37d89320d4fc2b3cc5167249a137f835cf80bc89c7911085c45fd72e11dbad2144","href":"https://static.observableusercontent.com/files/0892373c5732d9a1b0893fc560ae743544fa215444652cecdb6c5c0bd4ebff37d89320d4fc2b3cc5167249a137f835cf80bc89c7911085c45fd72e11dbad2144","name":"snowflake3@2.json","type":"application/json","lastModified":1598383319021.622,"size":1110785},{"id":"76f251233a9d9140982b972b4d10c08bb714fee296f6952956bcd9b5ad9a50baa4cbd989c6fb600aa1cdb45cb04d48248a11aefa9aede0bd1f9f806dbdf228b2","href":"https://static.observableusercontent.com/files/76f251233a9d9140982b972b4d10c08bb714fee296f6952956bcd9b5ad9a50baa4cbd989c6fb600aa1cdb45cb04d48248a11aefa9aede0bd1f9f806dbdf228b2","name":"polygon@2.json","type":"application/json","lastModified":1598567459448.367,"size":76361},{"id":"277ae543885c902bb64fde025f44a96e86df8c9de6d474ddcff033188c3fe8aa98f366e293faf81c7ee4897e7ea3a513d0100fed1354229e968de56a5588a2ef","href":"https://static.observableusercontent.com/files/277ae543885c902bb64fde025f44a96e86df8c9de6d474ddcff033188c3fe8aa98f366e293faf81c7ee4897e7ea3a513d0100fed1354229e968de56a5588a2ef","name":"lshape convergence.png","type":"image/png","lastModified":1601071666910.549,"size":36133},{"id":"a12337aa274670c3596379ef7c3fd31eb7a7cb941c1e835a64e9c2e61be934e465cdcdae44b675193a3c4a0bb1247150f3dcbff5a84ec0fbd490c800d993b5b5","href":"https://static.observableusercontent.com/files/a12337aa274670c3596379ef7c3fd31eb7a7cb941c1e835a64e9c2e61be934e465cdcdae44b675193a3c4a0bb1247150f3dcbff5a84ec0fbd490c800d993b5b5","name":"nail the singularities@2.jpg","type":"image/jpeg","lastModified":1598828183075.634,"size":114433},{"id":"9d9b24a66aa9ff2e624cc677be61dc96d97c818328705d4e659d2be810a6b9068b18146287e893ff26e9103f8cd6ed64288701fc6f1327f3d9bb836c83d9f734","href":"https://static.observableusercontent.com/files/9d9b24a66aa9ff2e624cc677be61dc96d97c818328705d4e659d2be810a6b9068b18146287e893ff26e9103f8cd6ed64288701fc6f1327f3d9bb836c83d9f734","name":"gon10.json","type":"application/json","lastModified":1598921245262.999,"size":64745},{"id":"0a8453d1468059bd8937b5808bbfe796dddc934f98d567224ccb63fdcd917c8b08a1d2a56822c26fa29e52f33e9772e41a29c32382e4b72b7c87df9df888e039","href":"https://static.observableusercontent.com/files/0a8453d1468059bd8937b5808bbfe796dddc934f98d567224ccb63fdcd917c8b08a1d2a56822c26fa29e52f33e9772e41a29c32382e4b72b7c87df9df888e039","name":"gon9.json","type":"application/json","lastModified":1598921248373.241,"size":60152},{"id":"60bce9f3b894f2270fc04d64216585abee527d710fff93a3e192fd622571e3793909b0b4c264045cfd68584f35b61b88b974c043175a3e0f87c43f359bd8304b","href":"https://static.observableusercontent.com/files/60bce9f3b894f2270fc04d64216585abee527d710fff93a3e192fd622571e3793909b0b4c264045cfd68584f35b61b88b974c043175a3e0f87c43f359bd8304b","name":"drum@2.png","type":"image/png","lastModified":1597218612673.299,"size":234627},{"id":"92f4187516ff79605e6664fa359cc5ee9041c0b7906894fe79d4afc7783391926a2282d64428008ba62a18f0a603e1ac7a90bd3cee055a546eb02cf27dc27888","href":"https://static.observableusercontent.com/files/92f4187516ff79605e6664fa359cc5ee9041c0b7906894fe79d4afc7783391926a2282d64428008ba62a18f0a603e1ac7a90bd3cee055a546eb02cf27dc27888","name":"snowflake.png","type":"image/png","lastModified":1597629414246.957,"size":651910}] 2e:{"title":"Evaluating Schwarz–Christoffel Maps in Pie Slices","cells":[{"id":0,"value":"Title = md`# Evaluating Schwarz–Christoffel Maps in Pie Slices`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1296,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":2056,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1801,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":2404,"value":"md`\nThe gist is that we have some polygon with vertices ${$`\\{w_k\\}`}, for each of which we can find a “prevertex” on the real line or unit circle ${$`z_k`}.\nAt the vertex ${$`w_k`}, we can call the [turning angle][] ${$`\\beta_k = \\tfrac{1}{\\pi}\\operatorname{arg}{\\left((w_{k+1} - w_k) / (w_k - w_{k+1})\\right)}`}, with the internal angle ${$`\\alpha_k = 1 - \\beta_k`}.\nFrom this data we can form the map as an integral involving the prevertices and angles, with arbitrary complex constants ${$`a`} and ${$`c`} chosen to fit the result into the polygon:\n\n${$$`\nw = a + c\\int^z { \\prod_k{\\left(\\zeta - z_k\\right)^{\\!-\\beta_k}} \\mathrm{d}\\zeta }\n`}\n\n [turning angle]: https://en.wikipedia.org/wiki/Internal_and_external_angles\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2393,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":3239,"value":"md`\nFor further background on Schwarz–Christoffel mapping, see [Toby Driscoll][] & [Nick Trefethen][] (2002) [*Schwarz–Christoffel Mapping*][SC Mapping].\nOr for a 19th century source in English, try [Andrew Forsyth][] (1893) *Theory of Functions of a Complex Variable* [§§ 267–275][forsyth 1893].\n\n [Toby Driscoll]: http://tobydriscoll.net\n [Nick Trefethen]: http://people.maths.ox.ac.uk/~trefethen/\n [SC Mapping]: https://www.cambridge.org/core/books/schwarzchristoffel-mapping/E8439E14F3CFE0F8CC2943217A5CA4C4\n [Andrew Forsyth]: https://en.wikipedia.org/wiki/Andrew_Forsyth\n [forsyth 1893]: https://archive.org/details/theoryoffunction028777mbp/page/n699/mode/2up\n\n* * *\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":2429,"value":"Pie_slices = md`\n## Pie slice approximation\n\nTo reduce the required degree of approximation, while improving accuracy, I propose to split the disk up into “pie slice” sectors, a pair for each edge of the polygon.\nThen we can separately approximate each sector near a prevertex ${$`z_k`} as the composition of 4 steps:\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4446,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":25,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":4688,"value":"Moebius = md`\n## Möbius map from disk to half-plane\n\nThe Möbius map ${$`M_k^+`} I propose to use for each pie slice after a prevertex should map ${$`z_k \\mapsto 0`}, ${$`m_k \\mapsto \\tfrac{1}{2}`}, ${$`z_{k+1} \\mapsto 1`}, and the unit disk to the upper half-plane, where ${$`m_k`} is the midpoint on the circle between prevertices ${$`z_k`} and ${$`z_{k+1}`}, and ${$`h_k = m_k/z_k`} is a unit complex number representing half the angle between ${$`z_k`} and ${$`z_{k+1}`}, i.e. the width of our current pie slice.\nThe map ${$`M_k^-`} for each pie slice *before* a prevertex is similar, but sends the unit disk to the *lower* half-plane.\n\nAs a notational convenience, picking out a specific index ${$`k`}, we will call ${$`\\tilde{z} = M_k^+(z)`} and ${$`1 - \\tilde{z} = M_{k+1}^-(z)`} the results of the maps corresponding to the polygon side between ${$`w_k`} and ${$`w_{k+1}`}, so ${$`\\tilde{z}_k = 0`} and ${$`\\tilde{z}_{k+1} = 1`}.\n\n`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":4692,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":4690,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":5276,"value":"md`\nTo compute the inside of the Schwarz-Christoffel integral (next section) at each point ${$`\\zeta`}, we take the quotients ${$`- \\zeta / \\tilde{z}_j`} for each Möbius-mapped prevertex ${$`\\tilde{z}_j`}.\nWe can save a tiny bit of work by directly precomputing these negative inverses, as we don’t ever actually need ${$`\\tilde{z}_j`} per se:\n\n${responsive_katex(invalidation, $$`\n\\dfrac{-1}{\\tilde{z}_j} = \\dfrac{\\bar{h}_k(z_j - z_{k+1}) - (z_j - z_k)}{z_j - z_k}\n`)}\n`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":192,"value":"md`\nIf the prevertices were computed in higher precision (e.g. using [double-double arithmetic][]), the differences in this Möbius map could also be computed at higher precision relatively cheaply, and the overall precision of the Schwarz–Christoffel map could potentially be improved (and the “crowding” effect which sometimes plagues SC maps might be ameliorated – cf. [Lehel Banjai][] (2008) [“Revisiting the Crowding Phenomenon in Schwarz–Christoffel Mapping”][crowding]), even if the rational approximation in the subsequent step were computed at ordinary double precision.\n\n [double-double arithmetic]: https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic\n [Lehel Banjai]: http://www.macs.hw.ac.uk/~lb138/\n [crowding]: http://user.math.uzh.ch/banjai/crowding.pdf\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4667,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":51,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":4643,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":4838,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"},{"id":275,"value":"$3d","pinned":false,"hidden":false,"mode":"ojs"},{"id":410,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs"},{"id":4648,"value":"$3f","pinned":false,"hidden":false,"mode":"ojs"},{"id":422,"value":"$40","pinned":false,"hidden":false,"mode":"ojs"},{"id":5115,"value":"$41","pinned":false,"hidden":false,"mode":"ojs"},{"id":4127,"value":"$42","pinned":false,"hidden":false,"mode":"ojs"},{"id":3360,"value":"$43","pinned":false,"hidden":false,"mode":"ojs"},{"id":4649,"value":"Results = md`\n## Results\n\nIn practice, it seems like we typically get accuracy comparable to the “apparent accuracy” claimed by SC Toolbox for the same maps (or at least within a digit or two), and under the assumption that the prevertices are accurate, accuracy of this method should be *better* than SC Toolbox very close to the prevertices.\nWithout a higher-precision reference known to be accurate, it’s hard to be completely sure about the error.\nEvaluation is typically about 100 times faster than SC Toolbox, though it varies a bit from map to map.\n\nFor the diskmaps I have tried, rational approximations of each sector with local singularities neutralized tend to have degree in the 10–30 range, substantially lower (and faster converging) than the degree required to make a single rational approximation for the whole disk.\n`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":5043,"value":"$44","pinned":false,"hidden":false,"mode":"ojs"},{"id":5365,"value":"lshape_convergence_plot = md`\n

\n`\n// This chart plots the log of the maximum residual error at each step of\n// the AAA algorithm while constructing the rational approximations for the\n// neutralized-singularity pie slices of scpie (in purple) vs. a single\n// rational approximation of the whole disk to the L shape (in teal)","pinned":false,"hidden":false,"mode":"ojs"},{"id":3316,"value":"$45","pinned":false,"hidden":false,"mode":"ojs"},{"id":4652,"value":"md`\nThe pie-slice method does inherently involve (slight) discontinuities along the seams, unlike a global rational approximation.\nIf this presents a problem for some specific application, adjacent approximations could perhaps be blended (and one additional approximation could be added for the central region, e.g. for all input points in the disk with absolute value less than ${$`\\tfrac{1}{2}`}).\n\nAlso, figuring out the tolerance to request from \\`aaa\\` to get the most accurate result takes a bit of trial and error.\nUsing the “apparent accuracy” reported by SC Toolbox seems to be a passable default choice, but in some cases (for instance, regular polygons) it results in approximations that are both less accurate and much higher degree than if we explicitly set the tolerance slightly looser.\n\nKeep scrolling down to see some Matlab code and more concrete examples.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1403,"value":"$46","pinned":false,"hidden":false,"mode":"ojs"},{"id":1313,"value":"Examples = md`\n* * *\n\n## Example maps\n\nHere are a few specific SC maps to play with:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1518,"value":"// An arbitrary looking polygon\npolygon = scpie_eval(await FileAttachment(\"polygon@2.json\").json())\n\n/*\n% Matlab code to generate polygon.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\npol = [ ...\n 0.82 - 1.01i; 1.11 + 0.38i; ...\n 0.49 + 0.46i; 0.80 + 1.03i; ...\n -0.31 + 1.12i; -0.65 + 0.59i; ...\n -0.36 + 0.10i; -1.07 + 0.34i; ...\n -0.88 - 0.04i; -1.05 - 0.88i; ...\n 0.15 - 0.40i; 0.34 - 0.86i; ...\n];\npolmap = dmap(pol, 0, 1e-13);\n[polpie, polpie_data] = scpie(polmap, 1e-11);\nscpie_json(polpie_data, 'polygon.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1222,"value":"// L shaped region, as seen in the Matlab logo\nlshape = scpie_eval(await FileAttachment(\"lshape@1.json\").json())\n\n/*\n% Matlab code to generate lshape.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nLshape = dmap([2; 2+1i; 1+1i; 1+2i; 2i; 0], 0.7+0.7i, 1e-14);\n[Lpie, Lpie_data] = scpie(Lshape, 1e-13);\nscpie_json(Lpie_data, 'lshape.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1185,"value":"// One of a pair of “isospectral drum” shapes which have the same eigenfrequencies.\n// See http://tobydriscoll.net/project/drums/\n// https://en.wikipedia.org/wiki/Hearing_the_shape_of_a_drum\ndrum = scpie_eval(await FileAttachment(\"drum@1.json\").json())\n\n/*\n% Matlab code to generate drum.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nDrum = dmap([1+1i; 2+1i; 2; 3+1i; 3+2i; 1+2i; 1+3i; 2i], 1.5+1.5i, 1e-12);\n[Dpie, Dpie_data] = scpie(Drum, 1e-11);\nscpie_json(Dpie_data, 'drum.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":5323,"value":"$47","pinned":false,"hidden":false,"mode":"ojs"},{"id":1220,"value":"// Rectangle with a 4:1 aspect ratio\nrectangle = scpie_eval(await FileAttachment(\"rectangle@1.json\").json())\n\n/*\n% Matlab code to generate rectangle.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nRectangle = dmap([2; 2+1i; -2+1i; -2], 0.5i, 1e-13);\n[Rpie, Rpie_data] = scpie(Rectangle, 1e-11);\nscpie_json(Rpie_data, 'rectangle.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1673,"value":"// Iteration 1 of the Koch snowflake\nsnowflake1 = scpie_eval(await FileAttachment(\"snowflake1.json\").json())\n\n/*\n% Matlab code to generate snowflake1.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nside = @(p) circshift(p, -1, 1) - p;\nsf_iter = @(p) reshape(([p,p,p,p] + bsxfun(@times, side(p), ...\n [0, 1/3, 0.5 - 0.5i/sqrt(3), 2/3])).', [], 1);\ntriangle = [1; -0.5 + 0.5i*sqrt(3); -0.5 - 0.5i*sqrt(3)];\nsnowflake1 = dmap(circshift(sf_iter(triangle), -1, 1), 0, 1e-14);\n[sfpie1, sfpie1_data] = scpie(snowflake1, 1e-13);\nscpie_json(sfpie1_data, 'snowflake1.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1218,"value":"// Iteration 2 of the Koch snowflake\nsnowflake2 = scpie_eval(await FileAttachment(\"snowflake2@2.json\").json())\n\n/*\n% Matlab code to generate snowflake2.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nside = @(p) circshift(p, -1, 1) - p;\nsf_iter = @(p) reshape(([p,p,p,p] + bsxfun(@times, side(p), ...\n [0, 1/3, 0.5 - 0.5i/sqrt(3), 2/3])).', [], 1);\ntriangle = [1; -0.5 + 0.5i*sqrt(3); -0.5 - 0.5i*sqrt(3)];\nsnowflake2 = dmap(circshift(sf_iter(sf_iter(triangle)), -1, 1), 0, 1e-11);\n[sfpie2, sfpie2_data] = scpie(snowflake2, 1e-11);\nscpie_json(sfpie2_data, 'snowflake2.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1214,"value":"// Iteration 3 of the Koch snowflake\nsnowflake3 = scpie_eval(await FileAttachment(\"snowflake3@2.json\").json())\n\n/*\n% Matlab code to generate snowflake3.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nside = @(p) circshift(p, -1, 1) - p;\nsf_iter = @(p) reshape(([p,p,p,p] + bsxfun(@times, side(p), ...\n [0, 1/3, 0.5 - 0.5i/sqrt(3), 2/3])).', [], 1);\ntriangle = [1; -0.5 + 0.5i*sqrt(3); -0.5 - 0.5i*sqrt(3)];\nsnowflake3 = dmap(circshift(sf_iter(sf_iter(sf_iter(triangle))), -1, 1), 0, 1e-8);\n[sfpie3, sfpie3_data] = scpie(snowflake3, 1e-8);\nscpie_json(sfpie3_data, 'snowflake3.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1594,"value":"md`Regular polygons – equilateral triangle, square, etc.:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1583,"value":"// Equilateral triangle\ngon3 = scpie_eval(await FileAttachment(\"gon3.json\").json())\n\n/*\n% Matlab code to generate gon3.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon3 = dmap(exp(2*pi*i*[1:3].' / 3), 0, 1e-15);\n[G3pie, G3pie_data] = scpie(Gon3, 1e-14)\nscpie_json(G3pie_data, 'gon3.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1580,"value":"// Square\ngon4 = scpie_eval(await FileAttachment(\"gon4.json\").json())\n\n/*\n% Matlab code to generate gon4.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon4 = dmap(exp(2*pi*i*[1:4].' / 4), 0, 1e-15);\n[G4pie, G4pie_data] = scpie(Gon4, 1e-14)\nscpie_json(G4pie_data, 'gon4.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1548,"value":"// Regular pentagon\ngon5 = scpie_eval(await FileAttachment(\"gon5.json\").json())\n\n/*\n% Matlab code to generate gon5.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon5 = dmap(exp(2*pi*i*[1:5].' / 5), 0, 1e-15);\n[G5pie, G5pie_data] = scpie(Gon5, 1e-14)\nscpie_json(G5pie_data, 'gon5.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1551,"value":"// Regular hexagon\ngon6 = scpie_eval(await FileAttachment(\"gon6.json\").json())\n\n/*\n% Matlab code to generate gon6.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon6 = dmap(exp(2*pi*i*[1:6].' / 6), 0, 1e-15);\n[G6pie, G6pie_data] = scpie(Gon6, 1e-14)\nscpie_json(G6pie_data, 'gon6.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":1569,"value":"// Regular heptagon\ngon7 = scpie_eval(await FileAttachment(\"gon7.json\").json())\n\n/*\n% Matlab code to generate gon7.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon7 = dmap(exp(2*pi*i*[1:7].' / 7), 0, 1e-14);\n[G7pie, G7pie_data] = scpie(Gon7, 1e-13)\nscpie_json(G7pie_data, 'gon7.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":5312,"value":"$48","pinned":false,"hidden":false,"mode":"ojs"},{"id":1572,"value":"// Regular octagon\ngon8 = scpie_eval(await FileAttachment(\"gon8@1.json\").json())\n\n/*\n% Matlab code to generate gon8.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon8 = dmap(exp(2*pi*i*[1:8].' / 8), 0, 1e-14);\n[G8pie, G8pie_data] = scpie(Gon8, 1e-13)\nscpie_json(G8pie_data, 'gon8.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":3678,"value":"// Regular nonagon\ngon9 = scpie_eval(await FileAttachment(\"gon9.json\").json())\n\n/*\n% Matlab code to generate gon9.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon9 = dmap(exp(2*pi*i*[1:9].' / 9), 0, 1e-14);\n[G9pie, G9pie_data] = scpie(Gon9, 1e-13)\nscpie_json(G9pie_data, 'gon9.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":3682,"value":"// Regular decagon\ngon10 = scpie_eval(await FileAttachment(\"gon10.json\").json())\n\n/*\n% Matlab code to generate gon10.json:\n\ndmap = @(vert, cen, tol) ...\n center(diskmap(polygon(vert), scmapopt('Tolerance', tol)), cen);\nGon10 = dmap(exp(2*pi*i*[1:10].' / 10), 0, 1e-14);\n[G10pie, G10pie_data] = scpie(Gon10, 1e-13)\nscpie_json(G10pie_data, 'gon10.json');\n*/","pinned":false,"hidden":false,"mode":"ojs"},{"id":4663,"value":"$49","pinned":false,"hidden":false,"mode":"ojs"},{"id":4665,"value":"$4a","pinned":false,"hidden":false,"mode":"ojs"},{"id":728,"value":"$4b","pinned":false,"hidden":false,"mode":"ojs"},{"id":848,"value":"$4c","pinned":false,"hidden":false,"mode":"ojs"},{"id":5253,"value":"md`\n\n* * *\n\nMatlab (especially older versions) is not my favorite tool for string manipulation.\nBut here’s my extremely hacky code for spitting out JSON from the SCPie data.\n(JSON is not a space-efficient format for big arrays of numerical data, but it is convenient.)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1079,"value":"$4d","pinned":false,"hidden":false,"mode":"ojs"},{"id":1264,"value":"Javacript = md`\n* * *\n\n## Javascript code\n\nOnce we have created an approximation in Matlab, and exported it as JSON, we can evaluate it in Javascript:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":958,"value":"$4e","pinned":true,"hidden":false,"mode":"ojs"},{"id":967,"value":"$4f","pinned":true,"hidden":false,"mode":"ojs"},{"id":1326,"value":"md`\n* * *\n\nHere’s some Javascript code for generating [quasirandom points][] in the disk:\n\n [quasirandom points]: https://observablehq.com/@jrus/plastic-sequence\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1229,"value":"diskpts = {\n const p1 = 0.7548776662466927; // inverse of plastic number\n const p2 = 0.5698402909980532;\n \n return (n) => {\n let squarepts = [...Array(2*n).keys()].map(i => [\n 2 * ((p1 * i) % 1) - 1,\n 2 * ((p2 * i) % 1) - 1]);\n return squarepts.slice(1)\n .filter(([x,y]) => x*x + y*y <= 1)\n .slice(0, n);\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1331,"value":"md`\n* * *\n\n## Miscellaneous imports and helper functions\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":244,"value":"import {animate, math_css, $, $$, responsive_katex} with {$css as $css} from \"@jrus/misc\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":243,"value":"$$css = html`${math_css}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1798,"value":"import {bezeval, bezpts_to_svgpath} from '@jrus/bezplot'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1970,"value":"polypts_to_svgpath = {\n const r = (x, d) => (+x).toFixed(d).replace(/\\.?0+$/,'') // round\n \n return function bezpts_to_svgpath(pts, digits) {\n const d = (digits == null) ? 4 : digits;\n const output = [];\n const n = pts.length;\n output.push(`M${r(pts[0],d)},${r(pts[1],d)}`);\n \n for (let i = 2; i < n; i += 2) {\n output.push(`L${r(pts[i],d)},${r(pts[i+1],d)}`);\n }\n return output.join('\\n');\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1820,"value":"import {circle_interp} from '@jrus/circle-arc-interpolation'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1824,"value":"import {linspace, chunk} from '@jrus/itertools'","pinned":false,"hidden":false,"mode":"ojs"},{"id":2946,"value":"import {rotate} from '@jrus/rotate'","pinned":false,"hidden":false,"mode":"ojs"},{"id":5198,"value":"// Chebyshev nodes of the 2nd kind in [-1, 1]\nchebpts = function chebpts (n) {\n let pts = [], scale = Math.PI * .5 / (n - 1);\n for (let i = n - 1; i >= 0; i--) {\n pts[i] = Math.sin((1 - n + 2*i) * scale);\n }\n return pts;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3906,"value":"$50","pinned":false,"hidden":false,"mode":"ojs"},{"id":2196,"value":"import {srgb_to_xyz, xyz_to_srgb, from_hex, to_hex} from '@jrus/srgb'","pinned":false,"hidden":false,"mode":"ojs"},{"id":2194,"value":"import {Jzazbz, Jzazbz_inverse, JzCzhz, JzCzhz_inverse} from '@jrus/jzazbz'","pinned":false,"hidden":false,"mode":"ojs"},{"id":3095,"value":"// bisect chroma until it is in gamut\nJCh_in_gamut = (jch) => {\n const [J, C, h] = jch;\n let cm = C / 2;\n \n return xyz_to_srgb(JzCzhz_inverse(jch).map(X => .5*X))\n \n // TODO: implement this\n \n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2200,"value":"hex_to_JCh = (hex) => JzCzhz(srgb_to_xyz(from_hex(hex)).map(X => 2*X))","pinned":false,"hidden":false,"mode":"ojs"},{"id":2667,"value":"JCh_to_hex = (jch) => to_hex(xyz_to_srgb(JzCzhz_inverse(jch).map(X => .5*X)))","pinned":false,"hidden":false,"mode":"ojs"},{"id":2664,"value":"JCha_to_rgba = {\n const rgb_to_8bit = function rgb_to_8bit([R, G, B]) {\n R = Math.round(0xff * R), G = Math.round(0xff * G), B = Math.round(0xff * B);\n if (Math.max(R, G, B) > 255 || Math.min(R, G, B) < 0) {\n throw new Error('Bad Input: R, G, and B must be in range [0, 1]'); }\n return [R, G, B];\n }\n \n return (jch, a) => `rgba(${\n rgb_to_8bit(xyz_to_srgb(\n JzCzhz_inverse(jch).map(X => .5*X))).join(',')\n },${a})`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5231,"value":"matlab_file = function (filename, content) {\n\n let URL=`\\nhttps://observablehq.com/@jrus/scpie\\n`\n \n let downloadData = new Blob([`%{${URL}${LICENSE}%}`, content], { type: \"text/plain\" });\n\n const size = (downloadData.size / 1024).toFixed(0);\n const button = DOM.download(\n downloadData,\n filename,\n `Download`\n );\n \n return md`\n### \\`${filename}\\` ${button}\n\n\\`\\`\\`matlab\n${content}\n\\`\\`\\`\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1288,"value":"d3 = require(\"d3-geo@1\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1290,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":2517,"value":"land50 = {\n const file = await fetch(\"https://cdn.jsdelivr.net/npm/world-atlas@2/land-50m.json\");\n const obj = await file.json();\n return topojson.feature(obj, obj.objects.land);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2533,"value":"MathJax = require('https://cdn.jsdelivr.net/npm/mathjax@3.1.0/es5/tex-svg.js')\n .catch(() => window['MathJax'])","pinned":false,"hidden":false,"mode":"ojs"},{"id":3330,"value":"youtube = (id, width, height) =>\n html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":5256,"value":"LICENSE = String.raw`\nCopyright (c) 2020 Jacob Rus\n\nPermission to use, copy, modify, and/or distribute this software for any\npurpose with or without fee is hereby granted, provided that the above\ncopyright notice and this permission notice appear in all copies.\n\nTHE SOFTWARE IS PROVIDED \"AS IS\" AND THE AUTHOR DISCLAIMS ALL WARRANTIES\nWITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF\nMERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR\nANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES\nWHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN\nACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF\nOR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.\n`","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W51","stdlib":"1","titling":"auto"}