in float time;
in vec3 pos;
layout (location = 0) out vec4 fragColor;

#define speed .05
#define iResolution vec2(1000,500)

vec2 toCarte(vec2 z) { return z.x*vec2(cos(z.y),sin(z.y)); }
vec2 toPolar(vec2 z) { return vec2(length(z),atan(z.y,z.x)); }

// All the following complex operations are defined for the polar form 
// of a complex number. So, they expect a complex number with the 
// format vec2(radius,theta) -> radius*eˆ(i*theta).
// The polar form makes the operations *,/,pow,log very light and simple
// The price to pay is that +,- become costly :/ 
// So I switch back to cartesian in those cases.
vec2 zmul(vec2 z1,vec2 z2) { return vec2(z1.x*z2.x,z1.y+z2.y); }
vec2 zdiv(vec2 z1, vec2 z2) { return vec2(z1.x/z2.x,z1.y-z2.y); }
vec2 zlog(vec2 z) { return toPolar(vec2(log(z.x),z.y)); }
vec2 zpow(vec2 z, float n) { return vec2(exp(log(z.x)*n),z.y*n); }
vec2 zpow(float n, vec2 z) { return vec2(exp(log(n)*z.x*cos(z.y)),log(n)*z.x*sin(z.y)); }
vec2 zpow(vec2 z1, vec2 z2) { return zpow(exp(1.),zmul(zlog(z1),z2)); }
vec2 zadd(vec2 z1, vec2 z2) { return toPolar(toCarte(z1) + toCarte(z2)); }
vec2 zsub(vec2 z1, vec2 z2) { return toPolar(toCarte(z1) - toCarte(z2)); }

//sinz, cosz and tanz came from -> https://www.shadertoy.com/view/Mt2GDV
vec2 zsin(vec2 z) {
   z = toCarte(z);
   float e1 = exp(z.y);
   float e2 = exp(-z.y);
   float sinh = (e1-e2)*.5;
   float cosh = (e1+e2)*.5;
   return toPolar(vec2(sin(z.x)*cosh,cos(z.x)*sinh));
}

vec2 zcos(vec2 z) {
   z = toCarte(z);
   float e1 = exp(z.y);
   float e2 = exp(-z.y);
   float sinh = (e1-e2)*.5;
   float cosh = (e1+e2)*.5;
   return toPolar(vec2(cos(z.x)*cosh,-sin(z.x)*sinh));
}

vec2 ztan(vec2 z) {
    z = toCarte(z);
    float e1 = exp(z.y);
    float e2 = exp(-z.y);
    float cosx = cos(z.x);
    float sinh = (e1 - e2)*0.5;
    float cosh = (e1 + e2)*0.5;
    return toPolar(vec2(sin(z.x)*cosx, sinh*cosh)/(cosx*cosx + sinh*sinh));
}

vec2 zeta(vec2 z) {
   vec2 sum = vec2(.0);
   for (int i=1; i<20; i++) 
       sum += toCarte(zpow(float(i),-z));
   return toPolar(sum);
}

vec2 lambert(vec2 z) {
   vec2 sum = vec2(.0);
   for (int i=1; i<15; i++)
      sum += toCarte(zdiv(zpow(z,float(i)),zsub(vec2(1.,.0),zpow(z,float(i)))));
   return toPolar(sum);
}

vec2 mandelbrot(vec2 z) {
   vec2 sum = vec2(.0);
   vec2 zc = toCarte(z);
   for (int i=1; i<11; i++) 
       sum += toCarte(zpow(toPolar(sum),2.)) + zc;
   return toPolar(sum);
}

vec2 julia(vec2 z) {
    vec2 sum = toCarte(zpow(z,2.));
    // the julia set is connected if C is in the mandelbrot set and disconnected otherwise
    // to make it interesting, C is animated on the boundary of the main bulb
    // the formula for the boundary is 0.5*eˆ(i*theta) - 0.25*eˆ(i*2*theta) and came from iq
    // https://iquilezles.org/articles/mset1bulb
    float theta = fract(speed)*2.*6.2830;
    vec2 c = toCarte(vec2(0.5,.5*theta)) - toCarte(vec2(0.25,theta)) - vec2(.25,.0);
    for (int i=0; i<7; i++) sum += toCarte(zpow(toPolar(sum),2.)) + c;
    return toPolar(sum);
}

vec2 map(vec2 uv) {
  
	float t = floor(mod(time * speed,10.0));
  
	//t = 7.;
	float s = t == 1.? 4.  : t==5.? .6: t==4.? 6. : t==5.? 2.5 : t== 9. ? 13. : 3.;  
   
    uv *= s + s*.2*cos(fract(speed)*6.2830);
    
    vec2 fz, z = toPolar(uv); 
    	
    	 // z + 1 / z - 1
	fz = t == 0. ? zdiv(zadd(z,vec2(1.0)),zsub(z,vec2(1.0,.0)) ) :
         // formula from wikipedia https://en.m.wikipedia.org/wiki/Complex_analysis
		 // fz = (zˆ2 - 1)(z + (2-i))ˆ2 / zˆ2 + (2+2i)
		 t == 1. ? zdiv(zmul(zsub(zpow(z,2.),vec2(1.,0)),zpow(zadd(z,toPolar(vec2(2.,-1.))),2.)),zadd(zpow(z,2.),toPolar(vec2(2.,-2.)))) :
		 // z^(3-i) + 1.
		 t == 2. ? zadd(zpow(z,vec2(3.,acos(-1.))),vec2(1.,.0)) :
		 // tan(z^3) / z^2
		 t == 3. ? zdiv(ztan(zpow(z,3.)),zpow(z,2.)) :
		 // tan ( sin (z) )
		 t == 4. ? ztan(zsin(z)) :
		 // sin ( 1 / z )
		 t == 5. ? zsin(zdiv(vec2(1.,.0),z)) :
		 // the usual coloring methods for the mandelbrot show the outside. 
		 // this technique allows to see the structure of the inside.
		 t == 6. ? mandelbrot(zsub(z,vec2(1.,.0))) : 
         // the julia set 
		 t == 7. ? julia(z) :
		 //https://en.m.wikipedia.org/wiki/Lambert_series
		 t == 8. ? lambert(z) :
		 // this is the Riemman Zeta Function (well, at least part of it... :P)
		 // if you can prove that all the zeros of this function are 
    	 // in the 0.5 + iy line, you will win:
		 // a) a million dollars ! (no, really...)
		 // b) eternal fame and your name will be worshiped in history books
		 // c) you will uncover the deep and misterious connection between PI and the primes
		 // https://en.m.wikipedia.org/wiki/Riemann_hypothesis
    	 // https://www.youtube.com/watch?v=rGo2hsoJSbo
         zeta(zadd(z,vec2(8.,.0)));
   
 
	return toCarte(fz);  
}

vec3 color(vec2 uv) {
    float a = atan(uv.y,uv.x);
    float r = length(uv);
    
    vec3 c = .5 * ( cos(a*vec3(2.,2.,1.) + vec3(.0,1.4,.4)) + 1. );

    return c * smoothstep(1.,0.,abs(fract(log(r)-time*.1)-.5)) // modulus lines
             * smoothstep(1.,0.,abs(fract((a*7.)/3.14+(time*.1))-.5)) // phase lines
             * smoothstep(11.,0.,log(r)) // infinity fades to black
             * smoothstep(.5,.4,abs(fract(speed)-.5)); // scene switch
}

void main() {//( out vec4 fragColor, in vec2 fragCoord ){
	 vec2 uv = pos.xy; //(fragCoord.xy - iResolution.xy *.5)/iResolution.y;
     fragColor = vec4( color(map(uv)), 1.0 );
}