31:T418,performSimulation = { // Trigger redraw on restart or zoombox change restart; zoomBox; const run = ~simulate.indexOf('Simulate'); const showPen = ~simulate.indexOf('Show pendulum at point under mouse'); let needsDraw = true; const frame = regl.frame(() => { // Configure a full-viewport shader blit(() => { // Step the simulation if running if (run) { // Run two steps from one data framebuffer to the other, // i.e. from y1=f(y0) and then back to y0=f(y1) update([{ dt, src: y[0], dst: y[1] }, { dt, src: y[1], dst: y[0] }]); } // Draw the simulation if something has happened if (needsDraw) { drawMap({ y: y[0] }); if (showPen && mouseHover.active && !mouseHover.boxActive) { drawPendulum([ { y, w1: [0, 0, 0], w2: [0, 1, 0] }, { y, w1: [0, 1, 0], w2: [0, 0, 1] } ]); } } // Trigger update needsDraw = run ? true : false; }); }); invalidation.then(() => frame.cancel()); }32:T6aa,update = regl({ frag: ` precision highp float; uniform sampler2D y; uniform vec2 resolution; uniform float dt; #define PI ${Math.PI} // From: https://en.wikipedia.org/wiki/Double_pendulum#Lagrangian vec4 derivative (vec4 state) { vec2 theta = state.xy; vec2 pTheta = state.zw; float threeCosTheta12 = 3.0 * cos(theta.x - theta.y); vec2 thetaDot = 6.0 * ( vec2( 2.0 * pTheta.x - threeCosTheta12 * pTheta.y, 8.0 * pTheta.y - threeCosTheta12 * pTheta.x ) / (16.0 - threeCosTheta12 * threeCosTheta12) ); float thetaDot12sinTheta12 = thetaDot.x * thetaDot.y * sin(theta.x - theta.y); vec2 pThetaDot = -0.5 * vec2( thetaDot12sinTheta12 + 3.0 * sin(theta.x), -thetaDot12sinTheta12 + sin(theta.y) ); return vec4(thetaDot, pThetaDot); } void main () { // Read the state (theta1, theta2, velocity1, velocity2) from the rgba texture vec4 yn = texture2D(y, gl_FragCoord.xy / resolution); // Convert from [0, 1] to [-PI, PI]: yn = yn * (2.0 * PI) - PI; // RK4 integration vec4 k1 = dt * derivative(yn); vec4 k2 = dt * derivative(yn + 0.5 * k1); vec4 k3 = dt * derivative(yn + 0.5 * k2); vec4 k4 = dt * derivative(yn + k3); yn += (k1 + k4 + 2.0 * (k2 + k3)) / 6.0; // Convert back from [-PI, PI] to [0, 1]: yn = (yn + PI) / (2.0 * PI); // Loop angles if they exceed the range yn.xy = fract(yn.xy); // Output the state to the four color channels gl_FragColor = yn; }`, framebuffer: regl.prop('dst'), uniforms: { y: regl.prop('src'), dt: regl.prop('dt') } })33:T69a,drawPendulum = regl({ vert: ` precision mediump float; attribute vec2 line; uniform vec2 resolution; uniform vec3 w1, w2; uniform vec2 hover; uniform float lineWidth; uniform sampler2D y; #define PI ${Math.PI} void main () { // Sample the state directly from the texture by passing the hovered // point as a uniform, so that the data never leaves the GPU. vec2 state = texture2D(y, hover).xy * (2.0 * PI) - PI; // Vectors representing the directions of the arms of the pendulum vec2 v1 = vec2(sin(state.x), -cos(state.x)); vec2 v2 = vec2(sin(state.y), -cos(state.y)); // Points representing the joints of the pendulum (and the end) vec2 p0 = vec2(0); vec2 p1 = v1; vec2 p2 = p1 + v2; // Select two of the three above (arm1 or arm2) as endpoints of // the line to be rendered vec2 a = w1.x * p0 + w1.y * p1 + w1.z * p2; vec2 b = w2.x * p0 + w2.y * p1 + w2.z * p2; // Normal and tangential vectors vec2 t = normalize(b - a); vec2 n = vec2(-t.y, t.x); // Select which end of line we're at vec2 p = mix(a, b, line.y) * 0.5; // Offset normal to the line by the line width p += n * line.x * lineWidth; gl_Position = vec4(p, 0, 1); }`, frag: ` precision lowp float; void main () { gl_FragColor = vec4(0, 0, 0, 1); }`, attributes: { line: [0.5, 0, 0.5, 1, -0.5, 0, -0.5, 1] }, uniforms: { y: regl.prop('y[0]'), w1: regl.prop('w1'), w2: regl.prop('w2'), hover: () => mouseHover.pt, lineWidth: ctx => (4.0 * ctx.pixelRatio) / ctx.framebufferWidth, resolution: ctx => [ctx.framebufferWidth, ctx.framebufferHeight] }, depth: { enable: false }, primitive: 'triangle strip', count: 4 })34:Tdca,viewof zoomBox = { viz; const dummy = html`

zoomBox

`; const zoomBox = { center: [0, 0], scale: [Math.PI * 2, Math.PI * 2] }; const canvas = regl._gl.canvas; let start = { x: NaN, y: NaN }; let end = { x: NaN, y: NaN }; let captured = false; function emitUpdate() { dummy.dispatchEvent(new CustomEvent('input')); } function onDoubleClick() { zoomBox.center[0] = 0; zoomBox.center[1] = 0; zoomBox.scale[0] = Math.PI * 2; zoomBox.scale[1] = Math.PI * 2; emitUpdate(); } canvas.addEventListener('dblclick', onDoubleClick, false); const box = document.getElementById('zoombox'); box.style.position = 'absolute'; box.style.border = '1px solid rgba(0,0,0,0.6)'; box.style.backgroundColor = 'rgba(255,255,255,0.2)'; box.style.pointerEvents = 'none'; function updateBox() { box.style.left = `${Math.min(start.x, end.x)}px`; box.style.top = `${Math.min(start.y, end.y)}px`; box.style.width = `${Math.abs(end.x - start.x)}px`; box.style.height = `${Math.abs(end.y - start.y)}px`; } function onPointerDown(event) { captured = true; start.x = event.offsetX; start.y = event.offsetY; end.x = event.offsetX; end.y = event.offsetY; canvas.setPointerCapture(event.pointerId); box.style.display = 'block'; updateBox(); mouseHover.boxActive = true; } function onPointerMove(event) { mouseHover.pt[0] = event.offsetX / dimension; mouseHover.pt[1] = 1 - event.offsetY / dimension; mouseHover.boxActive = false; if (!captured) return; end.x = event.offsetX; end.y = event.offsetY; updateBox(); mouseHover.boxActive = true; } function onPointerUp(event) { box.style.display = 'none'; captured = false; end.x = event.offsetX; end.y = event.offsetY; if (end.x === start.x || end.y === start.y) return; let curX0 = zoomBox.center[0] - zoomBox.scale[0] * 0.5; let curX1 = zoomBox.center[0] + zoomBox.scale[0] * 0.5; let curY0 = zoomBox.center[1] - zoomBox.scale[1] * 0.5; let curY1 = zoomBox.center[1] + zoomBox.scale[1] * 0.5; let xcen = (0.5 * (start.x + end.x)) / dimension; let ycen = (0.5 * (start.y + end.y)) / dimension; let xscale = Math.abs(end.x - start.x) / dimension; let yscale = Math.abs(end.y - start.y) / dimension; zoomBox.center[0] = lerp(curX0, curX1, xcen); zoomBox.center[1] = lerp(curY0, curY1, ycen); zoomBox.scale[0] *= xscale; zoomBox.scale[1] *= yscale; canvas.releasePointerCapture(event.pointerId); updateBox(); emitUpdate(); mouseHover.boxActive = false; } function onPointerEnter() { mouseHover.active = true; } function onPointerLeave() { mouseHover.active = false; } canvas.addEventListener('pointerdown', onPointerDown, false); canvas.addEventListener('pointermove', onPointerMove, false); canvas.addEventListener('pointerup', onPointerUp, false); canvas.addEventListener('pointerenter', onPointerEnter, false); canvas.addEventListener('pointerleave', onPointerLeave, false); invalidation.then(() => { canvas.addEventListener('dblclick', onDoubleClick, false); canvas.removeEventListener('pointerdown', onPointerDown); canvas.removeEventListener('pointermove', onPointerMove); canvas.removeEventListener('pointerup', onPointerUp); canvas.removeEventListener('pointerenter', onPointerEnter); canvas.removeEventListener('pointerleave', onPointerLeave); }); dummy.value = zoomBox; return dummy; }35:[{"id":"d20dde047e5675c3dc5d6ce17a4bd428f86de6ffdc8c88d8d22db8ff75495ef712b8065ef6d15f53bc9ffa71b97558adcaf0b0f2b47ae332f6627b6b5351fecd","href":"https://static.observableusercontent.com/files/d20dde047e5675c3dc5d6ce17a4bd428f86de6ffdc8c88d8d22db8ff75495ef712b8065ef6d15f53bc9ffa71b97558adcaf0b0f2b47ae332f6627b6b5351fecd","name":"1920px-Double-compound-pendulum-dimensioned.svg.png","type":"image/png","lastModified":1620968372572.889,"size":210679}] 2f:{"title":"The Double Pendulum Map","cells":[{"id":1024,"value":"# The Double Pendulum Map","pinned":false,"hidden":false,"mode":"md"},{"id":1031,"value":"*[Read this notebook ➡️](https://rreusser.github.io/notebooks/double-pendulum-map/) on my personal website via [Observable Notebook Kit](https://observablehq.com/notebook-kit/kit).*","pinned":false,"hidden":false,"mode":"md"},{"id":539,"value":"(*Edit: Previous title was *The Double Pendulum Fractal.* [Jacob Rus](https://observablehq.com/@jrus) has, I think correctly, questioned the extent to which this meets the formal definition of a fractal. I'm not certain, so I'll call it a map.*)\n\nA [double pendulum](https://en.wikipedia.org/wiki/Double_pendulum) consists of two connected pendulums, rotating under the force of gravity. A compound double pendulum has the weight distributed along the length of the arms rather than as point masses at the ends of strings.","pinned":false,"hidden":false,"mode":"md"},{"id":530,"value":"html`

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":363,"value":"md`[Sam Maksimovich](https://twitter.com/MaksimovichSam) recently created a lovely video in which he simultaneously plots all possible double pendulums starting with zero velocity, using color to illustrate their trajectories. The explanation is far better than anything you'll find here. You should watch it!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":526,"value":"tweet('1392357775666159618')","pinned":false,"hidden":false,"mode":"ojs"},{"id":547,"value":"md`At the end he mentions it took six hours to render the video. I'm confident I would not be able to create a lovelier or more effective explanation video, but the visualization is exactly the sort of thing GPUs do very well, so it seemed like it would be fun to plug this into WebGL and make a quick, interactive visualization of the pendulum. A major shortcoming is that WebGL only does single precision, and chaotic simulations definitely benefit from at least double precision, but maybe we'll see what we can get away with.\n\nThe visualization below uses the [regl](https://github.com/regl-project/regl) WebGL library to compute and visualize the angle of the of the first (x-axis) and second (y-axis) arms of all compound double pendulums starting with the specified velocities. This amounts to a two-dimensional slice (just the angular positions) of the four-dimensional phase space (two positions and two velocities).\n\nAdjust the initial velocities and zoom in on regions below and observe…… pretty colors! 🌈\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":610,"value":"{\n if (colorType === 'half float') {\n return html`

Warning! This simulation is using half-float precision and as a result won't be accurate. Try a device which supports the floating point WebGL textures for better results.

`;\n } else {\n return html``;\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":573,"value":"viewof omega0_0 = Range([-2, 2], {\n label: 'Initial angular velocity of arm 1',\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":576,"value":"viewof omega0_1 = Range([-2, 2], {\n label: 'Initial angular velocity of arm 2',\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":313,"value":"viewof dt = Range([0.001, 0.1], {\n label: 'Time step',\n value: 0.03,\n transform: Math.log\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":281,"value":"viewof simulate = Checkbox(['Simulate', 'Show pendulum at point under mouse'], {\n value: ['Simulate', 'Show pendulum at point under mouse']\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":300,"value":"viewof restart = Button('Restart')","pinned":false,"hidden":false,"mode":"ojs"},{"id":289,"value":"viz = html`

\n${regl._gl.canvas}\n\n

Hover to view pendulum at current point. Click and drag to zoom. Double click to return to initial view

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":814,"value":"md`${tex`\\theta_1 \\in [${(zoomBox.center[0] - zoomBox.scale[0] / 2).toFixed(\n 8\n)}, ${(zoomBox.center[0] + zoomBox.scale[0] / 2).toFixed(\n 8\n)}]`}
${tex`\\theta_2 \\in [${(\n zoomBox.center[1] -\n zoomBox.scale[1] / 2\n).toFixed(8)}, ${(zoomBox.center[1] + zoomBox.scale[1] / 2).toFixed(8)}]`}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":809,"value":"md`## Implementation`","pinned":false,"hidden":false,"mode":"ojs"},{"id":554,"value":"md`The code below sets up and runs the simulation. We start by initializing a regl WebGL context.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"regl = {\n var canvas = document.createElement('canvas');\n let regl = createREGL({\n optionalExtensions: [\n 'OES_texture_half_float',\n 'OES_texture_float',\n 'ANGLE_instanced_arrays'\n ],\n canvas,\n pixelRatio: pixelRatio,\n attributes: {\n antialias: true,\n alpha: false,\n depthStencil: false,\n preserveDrawingBuffer: false\n }\n });\n return regl;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"resizeContext = {\n // Depend on page width here to avoid recreating regl context when the page\n // width changes\n const canvas = regl._gl.canvas;\n var w = radius;\n var h = radius;\n canvas.style.width = `${Math.floor(w / pixelRatio)}px`;\n canvas.style.height = `${Math.floor(h / pixelRatio)}px`;\n canvas.width = w;\n canvas.height = h;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":556,"value":"md`The main simulation and initialization setup.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":158,"value":"performInitialization = {\n zoomBox;\n \n // Trigger reinitialization when 'restart' clicked\n restart;\n\n // Must call this when calling regl outside of frame loop\n regl.poll();\n\n // Call initialize in blit context, writing to y[0] framebuffer\n y[0].use(() =>\n blit(() => {\n initialize({ omega0_0, omega0_1, zoomBox });\n })\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":54,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":560,"value":"md`Configure two 2D textures to store (theta1, theta2, velocity1, velocity2) in the four channels of an (r, g, b, a) texture.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":78,"value":"colorType = canWriteToFBOOfType(regl, 'float')\n ? 'float'\n : canWriteToFBOOfType(regl, 'half float')\n ? 'half float'\n : 'uint8'","pinned":false,"hidden":false,"mode":"ojs"},{"id":58,"value":"y = {\n const fbos = [0, 1].map(() =>\n regl.framebuffer({\n colorType,\n radius,\n depth: false,\n stencil: false\n })\n );\n invalidation.then(() => fbos.forEach(fbo => fbo.destroy()));\n return fbos;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":564,"value":"md`Construct shader commands to initialize, update, and draw.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":113,"value":"blit = regl({\n // Configure a full-viewport vertex shader:\n vert: `\n precision mediump float;\n attribute vec2 xy;\n void main () {\n gl_Position = vec4(xy, 0, 1);\n }`,\n attributes: { xy: [-4, -4, 4, -4, 0, 4] },\n uniforms: {\n resolution: ctx => [ctx.framebufferWidth, ctx.framebufferHeight]\n },\n count: 3,\n depth: { enable: false }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":151,"value":"initialize = regl({\n frag: `\n precision highp float;\n uniform sampler2D colorscale;\n uniform vec2 resolution, omega0, center, scale;\n #define PI ${Math.PI}\n void main () {\n // Really just the full spectrum of angles for the two parts of the pendulum,\n // [0, 2pi] x [0, 2pi], with zero velocity, but we center it about 0.5 to\n // slightly improve how well it works with half float precision \n vec2 theta0 = center + (0.5 + scale * ((vec2(0, 1) + vec2(1, -1) * gl_FragCoord.xy / resolution) - 0.5));\n gl_FragColor = (vec4(theta0 * (2.0 * PI) - PI, omega0) + PI) / (2.0 * PI);\n }`,\n uniforms: {\n omega0: (ctx, props) => [props.omega0_0, props.omega0_1],\n center: (ctx, props) => props.zoomBox.center.map(x => (x * 0.5) / Math.PI),\n scale: (ctx, props) => [\n props.zoomBox.scale[0] / (Math.PI * 2),\n props.zoomBox.scale[1] / (Math.PI * 2)\n ]\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":176,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":95,"value":"drawMap = regl({\n frag: `\n precision highp float;\n uniform sampler2D colorscale, y;\n uniform vec2 resolution;\n #define TWOPI ${2 * Math.PI}\n \n vec3 colormap(vec2 uv) {\n uv *= TWOPI;\n return pow(0.5 + 0.5 * vec3(sin(uv.x) * vec2(-cos(uv.y), sin(uv.y)), -cos(uv.x)), vec3(0.75));\n }\n\n void main () {\n // Sample the state at this pixel\n vec4 state = texture2D(y, gl_FragCoord.xy / resolution);\n\n // Color it by the projection of the 4D state to 2D (theta1, theta2):\n gl_FragColor = vec4(colormap(state.xy), 1);\n }`,\n uniforms: {\n y: regl.prop('y')\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":886,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":889,"value":"mouseHover = ({ pt: [0, 0], active: false, boxActive: false })","pinned":false,"hidden":false,"mode":"ojs"},{"id":619,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":518,"value":"md`## Other setup and imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":766,"value":"dimension = width","pinned":false,"hidden":false,"mode":"ojs"},{"id":82,"value":"radius = Math.floor(dimension * pixelRatio) // (height/width in regl-speak)","pinned":false,"hidden":false,"mode":"ojs"},{"id":49,"value":"pixelRatio = window.devicePixelRatio","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"createREGL = require('regl')","pinned":false,"hidden":false,"mode":"ojs"},{"id":485,"value":"import { canWriteToFBOOfType } from '@rreusser/regl-tools'","pinned":false,"hidden":false,"mode":"ojs"},{"id":275,"value":"import { Button, Checkbox, Range } from '@observablehq/inputs'","pinned":false,"hidden":false,"mode":"ojs"},{"id":523,"value":"import { tweet } from '@mbostock/tweet'","pinned":false,"hidden":false,"mode":"ojs"},{"id":635,"value":"import { rangeSlider } from '@mootari/range-slider'","pinned":false,"hidden":false,"mode":"ojs"},{"id":791,"value":"function lerp(a, b, t) {\n return a + t * (b - a);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":651,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}