float pi = 4.*atan(1.);

mat2 rotate(float a) {
	float s = sin(a);
	float c = cos(a);
	return mat2(c, -s, s, c);
}

vec2 rand(vec2 n) { 
	float a = fract(sin(dot(n, vec2(12.9898, 4.1414))) * 43758.5453);
    float b = fract(cos(dot(n, vec2(12.9898, 4.1414))) * 43758.5453);
    return vec2(a,b);
}

float rand1(float n){return fract(sin(n) * 43758.5453123);}

vec3 mond(vec2 n) {
    float theta = rand1(atan(n.y, n.x));
    vec3 c1 = vec3(250., 201., 1.)/255.;
    vec3 c2 = vec3(1.);
    vec3 c3 = vec3(34., 80., 149.)/255.;
    vec3 c4 = vec3(1.);
    vec3 c5 = vec3(221., 1., 0.)/255.;
    float c = theta;
    if(c < 1./5.) { return c1; }
    if(c < 2./5.) { return c2; }
    if(c < 3./5.) { return c3; }
    if(c < 4./5.) { return c4; }
    return c5;
}

vec2 mult(in vec2 a, in vec2 b) {
    return vec2(a.x * b.x - a.y * b.y, a.y * b.x + a.x * b.y);
}

vec2 func(in vec2 z) {
    return mult(z - vec2(2.0*(1.0 + 0.5*sin(iTime/3.)), 0.0)*rotate(iTime/4.), z + vec2(2.0f*(1.0 + 0.5*sin(iTime/3.)), 0.0)*rotate(iTime/4.));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Normalized pixel coordinates (from 0 to 1)
    vec2 uv =  5.* ( 2.*fragCoord - iResolution.xy ) / iResolution.y;
    vec2 pos = uv;
    
    uv = func(uv);

    float r = length(uv)/4.;
    
    if (int(floor(r)) % 2 == 0) {
        float freq = 1.5*rand(vec2(floor(r),1)).y;
        uv *= rotate(freq*iTime);
    }
    else {
        float freq = 2.*rand(vec2(floor(r),-1)).y;
        uv *= rotate(-freq*iTime);
    }
    
    float theta = atan(uv.y,uv.x)+pi;

    float n = 7.;    
    float eps = 5./iResolution.x;
    float gEps = 0.01;
    float stretch = length(func(pos+vec2(gEps,0.))-func(pos))/gEps;
    float dt = abs(mod(theta+pi/n,2.*pi/n)-pi/n)/stretch;
    float ct = smoothstep(2.*eps,eps,dt*r);
    if (r < 1.)  { dt = dt/r; }
    
    //lr = ceil(log2(r+1.))-log2(r+1.)/(ceil(log2(r+1.))-floor(log2(r+1.)));
    float dr = abs(fract(r+.5)-.5)/stretch;
    float cr = smoothstep(2.*eps,eps,dr);
    float cc = max(ct,cr);
    
    vec2 dc = rand(vec2(floor(n*theta/(2.0*pi)),floor(r)));
    
    //vec2 rv = rand(vec2(floor(n*theta/(2.0*pi)),floor(r)));
    //vec3 pal = mond(rv);
    
    float val = dot(dc, vec2(1));
    // IQ's versatile cosine palette.
    vec3 col = .5 + .47*cos(6.2831*val + vec3(0, 1, 2));
    col = mix(col, vec3(0), cc);

    // Rough gamma correction.
    fragColor = vec4(sqrt(max(col, 0.)), 1);
    
    //vec3 col = mix(pal,vec3(0.),cc);

    // Output to screen
    //fragColor = vec4(col,1.0);
}