#extension GL_OES_standard_derivatives : enable

precision highp float;

uniform float time;
uniform vec2 mouse;
uniform vec2 resolution;

/*
Created by soma_arc - 2021
This work is licensed under Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported.
*/

// from Syntopia http://blog.hvidtfeldts.net/index.php/2015/01/path-tracing-3d-fractals/
vec2 rand2n(vec2 co, float sampleIndex) {
	vec2 seed = co * (sampleIndex + 1.0);
	seed+=vec2(-1,1);
	// implementation based on: lumina.sourceforge.net/Tutorials/Noise.html
	return vec2(fract(sin(dot(seed.xy ,vec2(12.9898,78.233))) * 43758.5453),
                 fract(cos(dot(seed.xy ,vec2(4.898,7.23))) * 23421.631));
}

struct Plane{
    vec3 p1;
    vec3 p2;
    vec3 p3;
    vec3 normal;
};


const vec3 BLACK = vec3(0);
const vec3 RED = vec3(1, 0, 0);
const vec3 LIGHT_POS = vec3(100, 100, 100);
const vec3 LIGHT_DIR = normalize(LIGHT_POS);
const float AMBIENT_FACTOR = 0.1;

const float EPSILON = 0.001;
const float PI = 3.14159265;
const float PI_2 = 3.14159265 / 2.;
const float RT_3 = sqrt(3.);
const float RT_3_INV = 1.0 / sqrt(3.);

const Plane PL1 = Plane(vec3(0, 5, RT_3_INV),
                vec3(1, 1, 0),
                vec3(2, 2, -RT_3_INV),
                normalize(vec3(RT_3 * 0.5, 0, 1.5)));

const Plane PL2 = Plane(vec3(0, 3, -RT_3_INV),
                        vec3(1, 3, 0),
                        vec3(2, 2, RT_3_INV),
                        normalize(vec3(RT_3 * 0.5, 0, -1.5)));
const Plane PL3 = Plane(
    vec3(-0.5, 0, 1),
    vec3(-0.5, 1, 0),
    vec3(-0.5, 2, 1),
    vec3(-1, 0, 0));

vec4 s2, s4, s6;
vec4 s2A, s4A, s6A;
vec4 s2B, s4B, s6B;
vec4 s2C, s4C, s6C;
vec4 gSpheres0, gSpheres1, gSpheres2, gSpheres3, gSpheres4, gSpheres5;
vec4 inversionSphere;

vec3 vertexes0, vertexes1, vertexes2, vertexes3, vertexes4,vertexes5, vertexes6, vertexes7;
Plane dividePlane;
vec4 convexSphere;

float gInvNum = 0.;
float gBoundingPlaneY = -9999999999.;

const float DISPLAY_GAMMA_COEFF = 1. / 2.2;
vec3 gammaCorrect(vec3 rgb) {
  return vec3((min(pow(rgb.r, DISPLAY_GAMMA_COEFF), 1.)),
              (min(pow(rgb.g, DISPLAY_GAMMA_COEFF), 1.)),
              (min(pow(rgb.b, DISPLAY_GAMMA_COEFF), 1.)));
}

vec3 Hsv2rgb(float h, float s, float v){
    vec3 c = vec3(h, s, v);
    const vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
    vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
    return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}


void computeCubeSphairahedronA(float zb, float zc) {
    float r2 = 0.5 + (zb * zc) / 3.0;
    float r4 = 0.5 + (zb * zb - zb * zc) / 3.0;
    float r6 = 0.5 + (-zb * zc + zc * zc) / 3.0;
    s2 = s2A = vec4(1. - r2, 0, 0, r2);
    s4 = s4A = vec4(-(1. - r4) * 0.5, zb, sqrt(3.) * (1. - r4) * 0.5, r4);
    s6 = s6A = vec4(-(1. - r6) * 0.5, zc, -sqrt(3.) * (1. - r6) * 0.5, r6);
    inversionSphere = vec4(-s6.x, -s6.y, s6.z, s6.w);
}

void computeCubeSphairahedronB(float zb, float zc) {
    float r2 = (3. * RT_3 + 2. * RT_3 * zb * zc) / 9.0;
    float r4 = (3. * zb * zb - 4. * zb * zc + 3.) / 9.0;
    float r6 = (3. * zc * zc - 2. * zb * zc + 6.) / 9.0;
    s2 = s2B = vec4((2. - RT_3 * r2) * 0.5, 0, r2 * 0.5, r2);
    s4 = s4B = vec4(-(1. - r4) * 0.5, zb, RT_3 * (1. - r4) * 0.5, r4);
    s6 = s6B = vec4(-(1. - r6) * 0.5, zc, -RT_3 * (1. - r6) * 0.5, r6);
    inversionSphere = vec4(-s6.x, -s6.y, s6.z, s6.w);
}

void computeCubeSphairahedronC(float zb, float zc) {
    float r2 = (zb * zb + 2. * zb * zc + 6.) / (5. * RT_3);
    float r4 = (3. * zb * zb - 4. * zb * zc + 3.) / (5. * RT_3);
    float r6 = (-zb * zb - 2. * zb * zc + 5. * zc * zc + 9.) / 15.0;
    s2 = s2C = vec4((2. - RT_3 * r2) * 0.5, 0, r2 * 0.5, r2);
    s4 = s4C = vec4(-0.5, zb, RT_3 / 2. - r4, r4);
    s6 = s6C = vec4(-(1. - r6) * 0.5, zc, -RT_3 * (1. - r6) * 0.5, r6);
    inversionSphere = vec4(-s6.x, -s6.y, s6.z, s6.w);
}


float MAX_FLOAT = 1e20;
const float THRESHOLD = 0.001;

bool intersectBoundingPlane(const vec3 n, const vec3 p,
							const vec3 rayOrigin, const vec3 rayDir,
							inout float t0, inout float t1) {
	float d = -dot(p, n);
    float v = dot(n, rayDir);
    float t = -(dot(n, rayOrigin) + d) / v;
    if(THRESHOLD < t){
		if(v < 0.) {
			t0 = max(t, t0);
			t1 = MAX_FLOAT;
		} else {
			t0 = t0;
			t1 = t;
		}
		return true;
    }
    t0 = t0;
    t1 = MAX_FLOAT;
	return (v < 0.);
}

/*
float[12] pivoting(float[12] mat, int n, int k){
    int col = k;
    float maxValue = abs(mat[k*4 + k]);
    for (int i = k + 1; i < n; i++) {
        if (abs(mat[i*4 + k]) > maxValue) {
            col = i;
            maxValue = abs(mat[i*4 + k]);
        }
    }
    if (k != col) {
        float tmp[4];
        tmp[0] = mat[col*4 + 0];
        tmp[1] = mat[col*4 + 1];
        tmp[2] = mat[col*4 + 2];
        tmp[3] = mat[col*4 + 3];
        mat[col*4 + 0] = mat[k*4 + 0];
        mat[col*4 + 1] = mat[k*4 + 1];
        mat[col*4 + 2] = mat[k*4 + 2];
        mat[col*4 + 3] = mat[k*4 + 3];
        mat[k*4 + 0] = tmp[0];
        mat[k*4 + 1] = tmp[1];
        mat[k*4 + 2] = tmp[2];
        mat[k*4 + 3] = tmp[3];
    }
    return mat;
}
*/

vec4 sphereFromPoints(vec3 p0, vec3 p1, vec3 p2, vec3 p3){
    /*
    float coefficient[12];
    for (int i = 0; i < 3; i++) {
        coefficient[i * 4 + 0] = 2. * (p[i + 1].x - p[i].x);
        coefficient[i * 4 + 1] = 2. * (p[i + 1].y - p[i].y);
        coefficient[i * 4 + 2] = 2. * (p[i + 1].z - p[i].z);
        coefficient[i * 4 + 3] = -(pow(p[i].x, 2.) + pow(p[i].y, 2.) + pow(p[i].z, 2.)) +
                pow(p[i + 1].x, 2.) + pow(p[i + 1].y, 2.) + pow(p[i + 1].z, 2.);
    }
    */
    float coefficient0, coefficient1, coefficient2, coefficient3, coefficient4;
    float coefficient5, coefficient6, coefficient7, coefficient8, coefficient9;
    float coefficient10, coefficient11;
    coefficient0 = 2. * (p1.x - p0.x);
    coefficient1 = 2. * (p1.y - p0.y);
    coefficient2 = 2. * (p1.z - p0.z);
    coefficient3 = -(pow(p0.x, 2.) + pow(p0.y, 2.) + pow(p0.z, 2.)) +
                     pow(p1.x, 2.) + pow(p1.y, 2.) + pow(p1.z, 2.);
    coefficient4 = 2. * (p2.x - p1.x);
    coefficient5 = 2. * (p2.y - p1.y);
    coefficient6 = 2. * (p2.z - p1.z);
    coefficient7 = -(pow(p1.x, 2.) + pow(p1.y, 2.) + pow(p1.z, 2.)) +
                     pow(p2.x, 2.) + pow(p2.y, 2.) + pow(p2.z, 2.);
    coefficient8 = 2. * (p3.x - p2.x);
    coefficient9 = 2. * (p3.y - p2.y);
    coefficient10 = 2. * (p3.z - p2.z);
    coefficient11 = -(pow(p2.x, 2.) + pow(p2.y, 2.) + pow(p2.z, 2.)) +
                      pow(p3.x, 2.) + pow(p3.y, 2.) + pow(p3.z, 2.);
    
    /*
    const int n = 3;
    for (int k = 0; k < n - 1; k++) {
        coefficient = pivoting(coefficient, n, k);

        float vkk = coefficient[k * 4 + k];
        for (int i = k + 1; i < n; i++) {
            float vik = coefficient[i * 4 + k];
            for (int j = k; j < n + 1; ++j) {
                coefficient[i * 4 + j] = coefficient[i*4 + j] - vik * (coefficient[k * 4 + j] / vkk);
            }
        }
    }
    
    coefficient[(n - 1) * 4 + n] = coefficient[(n - 1) * 4 + n] / coefficient[(n - 1) * 4 + n - 1];
    for (int i = n - 2; i >= 0; i--) {
        float acc = 0.0;
        for (int j = i + 1; j < n; j++) {
            acc += coefficient[i * 4 + j] * coefficient[j * 4 + n];
        }
        coefficient[i * 4 + n] = (coefficient[i * 4 + n] - acc) / coefficient[i * 4 + i];
    }
    */
        //coefficient = pivoting(coefficient, 3, 0);
    int col = 0;
    float maxValue = abs(coefficient0);
    if (abs(coefficient4) > maxValue) {
        col = 1;
        maxValue = abs(coefficient4);
    }
    if (abs(coefficient8) > maxValue) {
        col = 2;
        maxValue = abs(coefficient8);
    }
    
    if (col == 1) {
        float tmp0 = coefficient4;
        float tmp1 = coefficient5;
        float tmp2 = coefficient6;
        float tmp3 = coefficient7;
        coefficient4 = coefficient0;
        coefficient5 = coefficient1;
        coefficient6 = coefficient2;
        coefficient7 = coefficient3;
        coefficient0 = tmp0;
        coefficient1 = tmp1;
        coefficient2 = tmp2;
        coefficient3 = tmp3;
    }
    
    if (col == 2) {
        float tmp0 = coefficient8;
        float tmp1 = coefficient9;
        float tmp2 = coefficient10;
        float tmp3 = coefficient11;
        coefficient8 = coefficient0;
        coefficient9 = coefficient1;
        coefficient10 = coefficient2;
        coefficient11 = coefficient3;
        coefficient0 = tmp0;
        coefficient1 = tmp1;
        coefficient2 = tmp2;
        coefficient3 = tmp3;
    }

    
        float vkk = coefficient0;
        float vik = coefficient4;    
        coefficient4 = coefficient4 - vik * (coefficient0 / vkk);
        coefficient5 = coefficient5 - vik * (coefficient1 / vkk);
        coefficient6 = coefficient6 - vik * (coefficient2 / vkk);
        coefficient7 = coefficient7 - vik * (coefficient3 / vkk);
    
        vik = coefficient8;            
        coefficient8 = coefficient8 - vik * (coefficient0 / vkk);
        coefficient9 = coefficient9 - vik * (coefficient1 / vkk);
        coefficient10 = coefficient10 - vik * (coefficient2 / vkk);
        coefficient11 = coefficient11 - vik * (coefficient3 / vkk);

        
        //coefficient = pivoting(coefficient, 3, 1);
    col = 1;
    maxValue = abs(coefficient5);
    
        if (abs(coefficient9) > maxValue) {
            col = 2;
            maxValue = abs(coefficient9);
        }
    
    if (col == 2) {
        float tmp0 = coefficient8;
        float tmp1 = coefficient9;
        float tmp2 = coefficient10;
        float tmp3 = coefficient11;
        coefficient8 = coefficient4;
        coefficient9 = coefficient5;
        coefficient10 = coefficient6;
        coefficient11 = coefficient7;
        coefficient4 = tmp0;
        coefficient5 = tmp1;
        coefficient6 = tmp2;
        coefficient7 = tmp3;
    }

        vkk = coefficient5;
        vik = coefficient9;
            
        coefficient9 = coefficient9 - vik * (coefficient5 / vkk);
        coefficient10 = coefficient10 - vik * (coefficient6 / vkk);
        coefficient11 = coefficient11 - vik * (coefficient7 / vkk);


    
    coefficient11 = coefficient11 / coefficient10;

    float acc = 0.0;
    acc += coefficient6 * coefficient11;
        
    coefficient7 = (coefficient7 - acc) / coefficient5;
        
    acc = 0.0;
    acc += coefficient1 * coefficient7;
    acc += coefficient2 * coefficient11;

    coefficient3 = (coefficient3 - acc) / coefficient0;
    
    //vec3 center = vec3(coefficient[0 * 4 + 3], coefficient[1 * 4 + 3], coefficient[2 * 4 + 3]);
    vec3 center = vec3(coefficient3, coefficient7, coefficient11);
    float r = length(center - p0);
    return vec4(center, r);
}

vec3 invertOnPoint(vec4 sphere, vec3 p) {
    vec3 d = p - sphere.xyz;
    float len = length(d);
    return d * (sphere.r * sphere.r / (len * len)) + sphere.xyz;
}

vec4 invertOnSphere(vec4 invSphere, vec4 s) {
    float r = s.w;
    float coeffR = r * RT_3 / 3.;
    vec3 p1 = invertOnPoint(invSphere, s.xyz + vec3(coeffR, coeffR, coeffR));
    vec3 p2 = invertOnPoint(invSphere, s.xyz + vec3(-coeffR, -coeffR, -coeffR));
    vec3 p3 = invertOnPoint(invSphere, s.xyz + vec3(coeffR, -coeffR, -coeffR));
    vec3 p4 = invertOnPoint(invSphere, s.xyz + vec3(coeffR, coeffR, -coeffR));
    return sphereFromPoints(p1, p2, p3, p4);
}

vec4 invertOnPlane(vec4 invSphere, Plane p) {
    return sphereFromPoints(invertOnPoint(invSphere, p.p1),
                            invertOnPoint(invSphere, p.p2),
                            invertOnPoint(invSphere, p.p3),
                            invSphere.xyz);
}

void computeGSpheres(){
    gSpheres0 = invertOnPlane(inversionSphere, PL1);
    gSpheres1 = invertOnSphere(inversionSphere, s2);
    gSpheres2 = invertOnPlane(inversionSphere, PL2);
    gSpheres3 = invertOnSphere(inversionSphere, s4);
    gSpheres4 = invertOnPlane(inversionSphere, PL3);
    gSpheres5 = invertOnSphere(inversionSphere, s6);
    
}

vec3 computeVertex(vec4 a, vec4 b, vec4 c) {
    float AB = (dot(a.xyz, a.xyz) - dot(b.xyz, b.xyz) - a.w * a.w + b.w * b.w) * 0.5 -
               dot(a.xyz, a.xyz) + dot(a.xyz, b.xyz);
    float AC = (dot(a.xyz, a.xyz) - dot(c.xyz, c.xyz) - a.w * a.w + c.w * c.w) * 0.5 -
               dot(a.xyz, a.xyz) + dot(a.xyz, c.xyz);
    float x = -dot(a.xyz, a.xyz) - dot(b.xyz, b.xyz) + 2. * dot(a.xyz, b.xyz);
    float y = -dot(a.xyz, a.xyz) - dot(c.xyz, c.xyz) + 2. * dot(a.xyz, c.xyz);
    float z = -dot(a.xyz, a.xyz) + dot(a.xyz, b.xyz) +
               dot(a.xyz, c.xyz) - dot(b.xyz, c.xyz);
    float s = (AB * y - AC * z) / (x * y - z * z);
    float t = (AC * x - AB * z) / (x * y - z * z);
    return a.xyz + (b.xyz - a.xyz) * s + (c.xyz - a.xyz) * t;
}


/*
0, 1, 2, 
0, 3, 4, 
2, 4, 5,
0, 1, 3,
3, 4, 5,
1, 2, 5,
1, 3, 5,
0, 2, 4
*/
void computeVertexes(){
    vertexes0 = computeVertex(gSpheres0,
                                gSpheres1,
                                gSpheres2);
    vertexes1 = computeVertex(gSpheres0,
                                gSpheres3,
                                gSpheres4);
    vertexes2 = computeVertex(gSpheres2,
                                gSpheres4,
                                gSpheres5);
    vertexes3 = computeVertex(gSpheres0,
                                gSpheres1,
                                gSpheres3);
    vertexes4 = computeVertex(gSpheres3,
                                gSpheres4,
                                gSpheres5);
    vertexes5 = computeVertex(gSpheres1,
                                gSpheres2,
                                gSpheres5);
    vertexes6 = computeVertex(gSpheres1,
                                gSpheres3,
                                gSpheres5);
    vertexes7 = computeVertex(gSpheres0,
                                gSpheres2,
                                gSpheres4);
    /*
    for(int i = 0; i < 8; i++) {
        vertexes[i] = computeVertex(gSpheres[index[i*3 + 0]],
                                    gSpheres[index[i*3 + 1]],
                                    gSpheres[index[i*3 + 2]]);
    }
    */
}

Plane computePlane() {
        vec3 p1 = invertOnPoint(inversionSphere, vertexes0);
        vec3 p2 = invertOnPoint(inversionSphere, vertexes1);
        vec3 p3 = invertOnPoint(inversionSphere, vertexes2);

        vec3 v1 = p2 - p1;
        vec3 v2 = p3 - p1;
        vec3 normal = normalize(cross(v1, v2));
        if (normal.y < 0.) {
            normal = normal * -1.;
        }
        Plane p = Plane(p1, p2, p3, normal);
        
        return p;
}

vec4 computeConvexSphere(){
    return invertOnPlane(inversionSphere, dividePlane);
}

// p: center of the plane
// n: normal of the plane
bool intersectPlane(vec3 p, vec3 n, 
                    vec3 rayOrigin, vec3 rayDir,
                     inout float minDist,
                    inout vec3 intersection, inout vec3 normal){
    float d = -dot(p, n);
    float v = dot(n, rayDir);
    float t = -(dot(n, rayOrigin) + d) / v;
    if(EPSILON < t && t < minDist){
		intersection = rayOrigin + t * rayDir;
        normal = n;
        minDist = t;
        return true;
    }
    return false;
}

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

float distSphere(vec3 p, vec4 sphere){
	return distance(p, sphere.xyz) - sphere.w;
}

float distPlane(vec3 pos, vec3 p, vec3 n) {
    return dot(pos - p, n);
}

float distPrism(const vec3 pos) {
    float d = -1.;
	d = max(distPlane(pos, PL1.p1,
					  PL1.normal),
			d);
    d = max(distPlane(pos, PL2.p1,
					  PL2.normal),
			d);
    d = max(distPlane(pos, PL3.p1,
					  PL3.normal),
			d);
    return d;
}

float distInfSphairahedra(const vec3 pos) {
    float d = distPrism(pos);
    d = max(distPlane(pos, dividePlane.p1,
                           dividePlane.normal),
            d);
	d = max(-distSphere(pos, s2),
			d);
    d = max(-distSphere(pos, s4),
			d);
    d = max(-distSphere(pos, s6),
			d);
    return d;
}

void SphereInvert(inout vec3 pos, inout float dr, vec4 s) {
    vec3 diff = pos - s.xyz;
    float lenSq = dot(diff, diff);
    float k = (s.w * s.w) / lenSq;
    dr *= k; // (r * r) / lenSq
    pos = (diff * k) + s.xyz;
}

float distLimitSetTerrain(vec3 pos, out float invNum) {
    float dr = 1.;
    invNum = 0.;

    float d;
    for(int i = 0; i < 200; i++) {
        if(25 <= i) break;
        bool inFund = true;
		if(distance(pos, s2.xyz) < s2.w) {
            invNum ++;
			SphereInvert(pos, dr, s2);
			inFund = false;
		}
        if(distance(pos, s4.xyz) < s4.w) {
            invNum ++;
			SphereInvert(pos, dr, s4);
			inFund = false;
		}
        if(distance(pos, s6.xyz) < s6.w) {
            invNum ++;
			SphereInvert(pos, dr, s6);
			inFund = false;
		}

		pos -= PL1.p1;
		d = dot(pos, PL1.normal);
		if(d > 0.) {
            //invNum += 1.;
			pos -= 2. * d * PL1.normal;
			inFund = false;
		}
		pos += PL1.p1;
        
        pos -= PL2.p1;
		d = dot(pos, PL2.normal);
		if(d > 0.) {
            //invNum += 1.;
			pos -= 2. * d * PL2.normal;
			inFund = false;
		}
		pos += PL2.p1;
        
        pos -= PL3.p1;
		d = dot(pos, PL3.normal);
		if(d > 0.) {
            //invNum += 1.;
			pos -= 2. * d * PL3.normal;
			inFund = false;
		}
		pos += PL3.p1;

        if(inFund) break;
    }
    return distInfSphairahedra(pos) / abs(dr) * 0.1;
}

vec2 distFunc(vec3 p) {
	vec2 d = vec2(distLimitSetTerrain(p, gInvNum), 1);
    return d;
}

const vec2 NORMAL_COEFF = vec2(0.001, 0.);
vec3 computeNormal(const vec3 p){
  return normalize(vec3(distFunc(p + NORMAL_COEFF.xyy).x - distFunc(p - NORMAL_COEFF.xyy).x,
                        distFunc(p + NORMAL_COEFF.yxy).x - distFunc(p - NORMAL_COEFF.yxy).x,
                        distFunc(p + NORMAL_COEFF.yyx).x - distFunc(p - NORMAL_COEFF.yyx).x));
}

const int MAX_MARCH = 300;
int march (vec3 rayOrg, vec3 rayDir, inout float minDist,
           inout vec3 intersection, inout vec3 normal) {
    vec2 dist = vec2(-1);
    float rayLength = minDist;
    vec3 rayPos = rayOrg + rayDir * rayLength;

    for(int i = 0 ; i < MAX_MARCH ; i++){
        dist = distFunc(rayPos);
        rayLength += dist.x;
        rayPos = rayOrg + rayDir * rayLength;
        if(dist.x < EPSILON){
            int objId = int(dist.y);
            intersection = rayPos;
            normal = computeNormal(intersection);
            minDist = rayLength;
            return objId;
        }
    }
    return -1;
}

vec3 calcRay (const vec3 eye, const vec3 target, const vec3 up, const float fov,
              const float width, const float height, const vec2 coord){
  float imagePlane = (height * .5) / tan(fov * .5);
  vec3 v = normalize(target - eye);
  vec3 xaxis = normalize(cross(v, up));
  vec3 yaxis =  normalize(cross(v, xaxis));
  vec3 center = v * imagePlane;
  vec3 origin = center - (xaxis * (width  *.5)) - (yaxis * (height * .5));
  return normalize(origin + (xaxis * coord.x) + (yaxis * (height - coord.y)));
}

const vec4 K = vec4(1.0, .666666, .333333, 3.0);
vec3 hsv2rgb(const vec3 c){
  vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
  return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}


vec3 getMatColor(int objId, vec3 normal, vec3 intersection, float t){
    if(objId == 1){
        vec3 col = Hsv2rgb((1., -0.01 + (gInvNum) * 0.02), 1., 1.);;
        return col; 
    }
	return BLACK;
}

vec3 sky(vec3 rayDir){
	return clamp(vec3(.7, .8, 1.) + exp(dot(rayDir, LIGHT_DIR))*0.1, 0.0, 1.0);
}

float computeShadowFactor (vec3 rayOrg, vec3 rayDir,
                           float mint, float maxt, float k) {
    float shadowFactor = 1.0;
    float t = mint;
    for(int n = 0; n < 1000; n++){
        if(t > maxt) break;
        float d = distFunc(rayOrg + rayDir * t).x;
        if(d < 0.0001) return 0.;

        shadowFactor = min(shadowFactor, k * d / t);
        t += d;
    }
    return clamp(shadowFactor, 0.0, 1.0);
}


float specular(vec3 n,vec3 l,vec3 e,float s) {    
    return pow(clamp(dot(reflect(e,n),l),0.0, 1.),s);
}

const float FOG_START = 5.;
const float FOG_END = 100.;
const float FOG_END_START_RECIPROCAL = 1. / (FOG_END - FOG_START);
const vec3 FOG_F = vec3(1.);
vec3 calcColor(float time, vec3 eye, vec3 ray){
  	vec3 l = BLACK;

    float tmin = 0.;
    float t = 9999999.;
    const int maxLoop = 10;
    	vec3 intersection, normal;
    	bool isect;
    	int objId = -1;
        
        
    bool hit = intersectBoundingPlane(vec3(0, 1, 0), vec3(0, gBoundingPlaneY, 0),
                                 eye, ray,
                                 tmin, t);
    if(hit && tmin < 30.) {
        objId = march(eye, ray, tmin,
                      intersection, normal);
        if(objId != -1){
            vec3 matColor = getMatColor(objId, normal, intersection, t);
            float k = computeShadowFactor(intersection + 0.0001 * normal,
                                         LIGHT_DIR,
                                         0.001, 5., 100.);
            vec3 diffuse =  clamp(dot(normal, LIGHT_DIR), 0., 1.) * matColor;
            vec3 ambient = matColor * AMBIENT_FACTOR;
    		l += k * diffuse + ambient;
            
            l = mix(FOG_F, l, clamp((FOG_END - tmin) * FOG_END_START_RECIPROCAL, 0.1, 1.0));
        }
                      
    } else {
        l = mix( sky(ray), l, exp( -0.000000009*t * t ) );
    }
    
  	return l;
}

//w: start time
//s: duration
float scene(in float t, in float w, in float s){
    return clamp(t - w, 0.0, s) / s;
}

float expEasingIn(float t){
    return pow( 2., 13. * (t - 1.) );
}

float circEasingIn(float t){
	return -  (sqrt(1. - t*t) - 1.);
}

float circEasingOut(float t){
    return sqrt(1. - pow(t - 1., 2.));
}

float circEasingInOut(float t){
	t /= .5;
	if (t < 1.) return -.5 * (sqrt(1. - t*t) - 1.);
	t -= 2.;
	return .5 * (sqrt(1. - t*t) + 1.);
}

const vec3 target = vec3(0, 0, 0);
vec3 up = vec3(0, 1, 0);
float fov = radians(60.);

const float SAMPLE_NUM = 2.;
void main(  ){
    float t = mod(time, 31.);
    
    float rad = mix(0.0001, 1.0, circEasingInOut(scene(t, 1.0, 15.)));
    float rotationT = mix(0., 50., (scene(t, 0., 15.)));
    if(t < 21.) {
        
        float tb = rad * cos(rotationT);
        float tc = rad * sin(rotationT);
        tb = mix(tb, 0., circEasingInOut(scene(t, 15., 1.)));
        tc = mix(tc, 0., circEasingInOut(scene(t, 15., 1.)));
        
        tb = mix(tb, 1.2224464245160123, circEasingInOut(scene(t, 16., 1.)));
        tc = mix(tc, 0.612090457867328, circEasingInOut(scene(t, 16., 1.)));
        
        tb = mix(tb, -0.8661174154839888, circEasingInOut(scene(t, 17., 1.)));
        tc = mix(tc, 0., circEasingInOut(scene(t, 17., 1.)));
        
        tb = mix(tb, 0.61096041060052, circEasingInOut(scene(t, 18., 1.)));
        tc = mix(tc, -0.6110480831254523, circEasingInOut(scene(t, 18., 1.)));
        
        tb = mix(tb, -0.6144427845492827, circEasingInOut(scene(t, 19., 1.)));
        tc = mix(tc, 0.6155169219301033, circEasingInOut(scene(t, 19., 1.)));
        
        tb = mix(tb, -0.6492253932449349, circEasingInOut(scene(t, 20., 1.)));
        tc = mix(tc, -1.1670917737220714, circEasingInOut(scene(t, 20., 1.)));
        
        computeCubeSphairahedronA(tb, tc);
    } else if(t < 22.) {
        computeCubeSphairahedronA(-0.6492253932449349, -1.1670917737220714);
        computeCubeSphairahedronB(-1.4263771231354925, -0.4114778978502866);
        s2 = mix(s2A, s2B, circEasingInOut(scene(t, 21., 1.)));
        s4 = mix(s4A, s4B, circEasingInOut(scene(t, 21., 1.)));
        s6 = mix(s6A, s6B, circEasingInOut(scene(t, 21., 1.)));
        inversionSphere = vec4(-s6.x, -s6.y, s6.z, s6.w);
    } else if(t < 26.){
        float tb = -1.4263771231354925;
        float tc = -0.4114778978502866;
        tb = mix(tb, 0.9480876782578194, circEasingInOut(scene(t, 22., 1.)));
        tc = mix(tc, -0.27487295833102093, circEasingInOut(scene(t, 22., 1.)));
        
        tb = mix(tb, -0.7341647058410344, circEasingInOut(scene(t, 23., 1.)));
        tc = mix(tc, -0.8016177852310659, circEasingInOut(scene(t, 23., 1.)));
        
        tb = mix(tb, -1.1021579831206387, circEasingInOut(scene(t, 24., 1.)));
        tc = mix(tc, 0., circEasingInOut(scene(t, 24., 1.)));
        
        tb = mix(tb, 0.7341944231647511, circEasingInOut(scene(t, 25., 1.)));
        tc = mix(tc, 0.8010951070853616, circEasingInOut(scene(t, 25., 1.)));
        computeCubeSphairahedronB(tb, tc);
    } else if(t < 27.){
        computeCubeSphairahedronB(0.7341944231647511, 0.8010951070853616);
        computeCubeSphairahedronC(0.80897179163727, -0.6172090398213005);
        s2 = mix(s2B, s2C, circEasingInOut(scene(t, 26., 1.)));
        s4 = mix(s4B, s4C, circEasingInOut(scene(t, 26., 1.)));
        s6 = mix(s6B, s6C, circEasingInOut(scene(t, 26., 1.)));
        inversionSphere = vec4(-s6.x, -s6.y, s6.z, s6.w);

    } else if(t < 30.){
        float tb = 0.80897179163727;
        float tc = -0.6172090398213005;
        tb = mix(tb, -1.0724597039839807, circEasingInOut(scene(t, 27., 1.)));
        tc = mix(tc, 0.10239062161550339, circEasingInOut(scene(t, 27., 1.)));
        
        tb = mix(tb, -0.4626628191468804, circEasingInOut(scene(t, 28., 1.)));
        tc = mix(tc, -0.7985798919860417, circEasingInOut(scene(t, 28., 1.)));
        
        tb = mix(tb, 1.0842163989231046, circEasingInOut(scene(t, 29., 1.)));
        tc = mix(tc, -0.09944669059450531, circEasingInOut(scene(t, 29., 1.)));
        
        computeCubeSphairahedronC(tb, tc);
    } else {
        computeCubeSphairahedronC(1.0842163989231046, -0.09944669059450531);
        computeCubeSphairahedronA(0., 0.);
        s2 = mix(s2C, s2A, circEasingInOut(scene(t, 30., 1.)));
        s4 = mix(s4C, s4A, circEasingInOut(scene(t, 30., 1.)));
        s6 = mix(s6C, s6A, circEasingInOut(scene(t, 30., 1.)));
        inversionSphere = vec4(-s6.x, -s6.y, s6.z, s6.w);
    }

    computeGSpheres();
    computeVertexes();
    dividePlane = computePlane();
    convexSphere = computeConvexSphere();
    gBoundingPlaneY = max(gBoundingPlaneY, s2.y + .01);
    gBoundingPlaneY = max(gBoundingPlaneY, s4.y + .01);
    gBoundingPlaneY = max(gBoundingPlaneY, s6.y + .01);
    vec3 eye = vec3(2. , 2., 2. );
    eye = mix(eye, vec3(-2, 2, 2), circEasingInOut(scene(t, 15.0, 1.)));
    eye = mix(eye, vec3(-2, 2, -2), circEasingInOut(scene(t, 20.0, 1.)));
    eye = mix(eye, vec3(2, 2, -2), circEasingInOut(scene(t, 25.0, 1.)));
    eye = mix(eye, vec3(2, 2, 2), circEasingInOut(scene(t, 30.0, 1.)));
    
    vec3 sum = vec3(0);
  	for(float i = 0. ; i < SAMPLE_NUM ; i++){
    	vec2 coordOffset = rand2n(gl_FragCoord.xy, i);
          
    	vec3 ray = calcRay(eye, target, up, fov,
        	               resolution.x, resolution.y,
            	           gl_FragCoord.xy + coordOffset);
        
    	sum += calcColor(t, eye, ray);
  	}
 	vec3 col = (sum/SAMPLE_NUM);
  
    
  	gl_FragColor = vec4(gammaCorrect(col), 1.);
}