#version 300 es
precision highp float;
uniform float iTime;
uniform vec2  iResolution;
uniform float a0; // value=-0.14, min=-1.0, max=1.0, step=0.01
uniform float a1; // value=0.0, min=-1.0, max=1.0, step=0.01
uniform float b0; // value=0.02, min=-1.0, max=1.0, step=0.01
uniform float b1; // value=0.0, min=-1.0, max=1.0, step=0.01
uniform float c0; // value=-0.02, min=-1.0, max=1.0, step=0.01
uniform float c1; // value=0.0, min=-1.0, max=1.0, step=0.01
uniform float d0; // value=-0.08, min=-1.0, max=1.0, step=0.01
uniform float d1; // value=0.0, min=-1.0, max=1.0, step=0.01
uniform float mobius; // value=1.0, min=0.0, max=1.0, step=0.01
uniform float menger; // value=0.0, min=0.0, max=1.0, step=0.01
uniform float colour; // value=1, min=0, max=1, step=1 (No, Yes)
out vec4 fragColor;

// 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));
}
float hash(float n) {
    return fract(sin(n) * 43758.5453123);
}
float hash2(float n) {
    return fract(sin(n) * 4375345348.5453123);
}

// Raymarching code from Kishimisu
float fov = 1.1;
vec3 axis = vec3(1,0,0);

float  sdSphere(vec3 p, float s){
	return length(p) - s;
}
float max3 (vec3 v) {
  return max (max (v.x, v.y), v.z);
}
float sdBox(vec3 p, vec3 b){
	vec3 q = abs(p) - b;
	return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}
float opUnion(float d1, float d2){
	return min(d1, d2);
}
float opSubtraction(float d1, float d2){
	return max(-d1, d2);
}
float opIntersection(float d1, float d2){
	return max(d1, d2);
}
float opSmoothUnion(float d1, float d2, float k){
	float h = clamp(0.5 + 0.5 * (d2-d1)/k, 0.0, 1.0);
	return mix(d2, d1, h) - k*h*(1.0-h);
}
mat2 rot2D(float angle){
	float s = sin(angle);
	float c = cos(angle);
	return mat2(c, -s, s, c);
}
mat3 rot3D(vec3 axis, float angle){
	axis = normalize(axis);
	float s = sin(angle);
	float c = cos(angle);
	float oc = 1.0 - c;
	
	return mat3(
	oc*axis.x*axis.x+c,
	oc*axis.x*axis.y-axis.z*s,
	oc*axis.z*axis.x+axis.y*s,
	oc*axis.x*axis.y+axis.z*s,
	oc*axis.y*axis.y+c,
	oc*axis.y*axis.z-axis.x*s,
	oc*axis.z*axis.x-axis.y*s,
	oc*axis.y*axis.z+axis.x*s,
	oc*axis.z*axis.z+c
	);
}
vec3 rot3D(vec3 p, vec3 axis, float angle){
	return mix(dot(axis,p)*axis, p, cos(angle)) + cross(axis, p) * sin(angle);
}

// Menger sponge implementation from Inigo Quilezles
float sdMengerSponge(vec3 p) {

    float d = sdBox(p,vec3(1.0));
    
       float s = 1.0;
       for( int m=0; m<3; m++ )
       {
          vec3 a = mod( p*s, 2.0 )-1.0;
          s *= 3.0;
          vec3 r = abs(1.0 - 3.0*abs(a));
    
          float da = max(r.x,r.y);
          float db = max(r.y,r.z);
          float dc = max(r.z,r.x);
          float c = (min(da,min(db,dc))-1.0)/s;
    
          d = max(d,c);
       }
    
       return d;
}

float map(vec3 p){
	p = mod(p, 1.) - 0.5;
	
	float rightbox = sdBox(p, vec3(0.5,0.05 ,0.05));
	float topbox = sdBox(p, vec3(0.05,0.5 ,0.05));
	float backbox = sdBox(p, vec3(0.05,0.05 ,0.5));
	
	float sec = floor(iTime + 0.5);
	
	float r = step(0.5,  hash(sec));
	float dir = step(0.5, hash2(sec))*2. -1.;
    
    float x = mix(p.x,(mod(p.x + dir*iTime,1.)-0.5), r);
    float y = mix(p.y,(mod(p.y+ dir*iTime,1.)-0.5), 1.-r);
	
	float mainbox = sdBox(vec3(x,y,p.z), vec3(0.15,0.15 ,0.15));
	float mengerbox = sdMengerSponge(vec3(x,y,p.z) / 0.15) * 0.15;   
	

	mainbox = mix(mainbox, mengerbox, menger);
	
	
	return opSmoothUnion(min(min(rightbox, topbox),backbox),mainbox, 0.05);
}

vec3 palette( in float t, in vec3 a, in vec3 b, in vec3 c, in vec3 d )
{
    return a + b*cos( 6.28318*(c*t+iTime*0.1+d) );
}

void main()
{
	
	vec2 uv = (2.* gl_FragCoord.xy - iResolution.xy)/iResolution.y;
	
	vec2 a = vec2(a0, a1);
	vec2 b = vec2(b0, b1);
	vec2 c = vec2(c0, c1);
	vec2 d = vec2(d0, d1);
	
	vec2 numerator = cMul(a, uv) + b;
	vec2 denominator = cMul(c, uv) + d;
	vec2 transformed = cDiv(numerator, denominator);
	
	uv = mix(uv, transformed, mobius);
	
	
	// Initialisation
	vec3 ro = vec3(0, 0,mod(iTime,1.));
	vec3 rd = normalize(vec3(uv * fov, 1));
	vec3 col = vec3(0.);
	
// 	ro.xz *= rot2D(iTime*0.5);
 	//rd.xz *= rot2D(iTime*0.5);
    //ro.xyz *= rot3D(vec3(10.*cos(iTime*0.14),10.*cos(iTime*0.15),10.*cos(iTime*0.16)),iTime*0.5); 
    rd.xyz *= rot3D(vec3(10.*cos(iTime*0.13),10.*cos(iTime*0.12),10.*cos(iTime*0.11)),iTime*0.5);
	
	float t = 0.05;
	
	int i ;
	// Raymarching
	for(i = 0; i<150; i++){
		vec3 p = ro + rd*t; // position along the ray
		
		p.y += (cos(t + iTime*0.7)*0.2)*(0.5+0.5*cos(iTime*0.13));
		p.x += (sin(t + iTime*0.6)*0.2)*(0.5+0.5*cos(iTime*0.11));
		
		float d = map(p);
		
		t += d;
		
        if (d<.0001 || t>100.){
			break;
		}
	}
	
	col = vec3(t*1.);
	if(colour>0.){
	    col = palette((t*0.2 + float(i)*0.005)*(0.7 + 0.4*cos(iTime)), vec3(0.8, 0.5, 0.4), vec3(0.2, 0.4, 0.2), vec3(2.0, 1.0, 1.0), vec3(0.00, 0.25, 0.25));
	}
		
	fragColor = vec4(col,1.);
}