<!DOCTYPE html>
<html lang="en-us">
<head>
  <meta charset="utf-8">
  <meta name="robots" content="noindex">
  <meta name="viewport" content="width=device-width, initial-scale=1, viewport-fit=cover">
  <title>Newton vs Mobius</title>
  <style>
    html,
    body,
    canvas {
      margin: 0;
      padding: 0;
      width: 100%;
      height: 100%;
      overflow: hidden;
    }
  </style>
</head>
<body>
  <canvas id="canvas"></canvas>
  
  <script id="vertexShader" type="x-shader/x-vertex">
    #version 300 es
    in vec4 vertexPosition;
    void main() { // no-op vertex shader
      gl_Position = vertexPosition;
    }
  </script>
  
  <script id="fragmentShader" type="x-shader/x-fragment">
    #version 300 es
    precision highp float;
  
    uniform vec2 canvasSize;
    uniform float time;
    out vec4 fragColor;
  
    vec2 cmul(vec2 a, vec2 b) { // complex multiplication
      return vec2(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
    }
    vec2 cdiv(vec2 a, vec2 b) { // complex division
      float denom = (b.x*b.x+b.y*b.y);
      if (denom < 0.0000000001) denom = 0.0000000001; // avoid division by zero
      return vec2(a.x*b.x+a.y*b.y, -a.x*b.y+a.y*b.x) / denom;
    }
  
    vec2 fn(vec2 z) { // f(z) = z^3 - 1
      return cmul(z,cmul(z,z)) - vec2(1,0);
    }
    vec2 dfn(vec2 z) { // f'(z) = 3*z^2
      return cmul(vec2(3,0),cmul(z,z));
    }
  
    vec2 mobius(vec2 a, vec2 b, vec2 c, vec2 d, vec2 z) { // f(z) = (az + b)/(cz + d)
      return cdiv(cmul(a,z) + b, cmul(c,z) + d);
    }
  
    vec3 hsl2rgb(float h, float s, float l) {
      float hp = 6. * mod(h,1.);
      float c = s - s * abs(2.*l - 1.);
      float x = c - c * abs(mod(hp,2.) - 1.);
      float m = l - c/2.;
      if      (hp <= 1.) return vec3(c,x,0) + m;
      else if (hp <= 2.) return vec3(x,c,0) + m;
      else if (hp <= 3.) return vec3(0,c,x) + m;
      else if (hp <= 4.) return vec3(0,x,c) + m;
      else if (hp <= 5.) return vec3(x,0,c) + m;
      else if (hp <= 6.) return vec3(c,0,x) + m;
      else               return vec3(0,0,0);
    }
  
    float gain(float x, float k) { // https://www.iquilezles.org/www/articles/functions/functions.htm
      float a = 0.5*pow(2.0*((x<0.5)?x:1.0-x), k);
      return (x<0.5)?a:1.0-a;
    }
    float ease(float wave, float k) {
      // map cos/sin wave to [0,1], apply the gain function, and map back to [-1,1]
      return gain(wave * 0.5 + 0.5, k) * 2. - 1.;
    }
  
    vec3 newton(vec2 z) { // iteratively apply z = z-f(z)/f'(z)
      vec2 prevZ = z;
      float i = 0.;
      float intensity = 0.;
  
      for (i=1.0; i<100.; i++) {
          z -= cdiv(fn(z),dfn(z));
          if (length(z-prevZ) < 0.0001) break;
  
          // http://www.fractalforums.com/programming/smooth-colouring-of-convergent-fractals/msg33392/#msg33392
          intensity += exp(-length(z) - 0.5/(length(prevZ-z)));
  
          prevZ = z;
      }
  
      float theta = atan(z.y,z.x);
      float angle = mod(theta/6.2832+1.0, 1.0);
      float hue = mod(angle + time/50.0, 1.0);
  
      return hsl2rgb(hue, 0.7, intensity/3.-0.2);
    }
  
    vec3 draw(vec2 z) {
      float speed = 0.03;
      // The transformation is very fast-moving in the middle, so squash the animation curve to keep it under control.
      // squashFactor < 1.0 is a "rush in, slow in the middle, rush out" easing function.
      float squashFactor = 0.5;
  
      return newton(mobius(
          vec2(1,0),
          vec2(0,0), // vec2(0.005*sin(time/6.),0),// adds some rotational movement
          vec2(0, 300.*ease(cos(time*speed), squashFactor)),
          vec2(1,0),
          z
      ));
    }
  
    void main() {
      vec2 coord = (2.*gl_FragCoord.xy - canvasSize)/min(canvasSize.x, canvasSize.y); // [-1,1] in smaller dimension
      fragColor = vec4(draw(coord/10.), 1);
    }
  </script>

  <script>
    const gl = canvas.getContext("webgl2");
    if (!gl) {
      document.body.innerHTML = '<h2>Error: WebGL2 is <a href="https://get.webgl.org/webgl2/">not supported by your browser</a></h2>';
      throw "WebGL2 not supported";
    }
  
    function createShader(shaderType, sourceCode) {
      const shader = gl.createShader(shaderType);
      gl.shaderSource(shader, sourceCode);
      gl.compileShader(shader);
      if (!gl.getShaderParameter(shader, gl.COMPILE_STATUS)) throw gl.getShaderInfoLog(shader);
      return shader;
    }
  
    const program = gl.createProgram();
    gl.attachShader(program, createShader(gl.VERTEX_SHADER, vertexShader.textContent.trim()));
    gl.attachShader(program, createShader(gl.FRAGMENT_SHADER, fragmentShader.textContent.trim()));
    gl.linkProgram(program);
    if (!gl.getProgramParameter(program, gl.LINK_STATUS)) throw gl.getProgramInfoLog(program);
    gl.useProgram(program);
  
    const vertices = [[-1, -1], [1, -1], [-1, 1], [1, 1]];
    gl.bindBuffer(gl.ARRAY_BUFFER, gl.createBuffer());
    gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(vertices.flat()), gl.STATIC_DRAW);
  
    const vertexPosition = gl.getAttribLocation(program, "vertexPosition");
    gl.enableVertexAttribArray(vertexPosition);
    gl.vertexAttribPointer(vertexPosition, 2, gl.FLOAT, false, 0, 0);
  
    const canvasSizeUniform = gl.getUniformLocation(program, 'canvasSize');
    const timeUniform = gl.getUniformLocation(program, 'time');
  
    function draw() {
      const width = canvas.clientWidth;
      const height = canvas.clientHeight;
      canvas.width = width;
      canvas.height = height;
      gl.viewport(0, 0, width, height);
      gl.uniform2f(canvasSizeUniform, width, height);
  
      gl.uniform1f(timeUniform, performance.now() / 1000);
  
      gl.drawArrays(gl.TRIANGLE_STRIP, 0, vertices.length);
      requestAnimationFrame(draw);
    }
    draw();
  </script>
</body>
<!--
  © Adam Murray 2022
  https://adammurray.blog/

  Creative Commons License
  Attribution-NonCommercial-ShareAlike 4.0 International
  https://creativecommons.org/licenses/by-nc-sa/4.0/
-->