1d:{"title":"Periodic Planar Three-Body Orbits","cells":[{"id":8869,"value":"viewof title = md`\n# Periodic Planar Three-Body Orbits`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":9023,"value":"*[Read this notebook ➡️](https://rreusser.github.io/notebooks/periodic-planar-three-body-orbits/) on my personal website via [Observable Notebook Kit](https://observablehq.com/notebook-kit/kit).*","pinned":false,"hidden":false,"mode":"md"},{"id":1410,"value":"$1f","pinned":false,"hidden":false,"mode":"ojs"},{"id":8536,"value":"$20","pinned":false,"hidden":false,"mode":"ojs"},{"id":2902,"value":"$21","pinned":false,"hidden":false,"mode":"ojs"},{"id":5971,"value":"viewof integration = md`## Time integration`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1359,"value":"$22","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"import {ode45} from '@rreusser/integration'","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"trajectory = {\n let config = {\n tolerance: initialConditions.tolerance || 1e-9,\n tLimit: initialConditions.period\n };\n\n // An object to store the current state of the integration\n let state = { t: 0, y: initialConditions.y0.slice() };\n\n // An object into which to accumulate the full history\n let result = { position: [], t: [] };\n\n function storeStep(t, y) {\n result.t.push(t);\n // We store (x, y, t) for each step. That is, we extend time into\n // the z-dimension. This allows us to plot in the 2D plane and use\n // built-in GPU clipping to slide the orbit in the z-dimension and\n // animate the lines.\n result.position.push(y[0], y[1], t, y[4], y[5], t, y[8], y[9], t);\n }\n\n // Store the initial conditions\n storeStep(state.t, state.y);\n\n let step = 0;\n while (step++ < 1e6 && !state.limitReached) {\n // Perform a single adaptive step of arbitrary size\n ode45(state, planarThreeBodyDerivative, config);\n\n storeStep(state.t, state.y);\n }\n\n return result;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5963,"value":"md`## Initial conditions`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3695,"value":"$23","pinned":false,"hidden":false,"mode":"ojs"},{"id":2534,"value":"md`The initial conditions in this notebook are obtained from Milovan Šuvakov's [Three-body Gallery](http://three-body.ipb.ac.rs/) as well as from the listing of [Xiaoming Li and Shijun Liao](http://numericaltank.sjtu.edu.cn/three-body/three-body.htm). The initial conditions presented here fall into a few couple categories but all have zero net translational and angular momentum.\n\nFor more information as well as an excellent and approachable overview of the techniques used to locate the orbits, see Šuvakov and Dmitrašinović's [A guide to hunting periodic three-body orbits](https://ui.adsabs.harvard.edu/abs/2014AmJPh..82..609S).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5870,"value":"import {\n threeBodyInitialConditions,\n references,\n sequenceUrls\n} from '@rreusser/periodic-three-body-initial-conditions'","pinned":false,"hidden":false,"mode":"ojs"},{"id":8911,"value":"$24","pinned":false,"hidden":false,"mode":"ojs"},{"id":4364,"value":"$25","pinned":false,"hidden":false,"mode":"ojs"},{"id":105,"value":"$26","pinned":false,"hidden":false,"mode":"ojs"},{"id":3198,"value":"$27","pinned":false,"hidden":false,"mode":"ojs"},{"id":8330,"value":"$28","pinned":false,"hidden":false,"mode":"ojs"},{"id":2241,"value":"$29","pinned":false,"hidden":false,"mode":"ojs"},{"id":5845,"value":"$2a","pinned":false,"hidden":false,"mode":"ojs"},{"id":7196,"value":"$2b","pinned":false,"hidden":false,"mode":"ojs"},{"id":7445,"value":"$2c","pinned":false,"hidden":false,"mode":"ojs"},{"id":5836,"value":"function computeShape (out, x1, y1, x2, y2, x3, y3) {\n let m1 = initialConditions.m[0]\n let m2 = initialConditions.m[1]\n let m3 = initialConditions.m[2]\n\n // Compute Jacobi coordinates\n let r1x = (x1 - x2) / Math.sqrt(2)\n let r1y = (y1 - y2) / Math.sqrt(2)\n let r2x = ((m1 * x1 + m2 * x2) / (m1 + m2) - x3) * 2 / Math.sqrt(6)\n let r2y = ((m1 * y1 + m2 * y2) / (m1 + m2) - y3) * 2 / Math.sqrt(6)\n let r1sq = r1x * r1x + r1y * r1y\n let r2sq = r2x * r2x + r2y * r2y\n\n // Convert to shape space\n out[2] = 2.0 * (r1x * r2y - r1y * r2x) / (r1sq + r2sq)\n out[0] = 2.0 * (r1x * r2x + r1y * r2y) / (r1sq + r2sq)\n out[1] = (r2sq - r1sq) / (r1sq + r2sq)\n return out\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":7647,"value":"md`We return to the ${tex`${initialConditions.name}`} orbit and [plot it below](#projectedShapeSphere) on the shape sphere.\n\nSince all orbits here are collision-free, the trajectory meanders about on the surface of the sphere but never passes through any of the three pairwise collision points. In fact we can imagine continuously deforming a trajectory to find similar nearby orbits just as long as we don't ever pass over the collision points and create an invalid orbit. In topology, these forbidden missing points are called *punctures*, and a family of orbits we can obtain through continuous deformation without passing over the punctures corresponds the concept of the [*fundamental group*](https://en.wikipedia.org/wiki/Fundamental_group) from algebraic topology. Algebraic topology is a topic well beyond my expertise, so I'll tread lightly here and not pretend to know more about it.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5627,"value":"shape = {\n let result = []\n let y = trajectory.position;\n for (var i = 0, i9 = 0; i < trajectory.t.length; i++, i9 += 9) {\n result.push.apply(result, computeShape([],\n y[i9 + 0], y[i9 + 1],\n y[i9 + 3], y[i9 + 4],\n y[i9 + 6], y[i9 + 7],\n ))\n }\n return result\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5856,"value":"function computeStereographicProjection (out, xyz) {\n out[0] = xyz[0] / (1 - xyz[1]) * 2.0\n out[1] = -1\n out[2] = xyz[2] / (1 - xyz[1]) * 2.0\n return out\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5860,"value":"projectedShape = {\n let result = {xyt: [], xyz: []}\n let xy = []\n let n = shape.length\n for (let i = 0, i3 = 0; i3 < n; i++, i3 += 3) {\n let xy = computeStereographicProjection([], [shape[i3], shape[i3 + 1], shape[i3 + 2]])\n result.xyz.push(xy[0], xy[1], xy[2])\n \n let fisheye = computeFisheye([], [xy[0], xy[2]], shapeDrawing.fisheye);\n result.xyt.push(fisheye[0], fisheye[1], trajectory.t[i])\n }\n return result\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6271,"value":"$2d","pinned":false,"hidden":false,"mode":"ojs"},{"id":7684,"value":"md`Three dimensional space is still a bit much to work with, so we additionally turn the shape sphere on its side, select one of the punctures, and project the orbit stereographically about that point onto a plane. The figure above illustrates this process.\n\nThe figure below shows the orbit projected into the two-dimensional plane. Two punctures are visible while the the puncture about which we projected is removed to infinity. The shape trajectory now lies in a flat two dimensional plane, looping around the two projected punctures but never passing through them.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5467,"value":"$2e","pinned":false,"hidden":false,"mode":"ojs"},{"id":7797,"value":"viewof animateShape = checkbox({options: [{value: 'animate', label: 'animate shape trajectory'}]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":5505,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs"},{"id":7948,"value":"md`From here, our task is conceptually simple. We walk along the projected shape of the orbit and note every time it completes a loop, either clockwise or counterclockwise, around one of the two punctures. The result is a [_free group word_](https://en.wikipedia.org/wiki/Free_group) and (at last!) identifies our orbit in the topological sense. We record an ${tex`a`} for clockwise loops around the righthand puncture and ${tex`a^{-1}`} (which we abbreviate as ${tex`A`}) for counterclockwise loops. Similarly we record ${tex`b`} and ${tex`B`}, respectively, for counterclockwise and clockwise loops around the lefthand puncture.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":8418,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":8714,"value":"md`Though conceptually simple, reading these elements is somewhat tedious. In their paper [A guide to hunting periodic three-body orbits](https://ui.adsabs.harvard.edu/abs/2014AmJPh..82..609S), Šuvakov and Dmitrašinović present a \"free group word reading algorithm\". Following the algorithm, we track every crossing of the shape sphere equator, which corresponds to collinear configurations (*syzygies*). We track the direction of each crossing and which body is in the middle, then convert the information via a lookup table into the free group word.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":8207,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":8728,"value":"md`In this case, for the ${tex`${initialConditions.name}`} orbit, we compute the free group word *${freeGroupWord}*.\n\nMany patterns lie within. Figure-eight orbits, for example, follow a repeating *BabA* pattern. Other classes of orbits follow more complicated patterns. To preempt an issue you may have noticed, *yes*, there is indeed ambiguity in this classification. We could have chosen one of the other puncture points about which to project. For some families of orbits this makes a difference; for others it does not. Orbits are similarly equivalent under time reversal so that the free group word *Ab* is equivalent to *Ba*. Šuvakov and Dmitrašinović have a much more detailed discussion of these considerations in their paper [Linear stability of periodic three-body orbits with zero angular momentum and topological dependence of Kepler's third law: a numerical test](https://arxiv.org/abs/1705.03728). It's worth a read!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":8522,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":8848,"value":"md`## License\n\nWith the exception of code imported from other sources as stated above, the code in this notebook is MIT Licensed.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":8838,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":8861,"value":"md`The written content of this notebook is licensed under a ${html`Creative Commons Attribution-NonCommercial 4.0 International License
$\"Creative$ .`}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1474,"value":"md`## Imports and definitions`","pinned":false,"hidden":false,"mode":"ojs"},{"id":7530,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":244,"value":"d3 = require('d3@5')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1248,"value":"createREGL = require('regl')","pinned":false,"hidden":false,"mode":"ojs"},{"id":6489,"value":"import {createIcosphere} from '@rreusser/icosphere'","pinned":false,"hidden":false,"mode":"ojs"},{"id":241,"value":"import {\n getOrAttachReglFO,\n constrainLinearScaleAspectRatio,\n createReglViewportConfiguration,\n createReglLinearScaleConfiguration,\n persistentZoom,\n viewportAxes,\n createReglCamera,\n createInteractions\n} from '@rreusser/regl-tools'","pinned":false,"hidden":false,"mode":"ojs"},{"id":7785,"value":"import {vec3normalize, vec3scaleAndAdd, vec3create} from '@rreusser/gl-vec3'","pinned":false,"hidden":false,"mode":"ojs"},{"id":5074,"value":"import {\n mat4copy,\n mat4create,\n mat4rotateX,\n mat4rotateY,\n mat4rotateZ,\n mat4translate,\n mat4multiply\n} from '@rreusser/gl-mat4'","pinned":false,"hidden":false,"mode":"ojs"},{"id":7138,"value":"function createOrUpdateBuffer(existing, ctor, data) {\n let buffer;\n let dataSize = data.length;\n if (Array.isArray(data)) {\n let ptr = data;\n dataSize = 1;\n while (Array.isArray(ptr)) {\n dataSize *= ptr.length;\n ptr = ptr[0];\n }\n }\n if (!existing || dataSize > existing.capacity) {\n buffer = ctor(data);\n buffer.capacity = dataSize;\n } else {\n buffer = existing.subdata(data);\n buffer.capacity = dataSize;\n }\n buffer.count = dataSize;\n return buffer;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":7735,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":8318,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":5937,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":7657,"value":"function computeFisheye (out, x, amount) {\n let r = Math.pow(x[0] * x[0] + x[1] * x[1], 0.5 * amount)\n \n out[0] = x[0] / (1 + r)\n out[1] = x[1] / (1 + r)\n out[2] = x[2]\n \n return out\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":8105,"value":"import {\n createDrawBorderedLineCommand,\n createDrawLineCommand\n} from '@rreusser/quick-miterless-lines-in-webgl'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1735,"value":"import {binarySearch, decodeHash, isThumbnailDaemon, floatRgbToHex, arrowhead, hexRgbToFloat} from '@rreusser/utils'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1928,"value":"isNarrowScreen = width <= 640","pinned":false,"hidden":false,"mode":"ojs"},{"id":7794,"value":"import {checkbox} from '@jashkenas/inputs'","pinned":false,"hidden":false,"mode":"ojs"},{"id":2282,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W39","stdlib":"1"}