// inspired by:
// https://www.shadertoy.com/view/Xd3SDs
// https://www.shadertoy.com/view/Mtc3zS
// https://www.shadertoy.com/view/Xds3zN
// https://www.shadertoy.com/view/lsBXW3
// https://www.shadertoy.com/view/MtVGDt
// https://www.shadertoy.com/view/MsfGRr

//http://140.129.20.249/~jmchen/cg/docs/rendering%20pipeline/rendering/light_specular.html
//http://ogldev.atspace.co.uk/www/tutorial19/tutorial19.html

#define pi2            (2.*3.141593)
#define rotate(plane,a)  (plane=vec2(cos((a)*pi2)*plane.x+sin((a)*pi2)*plane.y,cos((a)*pi2)*plane.y-sin((a)*pi2)*plane.x))
#define pmod(a,b)    ( mod(mod((a),(b))+(b),(b)) )
#define rep(a,r)    ( pmod(((a)+(r)*.5),(r))-(r)*.5 )
#define repxz(a,r)    vec3( rep((a).x,(r)), (a).y + 0.5*sin((a).x+iTime)*cos((a).z+iTime), rep((a).z,(r)) )

#define EPSILON .05

vec2 opU(vec2 d1, vec2 d2)
{
	return (d1.x < d2.x) ? d1 : d2;
}

vec3 opTwist(vec3 p)
{
    float  c = cos(.05*sin(iTime)*p.z);
    float  s = sin(.05*cos(iTime)*p.z);
    mat2   m = mat2(c,-s,s,c);

    return vec3(m*p.xy,p.z);
}

float udBox(vec3 p, vec3 b)
{
    return length(max(abs(p)-b, 0.));
}

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

float collumns(vec3 p, vec3 c, float r)
{
	vec3 op = p;
    float d;

    d = udBox(p, vec3(20., 6., 20.));
    d = max(d, -sdBox(p, vec3(3., 3.5, 3.)));
	p = repxz(p, r);
    d = max(d, udBox(p, vec3(.25, 3., .25)));

    vec3 cell_center = op - p;
    float cut = 1000000.;

    if( textureLod(iChannel3, (cell_center.xz/r+.5)/64., 0.).x < .8 )
    {
       cut = sdBox(op-cell_center, vec3(.5, 3.5, .5));
    }

    return max(d,-cut);
}

float surface(vec3 p, vec3 s, float r)
{
    float f = udBox(p, vec3(20., 1., 20.));
    p = repxz(p, r);
    f = max(f, udBox(vec3(p.x, p.y, p.z), s));

    return f;
}

vec4 cloud(vec3 pos)
{
    rotate(pos.xy,iTime/50.);
    vec3 p = pos;
    p.xy += iTime * 5.;
    float noise = textureLod(iChannel3,(p.xy)/64., 0.).x;
    vec4 col = vec4(pos.z*.5+.5);
    col.a = noise * .25;

    return col;
}

float julia(vec3 p, in vec4 c)
{
	vec4 z = vec4(p, 0.);
    float md2 = 1.;

    for( int i=0; i<11; i++ )
    {
        // https://www.cs.cmu.edu/~kmcrane/Projects/QuaternionJulia/paper.pdf
        // http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/435/
        // |dz|^2 -> 4*|dz|^2
        // z'n+1 = 2zn * z'n
        md2 = 4.*dot(z, z)*md2;

        // z -> z2 + c
        z = vec4( z.x * z.x - z.y * z.y - z.z * z.z - z.w * z.w,
                  2.0 * z.x * z.y,
                  2.0 * z.x * z.z,
                  2.0 * z.x * z.w) + c;

        // Just like in 2D Julia set: https://www.shadertoy.com/view/XlSyDK
        if(dot(z, z)>4.) break;
    }

	// iq explained the following formula here:
	// http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/msg8505/#msg8505
    // http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/450/
    // https://en.wikipedia.org/wiki/Koebe_quarter_theorem
    // DE = (power/4) · |Z|·log |Z| / |dZ|
    return .25*sqrt(dot(z, z)/md2)*log(dot(z, z));
}

vec2 fn(vec3 p, in vec4 juliaFactor)
{
    vec2 f = vec2(1.);

    f = opU(f, vec2(julia(p, juliaFactor), 0.));

    p = opTwist(p);

    f = opU(f, vec2(surface(vec3(p.x, p.y + 3.3, p.z), vec3(.54, .25, .54), 1.15), .9));
    f = opU(f, vec2(surface(vec3(p.x, p.y - 2.3, p.z), vec3(.54, .25, .54), 1.15), .6));
    f = opU(f, vec2(collumns(p, vec3(0.), 1.15), .15));

    return f;
}

vec3 getNormal(vec3 pos, vec4 juliaFactor)
{
    vec2 e = vec2(EPSILON, 0.);    // delta (epsilon)

    vec3 n = normalize(vec3(
            fn(vec3(pos.x+e.x, pos.y, pos.z), juliaFactor).x - fn(vec3(pos.x-e.x, pos.y, pos.z), juliaFactor).x,
            fn(vec3(pos.x, pos.y+e.x, pos.z), juliaFactor).x - fn(vec3(pos.x, pos.y-e.x, pos.z), juliaFactor).x,
            fn(vec3(pos.x, pos.y, pos.z+e.x), juliaFactor).x - fn(vec3(pos.x, pos.y, pos.z-e.x), juliaFactor).x));

    return n;
}

vec3 getTexture(vec3 n, sampler2D tex)
{
	vec3 mat = vec3(0.);

	// material color from texture
	float u = atan(n.z, n.x) / 3.1415*2.;
	float v = asin(n.y) / 3.1415*2. + .5;
	mat = texture( tex, vec2(u,v)).xyz;

	return mat;
}

// https://iquilezles.org/articles/rmshadows
float softShadow( vec3 rayOrigin, vec3 rayDir, float tMin, float tMax, in vec4 juliaFactor)
{
	float res = 1.;
	float k = 5.;
    for( float t=tMin; t < tMax; )
    {
		float h = fn( rayOrigin + rayDir*t, juliaFactor).x;
        if( h<.001 || t>tMax)
            return 0.;
		res = min( res, k*h/t );
        t += clamp( h, .05, .1 );
    }
    return clamp( res, 0., 1. );
}

vec3 calcLighting(vec3 rayDir, vec3 rayPos, vec2 dist, in vec4 juliaFactor)
{
    vec3 material;
    vec3 texture[3];
    vec3 color = vec3(0.);
    vec3 factor = vec3(0.);

    vec3 normal = getNormal(rayPos, juliaFactor);
    vec3 sunCol = vec3(1., 1., 1.);
    vec3 ref = reflect(rayDir, normal);

    texture[0] = getTexture(normal, iChannel0);
    texture[1] = getTexture(ref, iChannel1);
    texture[2] = getTexture(ref, iChannel2);

    vec3 light = vec3(sin(iTime / 4.), -1., 1);

	// https://www.shadertoy.com/view/Xds3zN

    if(dist.y < 0.1f)
        material = 0.1 + 0.1*(sin( vec3(0.05*sin(rayPos.x),0.08*tan(rayPos.y),0.1*cos(rayPos.z))*100.0));
    else
    {
        float tmpCol = 1.f + sin(iTime) / 2. + dist.y;
        material = 0.1 + 0.1*(sin( vec3(0.05,0.08,0.1)*(100.0*tmpCol)));
    }

    if(dist.y < 0.1)
        color = clamp(material + .75*texture[0], 0., 1.);
    else if (dist.y > .1 && dist.y < .2)
        color = clamp(.5*material + .5*texture[1], 0., 1.);
   	else
    	color = clamp(.5*material + .5*texture[2], 0., 1.);

	float ambient = clamp(0.5+0.5*normal.y, 0., 1.);
    float diffuse = clamp(dot(normal, light), 0., 1.);
    float backlight = clamp(dot(normal, normalize(vec3(-light.x, 0.,-light.z))), 0., 1. )*clamp( 1.-rayPos.y, 0., 1.);
    float skydome = smoothstep( -.1, .1, ref.y);
    float fresnel = pow(clamp(1.+dot(normal, rayDir), 0., 1.), 2.);
	float specular = pow(clamp(dot(ref, light), 0., 2.), 16.);

    diffuse *= softShadow(rayPos, light, .02, 2.5, juliaFactor);
    skydome *= softShadow(rayPos, ref, .02, 2.5, juliaFactor);

    factor += 2.5*diffuse*vec3(.5, .5, .5);
	factor += 0.1*specular*vec3(.5, .5, .5)*diffuse;
    factor += 0.5*ambient*vec3(.5, .5, .5);
    factor += 0.5*skydome*vec3(.5, .5, .5);
    factor += 0.5*backlight*vec3(.5, .5, .5);
    factor += 1.0*fresnel*vec3(1., 1., 1.)*(5. + sin(iTime * .5)*5.);
	color *= factor;

    color = mix( color, vec3(1.,1.,1.), 1.-exp( -.0002*dist.x*dist.x*dist.x ));

	// https://iquilezles.org/articles/outdoorslighting
	// gamma correction
	color = pow( color, vec3(1./2.2));

    return color;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec4 juliaFactor = .5*cos(iTime*vec4(.25, .1, .6, .4));
    vec2 screenPos = -1. + 2. * fragCoord.xy / iResolution.xy; // screenPos can range from -1 to 1
	screenPos.x *= iResolution.x / iResolution.y; // Correct aspect ratio

    float angle = iTime * .1;

    vec3 front = normalize(vec3(sin(angle),0,cos(angle)));
    vec3 up = vec3(0., -1., 0.);
    vec3 right = cross(up, front);
    vec3 pos = vec3(0., 0., 0.)-front*2.5;

    vec3 rayDir = normalize(front + screenPos.x*right + screenPos.y*up);

    vec3 rayPos = pos;

    vec2 dist;

    float step = .2;
    vec4 color = vec4(0., 0., 0., 1.);

    rayPos += rayDir*abs(sin(screenPos.x*200.*cos(screenPos.y*200.+iTime)))*step;

    for(int i=0;i<100;i++)
    {
        dist = fn(rayPos, juliaFactor);
        rayPos += dist.x*rayDir;

        vec3 p = rayPos;

        p += vec3(0., .5, 0.);
        rotate(p.yz, sin(70.));
        vec4 c = cloud(p);

        color.r += color.a*c.a*sin(rayPos.x);
        color.g += color.a*c.a*cos(rayPos.y);
        color.b += color.a*c.a*sin(rayPos.z);
        color.a = 1.;

        if(dist.x<.001) break;
    }

    if (dist.x < EPSILON)
        fragColor = mix(color, vec4(calcLighting(rayDir, rayPos, dist, juliaFactor), 0.), .9);
    else
        fragColor = vec4(vec3(.85, .95, 1.)-pow(screenPos.y, 2.)/1.5, 0.);
}