#version 300 es
precision highp float;

uniform vec2  iResolution;
uniform float iTime;

uniform float MaxIter; // value=100, min=1, max=1000, step=1
uniform int AA; // value=4, min=1, max=9, step=1
out vec4 fragColor;

vec2 cMul(vec2 a, vec2 b) {
    // Somehow it becomes undefined for some radians.
    return vec2( a.x*b.x -  a.y*b.y,a.x*b.y + a.y * b.x);
}

vec2 cInverse(vec2 a) {
    return  vec2(a.x,-a.y)/dot(a,a);
}

vec2 cExp(in vec2 z){
    return vec2(exp(z.x)*cos(z.y),exp(z.x)*sin(z.y));
}

vec2 cDiv(vec2 a, vec2 b) {
    return cMul(a,cInverse(b));
}

vec2 newton(vec2 z) {
	return cMul(z - cDiv(cMul(z,cMul(z,z)) - vec2(1.0,0.0), 3.0*cMul(z,z)),cExp(vec2(0.,0.5*(0.5+0.5*cos(iTime)))));	
}

vec4 render(float xin, float yin){
    vec2 z = vec2(xin, yin);

    float j = 0.;
    
    while(j<MaxIter){
        z = newton(z);
        j++;
    }
    // Don't mind this mess. 
    vec3 col = min(max(z.x*sqrt(xin*xin+yin*yin), 0.)*10.,1.)*vec3(0.941, 0.635, 0.007);
    col += min(max(z.y*sqrt(xin*xin+yin*yin), 0.)*10.,1.)*vec3(0.945, 0.533, 0.019);
    col += min(max(z.x*z.y*sqrt(xin*xin+yin*yin), 0.)*10.,1.)*vec3(0.850, 0.364, 0.223);

    return vec4(col, 1.0); 
}

void main() {
    for(int m=0; m<AA; m++){
        for(int n=0; n<AA; n++){
            vec2 z   = (2.0*(gl_FragCoord.xy+vec2(float(m),float(n))/float(AA))-iResolution.xy) / iResolution.y;
            fragColor += render(z.x, z.y);
        }
    }
    fragColor/=float(AA*AA);

}