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