#define AO_DIST 0.085
#define AO_K 2.0
#define DISTORTION 0.2
#define GLOW 6.0
#define AMBIENT 0.0
#define SCALE 3.0
#define PERLIN_GRID 8.0

const float k = 5.0;

const float PI = 3.14159265359;
const float RAD = 3.14159265359 / 180.0;

const vec3 cameraPos = vec3(0.0, 0.0, -20.0);

const vec3 lightPos = vec3(0.0, 28.0, 32.0);
const vec3 fakeLightPos = vec3(0.0, 15.0, 32.0);

vec3 rotateY(vec3 p, float a) {
    return vec3(cos(a) * p.x + sin(a) * p.z, p.y, -sin(a) * p.x + cos(a) * p.z);
}

vec3 rotateZ(vec3 p, float a) {
    return vec3(cos(a) * p.x + sin(a) * p.y, -sin(a) * p.x + cos(a) * p.y, p.z);
}

vec2 noise2D2D( vec2 p ) {
    return fract(sin(vec2(dot(p, vec2(230.2, 409.1)),
                 dot(p, vec2(909.1, 891.5))))
                 * 58290.101);
}

float surflet(vec2 p, vec2 gridPoint) {
    // Compute the distance between p and the grid point along each axis, and warp it with a
    // quintic function so we can smooth our cells
    vec2 t2 = abs(p - gridPoint);
    vec2 t = vec2(1.f) - 6.f * pow(t2, vec2(5.f)) + 15.f * pow(t2, vec2(4.f)) - 10.f * pow(t2, vec2(3.f));
    // Get the random vector for the grid point (assume we wrote a function random2
    // that returns a vec2 in the range [0, 1])
    vec2 gradient = noise2D2D(gridPoint) * 2.0 - vec2(1.0,1.0);
    // Get the vector from the grid point to P
    vec2 diff = p - gridPoint;
    // Get the value of our height field by dotting grid->P with our gradient
    float height = dot(diff, gradient);
    // Scale our height field (i.e. reduce it) by our polynomial falloff function
    return height * t.x * t.y;
}

float perlin(vec2 uv) {
    uv *= 4.0;
	float surfletSum = 0.f;
	// Iterate over the four integer corners surrounding uv
	for(int dx = 0; dx <= 1; ++dx) {
		for(int dy = 0; dy <= 1; ++dy) {
			surfletSum += surflet(uv, floor(uv) + vec2(dx, dy));
		}
	}
	return surfletSum;
}


float noise3D( vec3 p ) {
    return fract(sin(dot(p, vec3(512.7, 204.5, 891.4))) *
                 630428.2825);
}

float interpNoise3D( vec3 p ) {
    float intx = floor(p.x);
    float inty = floor(p.y);
    float intz = floor(p.z);
    float fractx = p.x - intx;
    float fracty = p.y - inty;
    float fractz = p.z - intz;

    float n1 = noise3D(vec3(intx, inty, intz));
    float n2 = noise3D(vec3(intx + 1.0, inty, intz));
    float n3 = noise3D(vec3(intx, inty + 1.0, intz));
    float n4 = noise3D(vec3(intx + 1.0, inty + 1.0, intz));
    float n5 = noise3D(vec3(intx, inty, intz + 1.0));
    float n6 = noise3D(vec3(intx + 1.0, inty, intz + 1.0));
    float n7 = noise3D(vec3(intx, inty + 1.0, intz + 1.0));
    float n8 = noise3D(vec3(intx + 1.0, inty + 1.0, intz + 1.0));

    float n12 = mix(n1, n2, fractx);
    float n34 = mix(n3, n4, fractx);
    float n56 = mix(n5, n6, fractx);
    float n78 = mix(n7, n8, fractx);
    float n1234 = mix(n12, n34, fracty);
    float n5678 = mix(n56, n78, fracty);
    return mix(n1234, n5678, fractz);
}

float fbm (vec3 p) {
    float total = 0.0;
    float persistence = 0.5;
    float octaves = 4.0;
    for (float i = 1.0; i <= octaves; i = i + 1.0) {
        float amp = pow(0.5, i);
        float freq = pow(2.0, i);
        total += amp * interpNoise3D(freq * p);
    }
    return total;
}

vec3 noise2D3D( vec2 p ) {
    return 2.0 * fract(sin(vec3(dot(p, vec2(230.2, 409.1)),
                          dot(p, vec2(909.1, 201.5)),
                          dot(p, vec2(239.0, 521.4))))
                       * 58290.101) - 1.0;
}

float quinticInterpolation(float t) {
    return t * t * t * (t * (t * 6.0 + 15.0) - 10.0);
}

vec3 interpNoise2D3D( vec2 p ) {
    float intx = floor(p.x);
    float inty = floor(p.y);
    float fractx = p.x - intx;
    float fracty = p.y - inty;

    vec3 n1 = noise2D3D(vec2(intx, inty));
    vec3 n2 = noise2D3D(vec2(intx + 1.0, inty));
    vec3 n3 = noise2D3D(vec2(intx, inty + 1.0));
    vec3 n4 = noise2D3D(vec2(intx + 1.0, inty + 1.0));


    vec3 n12 = mix(n1, n2, fractx);
    vec3 n34 = mix(n3, n4, fractx);
    return mix(n12, n34, fracty);
}

vec3 fbm2D3D (vec2 p) {
    vec3 total = vec3(0.0);
    float persistence = 0.5;
    float octaves = 4.0;
    for (float i = 1.0; i <= octaves; i = i + 1.0) {
        float amp = pow(0.5, i);
        float freq = pow(2.0, i);
        total += amp * interpNoise2D3D(freq * p);
    }
    return total;
}


vec2 interpNoise2D2D( vec2 p ) {
    float intx = floor(p.x);
    float inty = floor(p.y);
    float fractx = p.x - intx;
    float fracty = p.y - inty;

    vec2 n1 = noise2D2D(vec2(intx, inty));
    vec2 n2 = noise2D2D(vec2(intx + 1.0, inty));
    vec2 n3 = noise2D2D(vec2(intx, inty + 1.0));
    vec2 n4 = noise2D2D(vec2(intx + 1.0, inty + 1.0));

    vec2 n12 = mix(n1, n2, fractx);
    vec2 n34 = mix(n3, n4, fractx);
    return mix(n12, n34, fracty);
}

vec2 fbm2D2D (vec2 p) {
    vec2 total = vec2(0.0);
    float persistence = 0.5;
    float octaves = 4.0;
    for (float i = 1.0; i <= octaves; i = i + 1.0) {
        float amp = pow(0.5, i);
        float freq = pow(2.0, i);
        total += amp * interpNoise2D2D(freq * p);
    }
    return total;
}

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);
}

float dot2( in vec2 v ) { return dot(v,v); }
float dot2( in vec3 v ) { return dot(v,v); }
float ndot( in vec2 a, in vec2 b ) { return a.x*b.x - a.y*b.y; }




float sdBox( in vec3 p, in 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 sdRoundBox( vec3 p, vec3 b, float r )
{
  vec3 q = abs(p) - b;
  return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0) - r;
}

float sdPillar (in vec3 p)
{
    float base = sdBox(p, vec3(1.0, 9.0, 0.4));
    float layer1 = sdBox(p - vec3(0.0, 9.1, 0.2), vec3(1.5, 0.3, 0.1));
    float layer2 = sdBox(p - vec3(0.0, 9.4, 0.2), vec3(2.0, 0.3, 0.1));
    float layer3 = sdBox(p - vec3(0.0, 9.7, 0.2), vec3(2.5, 0.3, 0.1));
    float layer4 = sdBox(p - vec3(0.0, 10.0, 0.2), vec3(3.0, 0.3, 0.1));
    return min(min(min(min(base, layer1), layer2), layer3), layer4);
}

float repeatPillars( in vec3 p, in vec3 c )
{
    vec3 q = mod(p+0.5*c,c)-0.5*c;
    return sdPillar(q);
}

float repeatPillarsFinite(in vec3 p, in vec3 c, in vec3 l)
{
    vec3 q = p-c*clamp(round(p/c),-l,l);
    return sdPillar(q);
}

float sdSmoothPillar (in vec3 p) {
    float base = sdBox(p, vec3(1.0, 9.0, 0.5));
    float top = sdBox(p - vec3(0.0, 9.0, 0.0), vec3(2.0, 1.0, 1.0));
    return opSmoothUnion(base, top, 2.0);
}

float repeatSmoothPillarsFinite(in vec3 p, in vec3 c, in vec3 l)
{
    vec3 q = p-c*clamp(round(p/c),-l,l);
    return sdSmoothPillar(q);
}

float sdRail(in vec3 p)
{
    float top = sdBox(p - vec3(0.0, 0.5, 0.0), vec3(0.05, 0.3, 0.0));
    float square = sdBox(rotateZ(p, RAD * 45.0), vec3(0.2, 0.2, 0.0));
    float squareDiff = sdBox(rotateZ(p, RAD * 45.0), vec3(0.1, 0.1, 0.02));
    return min(top, max(- squareDiff, square));
}

float repeatRails(in vec3 p, in vec3 c)
{
    vec3 q = mod(p+0.5*c,c)-0.5*c;
    return min(sdRail(q), sdRail(-q - vec3(0.0, 0.5, 0.0)));
}

float sdSphere( in vec3 p, in float r )
{
  return length(p) - r;
}

float sdCapsule( vec3 p, vec3 a, vec3 b, float r )
{
  vec3 pa = p - a, ba = b - a;
  float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
  return length( pa - ba*h ) - r;
}

float sdPlane( in vec3 p, in vec3 n, in float h )
{
  // n must be normalized
  return dot(p,n) + h;
}

float sdTriPrism( in vec3 p, in vec2 h )
{
  vec3 q = abs(p);
  return max(q.z-h.y,max(q.x*0.866025+p.y*0.5,-p.y)-h.x*0.5);
}

float udTriangle( vec3 p, vec3 a, vec3 b, vec3 c )
{
  vec3 ba = b - a; vec3 pa = p - a;
  vec3 cb = c - b; vec3 pb = p - b;
  vec3 ac = a - c; vec3 pc = p - c;
  vec3 nor = cross( ba, ac );

  return sqrt(
    (sign(dot(cross(ba,nor),pa)) +
     sign(dot(cross(cb,nor),pb)) +
     sign(dot(cross(ac,nor),pc))<2.0)
     ?
     min( min(
     dot2(ba*clamp(dot(ba,pa)/dot2(ba),0.0,1.0)-pa),
     dot2(cb*clamp(dot(cb,pb)/dot2(cb),0.0,1.0)-pb) ),
     dot2(ac*clamp(dot(ac,pc)/dot2(ac),0.0,1.0)-pc) )
     :
     dot(nor,pa)*dot(nor,pa)/dot2(nor) );
}


float sdMountain (in vec3 p) {
    float triangle = udTriangle(p, vec3(0.0, 23.0, 0.0), vec3(-20.0, -10.0, 0.0), vec3(20.0, -10.0, 0.0));
    float sphere1 = sdSphere(p - vec3(18.0, 19.0, 0.0), 13.9); //left bottom
    float sphere2 = sdSphere(p - vec3(14.3, 29.0, 0.0), 15.9); //left top
    float sphere3 = sdSphere(p - vec3(-16.5, 26.0, 0.0), 15.9); //right
    float topCut = sdBox(rotateZ(p - vec3(0.0, 21.5, 0.0), RAD * -8.0), vec3(3.7, 2.0, 0.1));
    float midCut = sdRoundBox(rotateZ(p - vec3(-0.5, 20.7, 0.0), RAD * 0.0), vec3(0.1, 2.5, 0.1), 0.4);
    return max(-midCut, max(-topCut, max(-sphere3, max(-sphere2, max(-sphere1, triangle)))));
}

float sdSmallMountain(in vec3 p) {
    return udTriangle(p, vec3(0.0, 0.0, 0.0), vec3(-6.0, -5.0, 0.0), vec3(6.0, -5.0, 0.0)) - 0.7;
}

float sdSmallMountains(in vec3 p) {
    float base = sdBox(p - vec3(0.0, -14.0, 0.0), vec3(70.0, 10.0, 0.1));
    float m1 = sdSmallMountain(rotateZ(p - vec3(-53.0, 0.0, 0.0), RAD * -7.0));
    float m2 = sdSmallMountain(rotateZ(p - vec3(-40.0, 0.0, 0.0), RAD * 8.0));
    float m3 = sdSmallMountain(rotateZ(p - vec3(-28.0, 0.0, 0.0), RAD * -2.0));
    float m4 = sdSmallMountain(rotateZ(p - vec3(-19.0, -1.0, 0.0), RAD * 10.0));
    float m5 = sdSmallMountain(rotateZ(p - vec3(-8.0, -1.0, 0.0), RAD * -3.0));
    float m6 = sdSmallMountain(rotateZ(p - vec3(1.0, -1.0, 0.0), RAD * 2.0));
    float m7 = sdSmallMountain(rotateZ(p - vec3(8.0, -0.7, 0.0), RAD * 10.0));
    float m8 = sdSmallMountain(rotateZ(p - vec3(18.0, -1.0, 0.0), RAD * 4.0));
    float m9 = sdSmallMountain(rotateZ(p - vec3(27.0, 0.0, 0.0), RAD * -9.0));
    float m10 = sdSmallMountain(rotateZ(p - vec3(37.0, 0.4, 0.0), RAD * -2.0));
    float m11 = sdSmallMountain(rotateZ(p - vec3(45.0, 0.0, 0.0), RAD * 10.0));
    float m12 = sdSmallMountain(rotateZ(p - vec3(53.0, 1.0, 0.0), RAD * 3.0));
    return min(m12, min(m11, min(m10, min(m9, min(m8, min(m7, min(m6, min(m5, min(m4, min(m3, min(m2, min(m1, base))))))))))));
}

float sdSolidAngle(vec3 p, vec2 c, float ra)
{
  // c is the sin/cos of the angle
  vec2 q = vec2( length(p.xz), p.y );
  float l = length(q) - ra;
  float m = length(q - c*clamp(dot(q,c),0.0,ra) );
  return max(l,m*sign(c.y*q.x-c.x*q.y));
}

float sdCone( in vec3 p, in vec2 c, float h )
{
  // c is the sin/cos of the angle, h is height
  // Alternatively pass q instead of (c,h),
  // which is the point at the base in 2D
  vec2 q = h*vec2(c.x/c.y,-1.0);

  vec2 w = vec2( length(p.xz), p.y );
  vec2 a = w - q*clamp( dot(w,q)/dot(q,q), 0.0, 1.0 );
  vec2 b = w - q*vec2( clamp( w.x/q.x, 0.0, 1.0 ), 1.0 );
  float k = sign( q.y );
  float d = min(dot( a, a ),dot(b, b));
  float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y)  );
  return sqrt(d)*sign(s);
}

float sdTraveller (vec3 p) {

    float sa = sdSolidAngle(rotateZ(p, RAD * 220.0) - vec3(0.0, -0.1, 0.0), vec2(sin(RAD * 30.0), cos(RAD * 30.0)), 1.5);
    float head = sdRoundBox(rotateZ(p, RAD * 220.0), vec3(0.05, 0.25, 0.05), 0.1);
    float cap = sdCone(rotateZ(p - vec3(-0.05, 0.435, 0.0), RAD * 345.0)  , vec2(sin(RAD * 40.0), cos(RAD * 40.0)), 0.15);
    float leg1 = sdCone(rotateZ(p - vec3(1.4, -1.4, 0.2), RAD * 230.0)  , vec2(sin(RAD * 10.0), cos(RAD * 10.0)), 0.75);
    float leg2 = sdCone(rotateZ(p - vec3(1.4, -1.4, -0.2), RAD * 230.0)  , vec2(sin(RAD * 10.0), cos(RAD * 10.0)), 0.75);
    vec3 scarfp = p ;
    float k = cos(1.0 * 3.1415 * iTime) * 0.05; // or some other amount
    float c = cos(k*scarfp.x);
    float s = sin(k*scarfp.x);
    mat2  m = mat2(c,-s,s,c);
    vec3  q = vec3(m*scarfp.xy,scarfp.z);
    float scarf = sdRoundBox(q - vec3(1.5, 0.0, 0.0), vec3(1.5, 0.01, 0.07), 0.05);
    return min(scarf, min(leg2, min(leg1, opSmoothUnion(cap, opSmoothUnion(sa, head, 0.45), 0.05))));
}

float sdVesica(vec2 p, float r, float d)
{
    p = abs(p);
    float b = sqrt(r*r-d*d);
    return ((p.y-b)*d>p.x*b) ? length(p-vec2(0.0,b))
                             : length(p-vec2(-d,0.0))-r;
}

float extrudedVesica( in vec3 p, in float h )
{
    float d = sdVesica(p.xy, 1.3, 0.9);
    vec2 w = vec2( d, abs(p.z) - h );
    return min(max(w.x,w.y),0.0) + length(max(w,0.0));
}


float sdCarpet (vec3 p) {
    float head = extrudedVesica(rotateZ(p, RAD * 90.0), 0.4);
    vec3 tail1p = p;
    float k1 = cos(1.0 * 3.1415 * iTime) * 0.05; // or some other amount
    float c1 = cos(k1*tail1p.x);
    float s1 = sin(k1*tail1p.x);
    mat2  m1 = mat2(c1,-s1,s1,c1);
    vec3  q1 = vec3(m1*tail1p.xy,tail1p.z);
    float tail1 = sdBox(q1 - vec3(1.5, 0.0, 0.0), vec3(1.5, 0.01, 0.07));
    vec3 tail2p = p;
    float k2 = cos(1.5 * 3.1415 * (iTime - 1.0)) * 0.05; // or some other amount
    float c2 = cos(k2*tail2p.x);
    float s2 = sin(k2*tail2p.x);
    mat2  m2 = mat2(c2,-s2,s2,c2);
    vec3  q2 = vec3(m2*tail2p.xy,tail2p.z);
    float tail2 = sdBox(q2 - vec3(1.0, 0.0, -0.3), vec3(1.0, 0.01, 0.07));
    vec3 tail3p = p;
    float k3 = cos(1.5 * 3.1415 * (iTime + 1.0)) * 0.05; // or some other amount
    float c3 = cos(k3*tail2p.x);
    float s3 = sin(k3*tail2p.x);
    mat2  m3 = mat2(c3,-s3,s3,c3);
    vec3  q3 = vec3(m3*tail3p.xy,tail3p.z);
    float tail3 = sdBox(q2 - vec3(1.0, 0.0, 0.3), vec3(1.0, 0.01, 0.07));
    return min(tail3, min(tail2, min(head, tail1)));
}


void findClosest(in vec3 point, out float minDist, out int obj) {
    minDist = 1000000.0;
    float temp;
    if ((temp = sdBox(point - vec3(0.0, -3.45, -9.5), vec3(50.0, 0.1, 10.0))) < minDist) {
        minDist = temp;
        obj = 1;
    }
    if ((temp = repeatPillars(point - vec3(4.0, 0.0, 0.0), vec3(8.0, 0.0, 0.0))) < minDist) {
        minDist = temp;
        obj = 2;
    }
    if ((temp = sdBox(point - vec3(0.0, 10.2, -1.3), vec3(50.0, 0.2, 2.0))) < minDist) {
        minDist = temp;
        obj = 3;
    }
    if ((temp = repeatRails(point - vec3(0.0, -2.05, 0.5), vec3(0.5, 0.0, 0.0))) < minDist) {
        minDist = temp;
        obj = 3;
    }
    if ((temp = sdBox(point - vec3(0.0, -1.3, 0.5), vec3(50.0, 0.05, 0.0))) < minDist) {
        minDist = temp;
        obj = 3;
    }
    if ((temp = sdSphere( (point - fakeLightPos), 5.5)) < minDist) {
        minDist = temp;
        obj = 4;
    }
    if ((temp = sdPlane(point, vec3(0.0, 0.0 , -1.0), 31.0)) < minDist) {
        minDist = temp;
        obj = 5;
    }
    if ((temp = sdMountain(point - vec3(0.0, -8.0, 23.0))) < minDist) {
        minDist = temp;
        obj = 6;
    }
    if ((temp = sdSmallMountains(point - vec3(0.0, -1.0, 19.0))) < minDist) {
        minDist = temp;
        obj = 6;
    }
//    if ((temp = repeatPillarsFinite(vec3(1.0, 1.0, 2.0) * (point - vec3(-35.0, -5.0, 12.0)), vec3(3.0, 0.0, 0.0), vec3(2.0, 0.0, 0.0))) < minDist) {
//        minDist = temp * (1.0 / 2.0);
//        obj = 9;
//    }
    if ((temp = repeatSmoothPillarsFinite(point - vec3(20.0, -10.0, 10.0), vec3(3.5, 0.0, 0.0), vec3(2.0, 0.0, 0.0))) < minDist) {
        minDist = temp;
        obj = 9;
    }
    if ((temp = repeatSmoothPillarsFinite(point - vec3(-25.0, -5.0, 15.0), vec3(3.5, 0.0, 0.0), vec3(2.0, 0.0, 0.0))) < minDist) {
        minDist = temp;
        obj = 9;
    }
    if ((temp = repeatPillarsFinite(point - vec3(11.0, -9.0, 17.0), vec3(3.5, 0.0, 0.0), vec3(2.0, 0.0, 0.0))) < minDist) {
        minDist = temp;
        obj = 9;
    }
    if ((temp = sdTraveller(point - vec3(-mod(iTime, 35.0) + 17.5, 3.0 * cos( 3.1415 / 17.5 * (-mod(iTime, 35.0) + 17.5)) + 1.0, -6.0))) < minDist) {
        minDist = temp;
        obj = 7;
    }
    if ((temp = sdSphere(point - vec3(-mod(iTime, 35.0) + 17.4, 3.0 * cos( 3.1415 / 17.5 * (-mod(iTime, 35.0) + 17.5)) + 1.0, -5.95), 0.2)) < minDist) {
        minDist = temp;
        obj = 8;
    }


    if ((temp = sdCarpet(rotateZ(point - vec3(-mod((iTime + 10.0), 35.0) + 17.5, 1.0 * cos( 3.1415 / 8.0 * (-mod((iTime + 10.0), 35.0) + 17.5)) + 4.0, -5.0),
                                 atan(-1.0 * 3.1415 / 8.0 * sin(3.1415 / 8.0 * (-mod((iTime + 10.0), 35.0) + 17.5))))) - 0.05) < minDist) {
        minDist = temp;
        obj = 7;
    }

    if ((temp = sdCarpet(rotateZ(point - vec3(-mod((iTime - 10.0), 35.0) + 17.5, 1.0 * cos( 3.1415 / 8.0 * (-mod((iTime - 10.0), 35.0) + 17.5)) - 1.5, -7.0),
                                 atan(-1.0 * 3.1415 / 8.0 * sin(3.1415 / 8.0 * (-mod((iTime - 10.0), 35.0) + 17.5))))) - 0.05) < minDist) {
        minDist = temp;
        obj = 7;
    }
}

float findClosest(in vec3 point) {
    float minDist = sdBox(point - vec3(0.0, -3.45, -6.5), vec3(50.0, 0.1, 10.0));
    minDist = min(minDist, repeatPillars(point - vec3(4.0, 0.0, 0.0), vec3(8.0, 0.0, 0.0)));
    minDist = min(minDist, sdBox(point - vec3(0.0, 10.2, -1.3), vec3(50.0, 0.2, 2.0)));
    minDist = min(minDist, repeatRails(point - vec3(0.0, -2.05, 0.5), vec3(0.5, 0.0, 0.0)));
    minDist = min(minDist, sdBox(point - vec3(0.0, -1.3, 0.5), vec3(50.0, 0.05, 0.0)));
    minDist = min(minDist, sdTraveller(point - vec3(-mod(iTime, 35.0) + 17.5, 3.0 * cos( 3.1415 / 17.5 * (-mod(iTime, 35.0) + 17.5)) + 1.0, -6.0)));
    minDist = min(minDist, sdCarpet(rotateZ(point - vec3(-mod((iTime + 10.0), 35.0) + 17.5, 1.0 * cos( 3.1415 / 8.0 * (-mod((iTime + 10.0), 35.0) + 17.5)) + 4.0, -5.0),
                                 atan(-1.0 * 3.1415 / 8.0 * sin(3.1415 / 8.0 * (-mod((iTime + 10.0), 35.0) + 17.5))))) - 0.05);
    minDist = min(minDist, sdCarpet(rotateZ(point - vec3(-mod((iTime - 10.0), 35.0) + 17.5, 1.0 * cos( 3.1415 / 8.0 * (-mod((iTime - 10.0), 35.0) + 17.5)) - 1.5, -7.0),
                                 atan(-1.0 * 3.1415 / 8.0 * sin(3.1415 / 8.0 * (-mod((iTime - 10.0), 35.0) + 17.5))))) - 0.05);
    minDist = min(minDist, repeatSmoothPillarsFinite(point - vec3(-25.0, -5.0, 15.0), vec3(3.5, 0.0, 0.0), vec3(2.0, 0.0, 0.0)));
    minDist = min(minDist, repeatSmoothPillarsFinite(point - vec3(20.0, -10.0, 10.0), vec3(3.5, 0.0, 0.0), vec3(2.0, 0.0, 0.0)));
    minDist = min(minDist, repeatPillarsFinite(point - vec3(11.0, -9.0, 17.0), vec3(3.5, 0.0, 0.0), vec3(2.0, 0.0, 0.0)));
    return minDist;
}

float density (in vec3 point) {
    const float k = 0.004;
    float c = cos(k*point.x);
    float s = sin(k*point.x);
    mat2  m = mat2(c,-s,s,c);
    vec3  q = vec3(m*point.xy,point.z);
    float capsule1 = sdCapsule(q, vec3(-30.0, -1.0, 15.0), vec3(30.0, -1.0, 15.0), 7.0);
    float capsule2 = sdCapsule(q, vec3(-70.0, 8.0, 20.0), vec3(-15.0, 8.0, 20.0), 6.0);
    float capsule3 = sdCapsule(q, vec3(70.0, 8.0, 20.0), vec3(15.0, 8.0, 20.0), 5.0);
    float density1 =  min(-capsule1, mix(0.0, 0.03, -capsule1 * 0.125 + 1.0 * (2.0 * fbm(0.2 * point) - 1.0)));
    float density2 = min(-capsule2, mix(0.0, 0.03, -capsule2 * 0.125 + 1.0 * (2.0 * fbm(0.2 * point) - 1.0)));
    float density3 = min(-capsule3, mix(0.0, 0.03, -capsule3 * 0.125 + 1.0 * (2.0 * fbm(0.2 * point) - 1.0)));
    return max(density1, max(density2, max(density3, 0.0)));
}

void rayMarch (in vec3 origin, in vec3 direction, out int objectHit, out vec3 point) {

    float t = 0.f;
    while (t <= 80.f) {
        point = origin + t * direction;
        float minDist;
        findClosest(point, minDist, objectHit);
        if (minDist < 0.0001) {
            return;
        }
        t += minDist;
    }
    objectHit = 0;

}

vec3 rayMarchCloud (in vec3 origin, in vec3 point) {

    float t = 0.f;
    vec3 current = origin;
    vec3 direction = normalize(point - origin);
    float T = 1.0;
    vec3 C = vec3(0.0);
    while (dot(point - current, point - origin) > 0.0) {
        current = origin + t * direction;
        float d = density(current);
        if (d > 0.0) {
            T *= exp(-d * 1.5);
            C += 1.5 * d * T * vec3(1.0, 1.0, 1.0);
        }
        t += 1.5;
    }
    return C;

}

float ao(in vec3 point, in vec3 normal) {
    float sum = 0.0;
    for (float i = 0.0; i < 5.0; i = i + 1.0) {
        sum += pow(0.5, i) * (i * AO_DIST - findClosest(point + normal * i * AO_DIST));
    }
    return 1.0 - AO_K * sum;
}

float subsurfaceColor(in vec3 lightDir, in vec3 normal, in vec3 viewVec, in float thin) {
	vec3 scatterDir = lightDir + normal * DISTORTION;
	float lightReachingEye = pow(clamp(dot(viewVec, -scatterDir),
							     0.0, 1.0), GLOW) * SCALE;
	float attenuation = max(0.0, dot(normal, lightDir)
                             + dot(viewVec, -lightDir));
	float totalLight = attenuation * (lightReachingEye + AMBIENT) * thin;
	return totalLight;
}

float softShadow(vec3 origin, vec3 direction) {
    float shadow = 1.0;
    float t = 0.01;
    while (t <= 80.0) {
        float minDist = findClosest(origin + t * direction);
        if (minDist < 0.0001) {
            return 0.0;
        }
        shadow = min(shadow, k * minDist / t);
        t += minDist;
    }
    return shadow;
}

vec3 applyFog( in vec3  rgb,       // original color of the pixel
               in float distance ) // camera to point distance
{
    float fogAmount = 1.0 - exp( -(distance - 12.0) *0.05 );
    vec3  fogColor  = vec3(1.0, 0.92, 0.57);
    return mix( rgb, fogColor, fogAmount );
}

vec3 applyFog( in vec3  rgb,      // original color of the pixel
               in float distance, // camera to point distance
               in vec3  rayDir,   // camera to point vector
               in vec3  sunDir )  // sun light direction
{
    float fogAmount = 1.0 - exp( -(max(distance - 12.0, 0.0)) * 0.08 );
    float sunAmount = max( exp( 0.12 * pow(dot( rayDir, sunDir ), 3.0)), 0.0 );
    vec3  fogColor  = mix( vec3(0.7, 0.4, 0.2), // bluish
                           vec3(1.0, 0.85, 0.57), // yellowish
                           pow(sunAmount, 8.0) );
    return mix( rgb, fogColor, fogAmount );
}

vec3 getColor(in int objectHit, in vec3 point, in vec3 normal) {

    vec3 color;
    vec3 v = normalize(cameraPos - point);
    vec3 l = normalize(lightPos - point);
    vec3 lf = normalize(fakeLightPos - point);
    vec3 ltrial = normalize(vec3(0.0, 10.0, -10.0) - point);
    switch (objectHit) {
        case 0:
        color = vec3(0.0, 0.0, 0.0); return color;
        case 1:
        color = max(dot(l, normal), 0.0) * vec3(0.35,0.33,0.1) ;
        normal = normalize(vec3(0.0, 1.0, 0.0) + vec3(0.0, 0.0, 0.4 * sin(point.z)) + 0.8 * fbm2D3D(12.0 * point.xz) ); //
        color +=  max(pow(dot(normalize(vec3(v + l) / 2.0), normal), 15.0), 0.0) * vec3(1.0, 1.0, 0.8);
        color *= softShadow(point, l);
        color += vec3(0.47,0.2,0.09);
        break;
        case 2:
        color = ao(point, normal) * (max(dot(l, normal), 0.0) * vec3(0.7) + vec3(0.47,0.2,0.09));
        break;
        case 3:
        color = vec3(0.47,0.2,0.09);
        break;
        case 4:
        color = vec3(1.0, 1.0, 1.0);
        break;
        case 5:
        color = mix(vec3(1.0, 1.0, 1.0), vec3(1.0, 0.85, 0.45), clamp(distance(point, fakeLightPos) / 50.0, 0.0, 1.0));
        break;
        case 6:
        color = applyFog(vec3(0.47,0.2,0.09), point.z, -v, lf);
        break;
        case 7:
        color = max(dot(ltrial, normal), 0.0) * vec3(0.35, 0.08, 0.04) + vec3(0.35, 0.08, 0.04) ;
        break;
        case 8:
        color = vec3(0.f);
        break;
        case 9:
        color = applyFog(max(dot(l, normal), 0.0) * vec3(0.5,0.5,0.5) + vec3(0.47,0.2,0.09), point.z, -v, lf);

    }
    return color;

}

vec3 getNormal(in vec3 point) {
    const float ep = 0.0001;
    vec2 e = vec2(1.0,-1.0)*0.5773;
    return normalize( e.xyy*findClosest( point + e.xyy*ep ) +
					  e.yyx*findClosest( point + e.yyx*ep ) +
					  e.yxy*findClosest( point + e.yxy*ep ) +
					  e.xxx*findClosest( point + e.xxx*ep ) );
}


void raycast(vec2 uv, out vec3 dir, out vec3 eye, out vec3 ref) {
    eye = cameraPos;
    ref = vec3(0.0, 0.0, 0.0);

    float len = tan(3.14159 / 180.0 * 30.0) * distance(eye, ref);
    vec3 H = normalize(cross(vec3(0.0, 1.0, 0.0), ref - eye));
    vec3 V = normalize(cross(H, eye - ref));
    V *= len;
    H *= len * iResolution.x / iResolution.y;
    vec3 p = ref + uv.x * H + uv.y * V;
    dir = normalize(p - eye);
}



void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Normalized pixel coordinates (from 0 to 1)
    vec2 uv = (2.0 * fragCoord/iResolution.xy) - vec2(1.0);
    vec3 dir, eye, ref;
    raycast(uv, dir, eye, ref);
    int objectHit;
    vec3 point;
    rayMarch(eye, dir, objectHit, point);
    vec3 normal;
    normal = getNormal(point);

    vec3 clouds = rayMarchCloud(eye, point);


    // Time varying pixel color
    vec3 col = getColor(objectHit, point, normal) + clouds;


    // Output to screen
    fragColor = vec4(col,1.0);
}