#version 300 es
precision highp float;

uniform float iTime;
uniform vec2  iResolution;

out vec4 fragColor;

#define PI 3.14159265359
#define gamma 22.5836845297


uniform float Inverse; // value = 0.0, min=0.0, max=1.0 step=1.0 (No, Yes)
uniform float Animate; // value = 0.0, min=0.0, max=1.0 step=1.0 (No, Yes)
uniform float DeformationScale; // value=1.0, min=0.0, max=1.0, step=0.01
uniform float Droste; // value=22.58368, min=1.0, max=256.0, step=0.00001
float deformationScale;
float d16 = 1./16.;
float td16 = 2./16.;



// Complex operators from Roy Wiggins
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 cInverse(vec2 a) {
    return  vec2(a.x,-a.y)/dot(a,a);
}
vec2 cDiv(vec2 a, vec2 b) {
    return cMul( a,cInverse(b));
}
vec2 cExp(in vec2 z){
    return vec2(exp(z.x)*cos(z.y),exp(z.x)*sin(z.y));
}

vec2 cLog(vec2 a) {
    float b =  atan(a.y,a.x);
    // 0.4 adjustment rotation
    return vec2(log(length(a)),b + 0.4*deformationScale/256.);
}

vec2 escher(vec2 z) {
    z = cLog(z);
    if(Animate>0.0){
	    z.y -= iTime/2.;
    }
    // Multiply by inverse of h formula paper
    vec2 denom = cDiv(vec2(log(deformationScale), 2.0*PI), vec2(0,2.0*PI));
    z = cMul(denom, z);
    //z.x = mod(z.x-iTime,log(256.));
    z = cExp(z);
    return z;
}

// From Reinder's prentententoonstellingd16, td16, d16
vec3 droste(vec2 z ) {
    float t = Droste;
    float changed = 0.0;
	for( int i=0; i<3; i++ ) {
		if(any(greaterThan(abs(z),vec2(1.)))) {
			z *= (1./t);
            changed = 1.0;
		}		
		if(all(lessThan(abs(z),vec2(1./t)))) {
			z *= t;
            changed = 1.0;
		}
	}
	return vec3(z, changed);
}

vec2 escher_inverse(vec2 z) {
    z = cLog(z);
    if(Animate>0.0){
	    z.y -= iTime/2.;
    }
    vec2 denom = cDiv(vec2(log(deformationScale), 2.0*PI), vec2(0,2.0*PI));
    z = cDiv(z, denom);
    //z.x = mod(z.x,log(256.));
    z = cExp(z);
    return z;
}

void main() {
    deformationScale = pow(DeformationScale, 2.) * 255. + 1.;
    float[] gaussian = float[] (d16, td16, d16, td16, 4./16., td16, d16, td16, d16);
    vec3 col = vec3(0.);
	for (int x=-1; x<=1; x++) {
    	for (int y=-1; y<=1; y++) {
            vec2 z   = (2. * gl_FragCoord.xy + vec2(x, y)  - iResolution) / iResolution.y;
			
            //z = escher(z);
            if(Inverse>0.0){
                z = escher_inverse(z);
            }
            else{
                z = escher(z);
            }
            
            
            vec3 dro;
            dro = droste(z);
			if(dro.z>0.0){
			    col *= vec3(1.0,0.0,0.0);
			}
			
            z = vec2(dro.x, dro.y);
            z = .5 + .5*z;

            //mod(z,1.0);

            z = sin(z*20.0)*10.;
            col += gaussian[(x+1)*3 + (y+1)]*vec3(z.x*z.y)*0.12;
        }
    }
    	
	
	fragColor = vec4( col, 1.0);
}