#define STEPS 30
#define FAR 60.0
#define PIXELR 0.5/iResolution.x

#define BOUNCES 3
#define SAMPLES 4.0

float CTIME = 0.0;

//Hash methods from https://www.shadertoy.com/view/4djSRW
//#define HASHSCALE3 vec3(.1031, .1030, .0973)
#define HASHSCALE3 vec3(443.897, 441.423, 437.195)
vec3 hash33(vec3 p3){
	p3 = fract(p3 * HASHSCALE3);
    p3 += dot(p3, p3.yxz+19.19);
    return fract((p3.xxy + p3.yxx)*p3.zyx);

}

#define HASHSCALE1 443.8975
float hash13(vec3 p3){
	p3  = fract(p3 * HASHSCALE1);
    p3 += dot(p3, p3.yzx + 19.19);
    return fract((p3.x + p3.y) * p3.z);
}


//Distance functions from Mercury's SDF library
//http://mercury.sexy/hg_sdf/

// Maximum/minumum elements of a vector
float vmax3(vec3 v) {
	return max(max(v.x, v.y), v.z);
}

// Box: correct distance to corners
float fBox(vec3 p, vec3 b) {
	vec3 d = abs(p) - b;
	return length(max(d, vec3(0))) + vmax3(min(d, vec3(0)));
}

// Rotate around a coordinate axis (i.e. in a plane perpendicular to that axis) by angle <a>.
// Read like this: R(p.xz, a) rotates "x towards z".
// This is fast if <a> is a compile-time constant and slower (but still practical) if not.
void pR(inout vec2 p, float a) {
	p = cos(a)*p + sin(a)*vec2(p.y, -p.x);
}


// 3D noise function (IQ)
float noise(vec3 p){
	vec3 ip = floor(p);
    p -= ip;
    vec3 s = vec3(7.0,157.0,113.0);
    vec4 h = vec4(0.0, s.yz, s.y+s.z)+dot(ip, s);
    p = p*p*(3.0-2.0*p);
    h = mix(fract(sin(h)*43758.5), fract(sin(h+s.x)*43758.5), p.x);
    h.xy = mix(h.xz, h.yw, p.y);
    return mix(h.x, h.y, p.z);
}

vec2 dist(vec3 p){
    vec3 pp = p;
    pR(pp.xy, CTIME);
    float tunnel = -fBox(pp, vec3(4.0, 4.0, 2.0*FAR));
    
    pp = p;
    pR(pp.xz, CTIME*0.75);
    pR(pp.yz, CTIME*0.25);
    float box = fBox(pp, vec3(0.5))-noise((pp+CTIME*0.5)*0.5)*2.0;
    
    float scene = min(tunnel, box);
    float id = 0.0;
    
    if(box < tunnel){
        id = 1.0;
    }
    
    return vec2(scene, id);
}

vec3 normals(vec3 p){
    vec3 eps = vec3(PIXELR, 0.0, 0.0);
    return normalize(vec3(
        dist(p+eps.xyy).x-dist(p-eps.xyy).x,
        dist(p+eps.yxy).x-dist(p-eps.yxy).x,
        dist(p+eps.yyx).x-dist(p-eps.yyx).x
    ));
}

//Enhanced sphere tracing algorithm introduced by Mercury

// Sign function that doesn't return 0
float sgn(float x) {
	return (x < 0.0)?-1.0:1.0;
}

vec2 march(vec3 ro, vec3 rd){
    float t = 0.001;//EPSILON;
    float step = 0.0;

    float omega = 1.0;//muista testata eri arvoilla! [1,2]
    float prev_radius = 0.0;

    float candidate_t = t;
    float candidate_error = 1000.0;
    float sg = sgn(dist(ro).x);

    vec3 p = vec3(0.0);

	for(int i = 0; i < STEPS; ++i){
		p = rd*t+ro;
		float sg_radius = sg*dist(p).x;
		float radius = abs(sg_radius);
		step = sg_radius;
		bool fail = omega > 1. && (radius+prev_radius) < step;
		if(fail){
			step -= omega * step;
			omega = 1.;
		}
		else{
			step = sg_radius*omega;
		}
		prev_radius = radius;
		float error = radius/t;

		if(!fail && error < candidate_error){
			candidate_t = t;
			candidate_error = error;
		}

		if(!fail && error < PIXELR || t > FAR){
			break;
		}
		t += step;
	}
    //discontinuity reduction
    float er = candidate_error;
    for(int j = 0; j < 6; ++j){
        float radius = abs(sg*dist(p).x);
        p += rd*(radius-er);
        t = length(p-ro);
        er = radius/t;

        if(er < candidate_error){
            candidate_t = t;
            candidate_error = er;
        }
    }
	if(t <= FAR || candidate_error <= PIXELR){
		t = candidate_t;
	}
    
    p = ro+rd*t;
    float id = dist(p).y;
    
	return vec2(t, id);
}

//returns material of the object hit
// emissive color is xyz, and reflectance w
vec4 getMaterial(float obj, vec3 p){
    vec3 base = vec3(0.0);
    float reflectance = 0.0;
    float m = mod(p.z-(CTIME*10.0), 8.0) - 4.0;
    vec3 col = vec3(0.4, 0.3, 0.8);
    
    if(obj == 0.0){
        if(m > 0.0 && m > 2.0){
            base = col;
        }
        else if( m < 0.0 && m > -2.0){
            base = col.bgb;
        }
        reflectance = m > 0.0 ? 0.2 : 0.5;
    }
    else if(obj == 1.0){
        base = col.brg;
        reflectance = 0.8;
    }
    
    
    return vec4(base, reflectance);
}

vec3 render(vec3 o, vec3 d, vec2 uv){
    
    vec3 ro = o;
    vec3 rd = d;
    
    vec3 pixel_color = vec3(0.0);
    vec3 absorption_factor = vec3(1.0);
    
    for(int i = 0; i < BOUNCES; ++i){
        vec2 t = march(ro, rd);
        vec3 p = ro+rd*t.x;
        
        if(t.y < 0.0 || t.x > FAR){
            break;
        }
        
        //material.xyz == emissive
        //material.w == reflectance
        vec4 material = getMaterial(t.y, p);
        pixel_color += material.xyz * absorption_factor;
        absorption_factor *= material.w;
        
        vec3 n = normals(p);
        ro = p+(n*0.02);
        if(t.y == 0.0){
            rd = reflect(rd,n);
            //Thanks to fizzer to introducing this skew thing! :)
            rd = normalize(rd + (hash33(vec3(uv, float(i))) - 0.5)*0.1); 
        }
        else if(t.y == 1.0){
            rd = reflect(rd,n);
        }
        
    }
    
    return pixel_color;
}


void mainImage( out vec4 fragColor, in vec2 fragCoord ){
    vec2 uv = fragCoord.xy / iResolution.xy;
    vec2 q = -1.0+2.0*uv;
    q.x *= iResolution.x/iResolution.y;
    
    vec3 ro = vec3(0.0, 0.0, 4.0);
    vec3 rt = vec3(0.0, 0.0, -2.0);

    vec3 z = normalize(rt-ro);
    vec3 x = normalize(cross(z, vec3(0.0, 1.0, 0.0)));
    vec3 y = normalize(cross(x, z));
    
    vec3 color = vec3(0.0);
    
    for(float i = 0.0; i < SAMPLES; ++i){
        //from https://iquilezles.org/articles/simplepathtracing
        CTIME = (iTime-iTimeDelta) + 0.6*(1.0/24.0)*hash13(vec3(uv, iTime*0.01));
        
    	vec3 rd = normalize(mat3(x, y, z)*vec3(q, radians(60.0)));
    	color += render(ro, rd, uv);
    }
    color /= SAMPLES;
    color = smoothstep(0.2, 0.9, color);
    
    color = pow(color, 1.0/vec3(2.2));

	fragColor = vec4(color, 1.0);
}