2e:null 2f:["$","title",null,{"children":"A Julia set on the Riemann sphere | Mark McClure | Observable"}] 34:T55e,riemann_sphere = { yield html` `; x3dom.reload(); }35:T81c,standard_pic = { let xmin = -3.5; let xmax = 3.5; let ymin = -2; let ymax = 2; // An interesting zoom // xmin = -0.1; // xmax = 0.3; // ymin = 0.9; // ymax = 1.2; let aspect = (ymax - ymin) / (xmax - xmin); let w = 800; let h = w * aspect; let bail = 100; let x_scale = d3 .scaleLinear() .domain([xmin, xmax]) .range([0, w]); let y_scale = d3 .scaleLinear() .domain([ymin, ymax]) .range([h, 0]); let canvas = d3 .create('canvas') .attr('width', w) .attr('height', h); let context = canvas.node().getContext("2d"); let canvasData = context.createImageData(w, h); for (let i = 0; i < w; i++) { for (let j = 0; j < h; j++) { let x = x_scale.invert(i); let y = y_scale.invert(j); let info = P_orbitInfo([x, y], bail); let idx = (i + j * w) * 4; let c = color(info, bail); canvasData.data[idx + 0] = c[0]; canvasData.data[idx + 1] = c[1]; canvasData.data[idx + 2] = c[2]; canvasData.data[idx + 3] = 255; } } context.putImageData(canvasData, 0, 0); // Place the two, finite fixed points context.fillStyle = 'red'; context.beginPath(); context.arc(x_scale(1), y_scale(0), 5, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(0), y_scale(0), 5, 0, 2 * Math.PI); context.fill(); // Place the critical points context.fillStyle = '#00ff00'; context.beginPath(); context.arc(x_scale(1), y_scale(0), 3, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(0.477), y_scale(0), 3, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(0.154), y_scale(0.59), 3, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(0.154), y_scale(-0.59), 3, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(0.107), y_scale(1.67), 3, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(0.107), y_scale(-1.67), 3, 0, 2 * Math.PI); context.fill(); return canvas.node(); }36:T4e0,conjugate_pic = { let xmin = -1; let xmax = 1.6; let ymin = -1.2; let ymax = 1.2; let aspect = (ymax - ymin) / (xmax - xmin); let w = 500; let h = w * aspect; let bail = 50; let x_scale = d3 .scaleLinear() .domain([xmin, xmax]) .range([0, w]); let y_scale = d3 .scaleLinear() .domain([ymin, ymax]) .range([h, 0]); let canvas = d3 .create('canvas') .attr('width', w) .attr('height', h); let context = canvas.node().getContext("2d"); let canvasData = context.createImageData(w, h); for (let i = 0; i < w; i++) { for (let j = 0; j < h; j++) { let x = x_scale.invert(i); let y = y_scale.invert(j); let info = p_orbitInfo([x, y], bail); let idx = (i + j * w) * 4; let c = color(info, bail); canvasData.data[idx + 0] = c[0]; canvasData.data[idx + 1] = c[1]; canvasData.data[idx + 2] = c[2]; canvasData.data[idx + 3] = 255; } } context.putImageData(canvasData, 0, 0); context.fillStyle = 'red'; context.beginPath(); context.arc(x_scale(1), y_scale(0), 4, 0, 2 * Math.PI); context.fill(); context.beginPath(); context.arc(x_scale(0), y_scale(0), 4, 0, 2 * Math.PI); context.fill(); return canvas.node(); }37:T440,// A version of the Julia set appropriate to use as a texture. texture = { let n = 400; let step = Math.PI / n; let theta_min = -Math.PI; let theta_max = Math.PI; let phi_min = step; let phi_max = Math.PI; let w = 2 * n; let h = n; let bail = 100; let canvas = d3 .create('canvas') .attr('width', w) .attr('height', h); let context = canvas.node().getContext("2d"); let canvasData = context.createImageData(w, h); let j = 0; for (let theta = theta_min; theta < theta_max; theta = theta + step) { let i = 0; for (let phi = phi_min; phi < phi_max; phi = phi + step) { let d = Math.tan(phi / 2); let x = Math.cos(theta) / d; let y = Math.sin(theta) / d; let info = P_orbitInfo([x, y], bail); let idx = (j + i * w) * 4; let c = color(info, bail); canvasData.data[idx + 0] = c[0]; canvasData.data[idx + 1] = c[1]; canvasData.data[idx + 2] = c[2]; canvasData.data[idx + 3] = 255; i = i + 1; } j = j + 1; } context.putImageData(canvasData, 0, 0); return canvas; }38:T4b2,// The conjugate p(z)= 1/P(1/z) of the original function // P(z) = P(x+iy) expressed in terms of x and y function p([x, y]) { let denominator = 4 - 8 * x + 12 * x ** 2 - 12 * x ** 3 + 8 * x ** 4 - 4 * x ** 5 + x ** 6 - 4 * y ** 2 + 4 * x * y ** 2 + 8 * x ** 2 * y ** 2 - 8 * x ** 3 * y ** 2 + 3 * x ** 4 * y ** 2 - 4 * x * y ** 4 + 3 * x ** 2 * y ** 4 + y ** 6; let numerator_re = 4 * x - 8 * x ** 2 + 14 * x ** 3 - 16 * x ** 4 + 12 * x ** 5 - 7 * x ** 6 + 2 * x ** 7 - 10 * x * y ** 2 + 16 * x ** 2 * y ** 2 + 8 * x ** 3 * y ** 2 - 13 * x ** 4 * y ** 2 + 6 * x ** 5 * y ** 2 - 4 * x * y ** 4 - 5 * x ** 2 * y ** 4 + 6 * x ** 3 * y ** 4 + y ** 6 + 2 * x * y ** 6; let numerator_im = 4 * y - 8 * x * y + 18 * x ** 2 * y - 24 * x ** 3 * y + 16 * x ** 4 * y - 8 * x ** 5 * y + 2 * x ** 6 * y - 6 * y ** 3 + 8 * x * y ** 3 + 16 * x ** 2 * y ** 3 - 16 * x ** 3 * y ** 3 + 6 * x ** 4 * y ** 3 - 8 * x * y ** 5 + 6 * x ** 2 * y ** 5 + 2 * y ** 7; return [numerator_re / denominator, numerator_im / denominator]; }39:T4b1,// The original function P(z) = P(x+iy) // expressed in terms of x and y function P([x, y]) { let denominator = 4 - 12 * x + 17 * x ** 2 - 20 * x ** 3 + 16 * x ** 4 - 8 * x ** 5 + 4 * x ** 6 + y ** 2 + 12 * x * y ** 2 + 8 * x ** 2 * y ** 2 - 16 * x ** 3 * y ** 2 + 12 * x ** 4 * y ** 2 - 8 * y ** 4 - 8 * x * y ** 4 + 12 * x ** 2 * y ** 4 + 4 * y ** 6; let numerator_re = 2 * x - 7 * x ** 2 + 12 * x ** 3 - 16 * x ** 4 + 14 * x ** 5 - 8 * x ** 6 + 4 * x ** 7 + y ** 2 - 4 * x * y ** 2 + 16 * x ** 2 * y ** 2 + 4 * x ** 3 * y ** 2 - 16 * x ** 4 * y ** 2 + 12 * x ** 5 * y ** 2 - 10 * x * y ** 4 - 8 * x ** 2 * y ** 4 + 12 * x ** 3 * y ** 4 + 4 * x * y ** 6; let numerator_im = 2 * y - 8 * x * y + 16 * x ** 2 * y - 24 * x ** 3 * y + 18 * x ** 4 * y - 8 * x ** 5 * y + 4 * x ** 6 * y + 8 * x * y ** 3 + 12 * x ** 2 * y ** 3 - 16 * x ** 3 * y ** 3 + 12 * x ** 4 * y ** 3 - 6 * y ** 5 - 8 * x * y ** 5 + 12 * x ** 2 * y ** 5 + 4 * y ** 7; return [numerator_re / denominator, numerator_im / denominator]; }3a:[] 32:{"id":"0f791990b240acc4","slug":"a-julia-set-on-the-riemann-sphere","title":"A Julia set on the Riemann sphere","thumbnail":"ddcda0e00945bb8b6aa330a678c80f0897af19f65f148a87eee5b4dddb147ed1","access_level":"public","create_time":"$D2020-12-17T09:48:07.210Z","update_time":"$D2020-12-17T17:02:08.654Z","view_time":null,"first_public_time":"$D2020-12-17T17:02:24.352Z","has_unpublished_edits":null,"fork":false,"stars":1,"starred":false,"owner":{"type":"team","id":"d771c8ae812d5130","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":"starter_2024"},"version":256,"publish_version":256,"role":null,"authored":false,"authors":[{"type":"individual","id":"c3301e515dd08a41","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":null}],"members":[{"type":"individual","id":"c3301e515dd08a41","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":2,"collection":{"id":"8acc996ecd7d2157","slug":"complex-dynamics","title":"Complex Dynamics","owner":{"id":"d771c8ae812d5130","login":"mcmcclur","name":"Mark McClure"}},"type":"notebook","body":{"title":"A Julia set on the Riemann sphere","cells":[{"id":0,"value":"md`# A Julia set on the Riemann sphere\n\nHere are some images of the Julia set of\n${tex.block`P(z) = \\frac{2z^{4}-2z^{3}+2z^2-z}{2z^{3}-2z^{2}+3z-2}`}\nin response to [a question on StackExchage](https://math.stackexchange.com/questions/3944442/attracting-or-parabolic-cycles-other-than-fixed-points). I suppose the highlight is this projection onto the Riemann sphere:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":130,"value":"viewof rotation = Scrubber(rotation_values)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":128,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":190,"value":"md`As conjectured in the post, every initial seed with a stable orbit eventually lands in the immediate basin of one of the three fixed points, which are:\n\n * an attractive fixed point at zero,\n * a super-attractive fixed point at one, and\n * a neutral fixed point at infinity.\n\nThe standard way to visualize the dynamics would be to choose a rectangle in the complex plane and iterate the function from each point in that rectangle. We then shade those points according to which basin they land in. The resulting picture, with the fixed points marked in red and critical points marked in green, looks something like so:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":194,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":236,"value":"md`The black region is the basin of attraction of the point at infinity. \"Attraction\" is really the wrong word here, since point is really a *neutral* fixed point, i.e. ${tex`P'(\\infty)=1`}. A good way to understand what that means is to examine the conjugate function \n${tex.block`p(z) = 1/P(1/z) = \\frac{2 z^4-3 z^3+2 z^2-2 z}{z^3-2 z^2+2 z-2}.`}\nFrom there we can compute ${tex`p(0)=0`} and ${tex`p'(0)=1`}. The Julia set of ${tex`p`} looks like so:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":217,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":245,"value":"md`## Code`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":90,"value":"$37","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":40,"value":"function color(info, bail) {\n let v = 255 - (255 * info.cnt) / (bail + 1);\n let lim = info.lim;\n if (lim == 0) {\n return [0, 0, v];\n } else if (lim == 1) {\n return [v, v, 0];\n // } else if (lim == 2) {\n // return [255 - v, 255 - v, 255 - v];\n } else {\n return [0, 0, 0];\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":215,"value":"function p_orbitInfo([x0, y0], bail) {\n let cnt = 0;\n let x = x0;\n let y = y0;\n let lim = p_classify([x, y]);\n while (lim == 3 && cnt < bail) {\n [x, y] = p([x, y]);\n lim = p_classify([x, y]);\n cnt = cnt + 1;\n }\n return { lim: lim, cnt: cnt };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":22,"value":"function P_orbitInfo([x0, y0], bail) {\n let cnt = 0;\n let x = x0;\n let y = y0;\n let lim = P_classify([x, y]);\n while (lim == 3 && cnt < bail) {\n [x, y] = P([x, y]);\n lim = P_classify([x, y]);\n cnt = cnt + 1;\n }\n return { lim: lim, cnt: cnt };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":213,"value":"function p_classify([x, y]) {\n if ((x - 0.3) ** 2 + y ** 2 > 2) {\n return 0;\n } else if ((x - 1) ** 2 + y ** 2 < 0.01) {\n return 1;\n } else {\n return 3;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":11,"value":"function P_classify([x, y]) {\n if (x ** 2 + y ** 2 < 0.25) {\n return 0;\n } else if ((x - 1) ** 2 + y ** 2 < 0.015625) {\n return 1;\n } else if (Math.abs(y) > 1.2) {\n return 2;\n } else {\n return 3;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":211,"value":"$38","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2,"value":"$39","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":133,"value":"rotate = d3\n .select(riemann_sphere)\n .select('#rotation')\n .attr('rotation', `0 1 0 ${rotation}`)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":126,"value":"rotation_values = d3.range(0, 2 * Math.PI, Math.PI / 400)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":124,"value":"import { Scrubber } from \"@mbostock/scrubber\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":26,"value":"d3 = require('d3@6')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":121,"value":"x3dom = require('x3dom').catch(() => window['x3dom'])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":118,"value":"html`\n \n`","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W3a","stdlib":"1","titling":"auto"}} 33:"$32:body"