/*{
  "DESCRIPTION": "Newton's method made visible: every pixel is a starting guess for solving z^3 = a, coloured by which of the three roots it falls into and shaded by how fast it gets there. The target a sweeps once around the unit circle, rotating the three roots by 120 deg each, so the fractal basins spin and lock back onto themselves for a seamless loop. Jewelled magenta / cyan / gold basins divided by infinitely-fine filigree.",
  "CREDIT": "glslop agent (Claude)",
  "CATEGORIES": ["newton-fractal", "complex", "fractal", "basins", "math"],
  "INPUTS": [
    { "NAME": "speed", "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.0, "MAX": 3.0 },
    { "NAME": "scale", "TYPE": "float", "DEFAULT": 1.25, "MIN": 0.6, "MAX": 3.0 },
    { "NAME": "glow",  "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.3, "MAX": 2.5 },
    { "NAME": "spin",  "TYPE": "float", "DEFAULT": 0.15, "MIN": 0.0, "MAX": 1.0 }
  ]
}*/

#define PI  3.14159265359
#define TAU 6.28318530718
#define ITER 42
const float PERIOD = 16.0;

vec2 cmul(vec2 a, vec2 b){ return vec2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); }
vec2 cdiv(vec2 a, vec2 b){ float d = dot(b,b); return vec2(a.x*b.x + a.y*b.y, a.y*b.x - a.x*b.y)/d; }

void main(){
  vec2 uv = (2.0*gl_FragCoord.xy - RENDERSIZE.xy) / RENDERSIZE.y;

  float ph = fract(TIME * speed / PERIOD);
  float th = TAU * ph;

  // slow overall spin (must complete an integer number of turns over the loop)
  float rot = th * floor(spin * 6.0 + 0.5);
  float cr = cos(rot), sr = sin(rot);
  vec2 z = mat2(cr, -sr, sr, cr) * (uv * scale);

  // target a = e^{i th}, roots are its three cube roots (they rotate as th advances)
  vec2 a  = vec2(cos(th), sin(th));
  vec2 r0 = vec2(cos(th/3.0),             sin(th/3.0));
  vec2 r1 = vec2(cos((th + TAU)/3.0),     sin((th + TAU)/3.0));
  vec2 r2 = vec2(cos((th + 2.0*TAU)/3.0), sin((th + 2.0*TAU)/3.0));

  float it  = float(ITER);
  int   idx = 0;
  for(int i = 0; i < ITER; i++){
    vec2 z2 = cmul(z, z);
    vec2 z3 = cmul(z2, z);
    z = z - cdiv(z3 - a, 3.0 * z2);
    float d0 = length(z - r0), d1 = length(z - r1), d2 = length(z - r2);
    float dm = min(d0, min(d1, d2));
    if(dm < 0.012){
      it  = min(it, float(i) + dm / 0.012);              // smooth iteration count
      idx = d0 < d1 ? (d0 < d2 ? 0 : 2) : (d1 < d2 ? 1 : 2);
    }
  }

  float t = it / float(ITER);
  float shade = pow(1.0 - t, 1.4);                        // fast convergence = bright

  vec3 c0 = vec3(1.0, 0.22, 0.55);   // magenta
  vec3 c1 = vec3(0.20, 0.90, 1.0);   // cyan
  vec3 c2 = vec3(1.0,  0.78, 0.26);  // gold
  vec3 rc = idx == 0 ? c0 : (idx == 1 ? c1 : c2);

  vec3 col = rc * (0.12 + 0.95 * shade);

  // iteration-count contour bands -> topographic texture inside each basin
  float band = 0.5 + 0.5 * cos(it * 0.85 - 1.0);
  col *= 0.72 + 0.40 * band;
  // a brighter ridge on every few bands
  col += rc * pow(band, 8.0) * 0.5 * shade;

  // bright filigree where convergence is slow (the fractal boundary)
  float edge = smoothstep(0.28, 0.62, t) * smoothstep(1.0, 0.80, t);
  col += vec3(1.0, 0.96, 0.9) * edge * 1.3 * glow;

  col *= 1.0 - 0.18 * dot(uv, uv);                        // vignette
  col = col / (1.0 + col);
  col = pow(col, vec3(0.85));
  gl_FragColor = vec4(col, 1.0);
}