30:T7bd,{ const height = a => (a < 0 ? -1 : a > 1 ? 1 : chordAreaInverse(a)); const context = DOM.context2d(440, 330); context.translate(200, 165); d3.select(context.canvas) .on("pointermove", event => { step = 0; const y = 1.05 - d3.pointers(event)[0][1] / 300; a = chordArea(y < -1 ? -1 : y > 1 ? 1 : y); }) .on("pointerout", () => setTimeout(() => (step = .002), 1500)) .style("cursor", "crosshair"); let step = .002; let a = 0; let style; while (true) { if (step > 0) style = d3.hsl(Math.random() * 360, 0.6, 0.7); do { a += step; context.lineWidth = 1; context.save(); context.clearRect(-200, -165, 440, 330); context.beginPath(); context.arc(0, 0, 148, 0, 2 * Math.PI); context.stroke(); context.beginPath(); context.arc(0, 0, 145, 0, 2 * Math.PI); context.clip(); const y = height(a); context.fillStyle = style; context.fillRect(-150, -145 * y, 300, 300); context.restore(); context.lineWidth = .25; context.beginPath(); context.moveTo(-45 * sqrt(1 - y * y), -45 * y); context.lineTo(0, 0); context.lineTo(45 * sqrt(1 - y * y), -45 * y); context.stroke(); const t = 2 * acos(y); context.beginPath(); context.arc(0, 0, 10, -Math.PI / 2 - t / 2, -Math.PI / 2 + t / 2); context.stroke(); context.fillStyle = -145 * y < 80 ? "white" : "black"; context.fillText(`a = ${(Math.max(a, 0) * 100).toFixed(0)}%`, -20, 80); context.fillStyle = -145 * y < 0 ? "white" : "black"; context.fillText(`t = ${((t * 180) / Math.PI).toFixed(0)}°`, 20, 0); context.fillStyle = "black"; context.fillText( `y = ${y.toFixed(2)}`, 35 + 120 * sqrt(1 - y * y), -152 * y ); yield context.canvas; } while (0 <= a && a <= 1); await Promises.delay(500); step = -step; a += 2 * step; } yield context.canvas; }31:T47d,chordAreaInverse2 = { const nodes = [ 0.5, 0, 2.3584407814e-1, 2.0697979931e-2, 2.0894975402e-3, 1.1722274484e-1, 2.1236567413e-4, 4.6305986706e-1, 5.3748971858e-5, 8.1099031276e-3, 5.8743883731e-2], values = [ 1.5874010464, 1.4053918380, 1.4973407222, 1.4208001108, 1.4086435353, 1.4585411281, 1.4060959864, 1.5737389975, 1.4056733185, 1.4135101277, 1.4374270332], weights = [ 2.5236473233e-1, 3.2702933141e-6, 2.1134869853e-1, -4.5296987518e-2, -1.8210616414e-3, -2.2914836584e-1, 2.1728601712e-4, -3.3969424647e-1, -6.8763426681e-5, 1.3707998140e-2, 1.3854310301e-1]; return function chordAreaInverse(a) { var s = (a <= 0.5) - (a > 0.5), b = (a > 0.5) + s * a, c = Math.cbrt(b); return s * (reval(nodes, values, weights, b) * c * c - 1); } function reval(nodes, values, weights, x) { var n = nodes.length, p = 0, q = 0, f, j, xj; for (j = 0; j < n; j++) { q += xj = weights[j] / (x - nodes[j]); p += xj * values[j]; } f = p / q; if (f !== f) { for (j = 0; j < n; j++) { if (x - nodes[j] === 0) { f = values[j]; break; }}} return f; }; }32:[] 2e:{"title":"An approximate formula for the chord area inverse problem","cells":[{"id":0,"value":"md`# An approximate formula for the chord area inverse problem`","pinned":false,"hidden":false,"mode":"ojs"},{"id":566,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":730,"value":"md`_**LEGEND:** The chord area ${tex`a`} is painted in solid color and gets filled (or emptied) at a constant speed. The angle ${tex`t`} is indicated by the compas. ${tex`y`} is the vertical coordinate of the chord relative to the center._`","pinned":false,"hidden":false,"mode":"ojs"},{"id":228,"value":"md`The chord area inverse problem considers the following question: this tank (which, since we’re in two dimensions, is a disc) contains a certain quantity of water; what level ${tex`y`} does the water reach?\n\nWe’ll consider ${tex`y=0`} to be the center of the disc, with ${tex`y=1`} at the top (full tank), and ${tex`-1`} the level of the empty tank.\n\nIt can be demonstrated that the relative area ${tex`a`} under the chord is:\n\n> ${tex`a = {\\small \\dfrac{1}{\\pi}}\\Bigl( \\arccos(-y) + y\\sqrt{1-y^2} \\Bigr)`}.\n\nTo answer our initial question (given ${tex`a`}, what is the value of ${tex`y`}?), we want the *inverse* of this function. However, there exists [no closed-form solution](https://en.wikipedia.org/wiki/Kepler's_equation). We are left to create an approximations via polynoms or iterative algorithms.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":293,"value":"md`**Here I suggest that we compute:**\n> ${tex`\\bm {y \\approx a^{3/5} - (1-a)^{3/5}}`}\n\nThis approximation appears to have a maximal error of 0.17% and a median error of 0.08% over the domain. Very crude, but simple and fast enough [for some purposes](/@fil/disc-transport).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":125,"value":"chordAreaInverse = a => Math.pow(a, .6) - Math.pow(1 - a, .6)","pinned":true,"hidden":false,"mode":"ojs"},{"id":945,"value":"md`**Alternative method.** If you’re looking for a more accurate method (9 decimal digits), [Jacob Rus](/@jrus) suggested the following, which is the product of a power function ${tex`a^{2/3}`} and a minimax rational approximation ${tex`r(a)`},\n> ${tex`y \\approx \\begin{cases}\n-1 + a^{2/3}r(a), & \\text{if } 0 \\leq a \\leq \\tfrac{1}{2} \\\\[.2em]\n1 - (1-a)^{2/3}r(1-a), & \\text{if } \\tfrac{1}{2} \\lt a \\leq 1\n\\end{cases}`}\n\nusing [Chebfun’s \\`aaa\\` tool](https://arxiv.org/pdf/1908.06001) to find the coefficients for ${tex`r`}. If this looks like high art, it’s because it is:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":944,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":601,"value":"md`----\n\n_The following are various computations to validate both approximations._`","pinned":false,"hidden":false,"mode":"ojs"},{"id":967,"value":"time(chordAreaInverse)","pinned":true,"hidden":false,"mode":"ojs"},{"id":969,"value":"time(chordAreaInverse2)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1007,"value":"error(chordAreaInverse)","pinned":true,"hidden":false,"mode":"ojs"},{"id":986,"value":"error(chordAreaInverse2)","pinned":true,"hidden":false,"mode":"ojs"},{"id":965,"value":"time = func => {\n const N = 1e6;\n const a = Array.from({length: N}, Math.random);\n const t = performance.now();\n for (let i = 0; i < N; i++) func(a[i]);\n return md`${Math.round(performance.now() - t)} μs per evaluation`;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1011,"value":"// worst case round-trip error\nerror = func => {\n const N = 1e6;\n let err = 0;\n for (let i = 0; i <= N; i++) {\n const a = i / N;\n err = Math.max(Math.abs(chordArea(func(a) / 2 + 0.5) - a), err);\n }\n return err;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":808,"value":"md`*Note:* all the solutions I could find online work with the angle ${tex`t`}, formed by a compass that joins the center of the circle to the chord’s intersections with the circle: ${tex`t = 2 \\cos^{-1}(y)`}. The area is then expressed as ${tex`a = \\frac1{2π}(t - \\sin t) `}. Whatever the notation, the approximation is a more or less involved formula —and some with a great precision—, but that risks to be slow.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":329,"value":"A = th => (th - sin(th)) / π / 2 // 0 to 1","pinned":true,"hidden":false,"mode":"ojs"},{"id":336,"value":"chordArea = x => A(2 * acos(1 - 2 * x))","pinned":true,"hidden":false,"mode":"ojs"},{"id":333,"value":"π = Math.PI","pinned":true,"hidden":false,"mode":"ojs"},{"id":341,"value":"chordArea(.5)","pinned":true,"hidden":false,"mode":"ojs"},{"id":354,"value":"DIFF = (x, e) => Math.abs(x - (1 + approx2(chordArea(x), e)) / 2)","pinned":true,"hidden":false,"mode":"ojs"},{"id":455,"value":"approx2 = (t, e) => Math.pow(t, e) - Math.pow(1 - t, e)","pinned":true,"hidden":false,"mode":"ojs"},{"id":415,"value":"d3.max(\n Float64Array.from({ length: 100000 }, (_, i) => DIFF(i / (100000 - 1), 0.6))\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":488,"value":"d3.median(\n Float64Array.from({ length: 100000 }, (_, i) => DIFF(i / (100000 - 1), 0.6))\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":296,"value":"md`The curves:\n- in blue, a count of points\n- in red, the function above\n- in black, the absolute error (magnified ×500)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":107,"value":"{\n const height = 500;\n const svg = d3.create(\"svg\").attr(\"viewBox\", [0, 0, width, height]);\n\n svg\n .append(\"path\")\n .attr(\"fill\", \"none\")\n .attr(\"stroke\", \"steelblue\")\n .attr(\"stroke-width\", \"3\")\n .attr(\n \"d\",\n d3.line()(\n Array.from({ length: 400 }, (_, i) => [\n (width * i) / 399,\n (height * (1 - CDF.invert(i / 399))) / 2\n ])\n )\n );\n\n svg\n .append(\"path\")\n .attr(\"fill\", \"none\")\n .attr(\"stroke\", \"red\")\n .attr(\n \"d\",\n d3.line()(\n Array.from({ length: 400 }, (_, i) => [\n (width * i) / 399,\n (height * (1 - chordAreaInverse(i / 399))) / 2\n ])\n )\n );\n\n svg\n .append(\"path\")\n .attr(\"fill\", \"none\")\n .attr(\"stroke\", \"black\")\n .attr(\n \"d\",\n d3.line()(\n Array.from({ length: 400 }, (_, i) => [\n (width * i) / 399,\n height * (1 - 500 * DIFF(i / 399, 0.6))\n ])\n )\n );\n\n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":302,"value":"md`\n\n---\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"f = x => sqrt(1 - x * x)","pinned":true,"hidden":false,"mode":"ojs"},{"id":43,"value":"CDF = integral(f, -1, 1)","pinned":true,"hidden":false,"mode":"ojs"},{"id":60,"value":"CDF(.3)","pinned":true,"hidden":false,"mode":"ojs"},{"id":62,"value":"CDF.invert(0.6893155371739725)","pinned":true,"hidden":false,"mode":"ojs"},{"id":75,"value":"CDF(0.1)","pinned":true,"hidden":false,"mode":"ojs"},{"id":87,"value":"CDF.invert(CDF(0.1))","pinned":true,"hidden":false,"mode":"ojs"},{"id":15,"value":"function integral(f, a, b) {\n const n = 10000;\n const step = (1 / (n - 1)) * (b - a);\n const v = d3.cumsum(\n Float32Array.from({ length: n }, (_, i) => f(a + i * step) * step)\n );\n const C = 1 / (v[v.length - 1] - v[0]);\n for (let i = 0; i < n; ++i) v[i] = (v[i] - v[0]) * C;\n const g = d3.scaleSequential(d3.piecewise(v)).domain([a, b]);\n g.invert = x => {\n const i = d3.bisect(v, x),\n delta = (x - v[i]) / (v[i + 1] - v[i]);\n return ((i + (isNaN(delta) ? 0 : delta)) / (n - 1)) * (b - a) + a;\n };\n return g;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":7,"value":"import { acos, cos, sin, sqrt } from \"@fil/math\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":16,"value":"d3 = require(\"d3@6\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":837,"value":"md`And let’s end with a useless disc:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":219,"value":"{\n const context = DOM.context2d(500, 500);\n context.translate(250, 250);\n\n context.beginPath();\n context.arc(0, 0, 248, 0, 2 * Math.PI);\n context.stroke();\n for (let x = -1; x <= 1; x += .017) {\n const A = sqrt(1 - x * x);\n for (let y = -A; y <= A; y += .017) {\n context.beginPath();\n context.arc(x * 245, y * 245, 1, 0, 2 * Math.PI);\n context.fill();\n }\n yield context.canvas;\n }\n yield context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W32","stdlib":"1","titling":"auto"}