#define MAX_STEPS 100
#define MAX_DIST 100.
#define SURF_DIST .001
#define TAU 6.283185
#define PI 3.141592
#define S smoothstep
#define T iTime

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

float sdBox(vec3 p, vec3 s) {
    p = abs(p)-s;
	return length(max(p, 0.))+min(max(p.x, max(p.y, p.z)), 0.);
}


float GetDist(vec3 p, inout vec3 orbit) {
    vec3 v = vec3(cos(iTime*.12)*sin(iTime*0.73)*4. + 5., 5.+(sin(iTime*1.23)+0.5)*5., 3.+(cos(iTime*0.24)+0.5)*5.);
    vec3 n = normalize(v);
    float scale = 2.;
    p *= scale;
    
    for(int i = 0; i < 10; i++){
    
        p = abs(p);
        float lp = length(p);
        if(orbit.x > lp) orbit.x = lp;
        if(orbit.y < lp/scale) orbit.y = lp/scale;
        
        
        p -= 2.*n*max(0., dot(n, p-n*0.6));
        scale *= 2.;
        p *= 2.;
        p += n;
        
    }
    
    float d = sdBox(p, vec3(10, 10, 10))/scale;
    
    return d;
}

float GetDist(vec3 p){
    vec3 fake;
    return GetDist(p, fake);

}

float RayMarch(vec3 ro, vec3 rd, inout int i, inout float mindist, inout vec3 orbit) {
	float dO=0.;
    mindist = length(ro);
    orbit += mindist;
    for(i=0; i<MAX_STEPS; i++) {
    	vec3 p = ro + rd*dO;
        float dS = GetDist(p, orbit);
        orbit.z = fract(orbit.z + length(p));
        if(dS < mindist) mindist = dS;
        dO += dS;
        if(dO>MAX_DIST || abs(dS)<SURF_DIST) break;
    }
    orbit = fract(orbit);
    return dO;
}

vec3 GetNormal(vec3 p) {
    vec2 e = vec2(.001, 0);
    vec3 n = GetDist(p) - 
        vec3(GetDist(p-e.xyy), GetDist(p-e.yxy),GetDist(p-e.yyx));
    
    return normalize(n);
}

vec3 GetRayDir(vec2 uv, vec3 p, vec3 l, float z) {
    vec3 
        f = normalize(l-p),
        r = normalize(cross(vec3(0,1,0), f)),
        u = cross(f,r),
        c = f*z,
        i = c + uv.x*r + uv.y*u;
    return normalize(i);
}

vec3 square_bezier(vec3 a, vec3 b, vec3 c, float t){
    return mix(mix(a, b, t), mix(a, b, t), t);

}

vec3 BG(vec3 rd){
    vec3 a = vec3(0.1, 0.6, 0.2);
    vec3 b = vec3(0.2, 0.1, 0.8);
    return mix(a, b, dot(rd, vec3(0, 1, 0)) * 0.5 + 0.5);

}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (fragCoord-.5*iResolution.xy)/iResolution.y;
	vec2 m = (iMouse.xy-.5*iResolution.xy)/iResolution.xy;

    vec3 ro = vec3(cos(iTime)*3., 1., 3.*sin(iTime));
    ro.yz *= Rot(-m.y*PI+1.);
    ro.xz *= Rot(-m.x*TAU);
    
    vec3 rd = GetRayDir(uv, ro, vec3(0,0.,0), 1.);
    vec3 col = vec3(0);
    int i = 0;
    float mindist;
    vec3 orbit_trap;
    float d = RayMarch(ro, rd, i, mindist, orbit_trap);
    float iter = float(i)/float(MAX_STEPS);
    vec3 a = vec3(0.7, 0.1, 0.1);
    vec3 b = vec3(0.1, 0.1, 0.7);
    vec3 c = vec3(0.1, 0.5, 0.1);
    col += BG(rd);
    col *= S(SURF_DIST*10., SURF_DIST*20., mindist);
    
    if(d<MAX_DIST) {
        vec3 p = ro + rd * d;
        vec3 n = GetNormal(p);
        vec3 r = reflect(rd, n);

        float dif = dot(n, normalize(vec3(1,2,3)))*.5+.5;
        col = vec3(1)*pow(1.-iter, 5.)*(sin(orbit_trap*TAU+iTime)*0.3+0.5);
    }
    
    col = pow(col, vec3(.4545));	// gamma correction
    
    fragColor = vec4(col,1.0);
}