#define ALTITUDE 0.6
#define RHOMBICTRICONTAHEDRON
#ifdef RANDOM
	#define SHARD_SEED 6
    #define SHARD_PLANE_NUM 32
    vec3 SHARD_SHAPE = vec3(.5, .5, 1.2);
#elif defined TETRAHEDRON
	#define SHARD_PLANE_NUM 4
#elif defined CUBE
	#define SHARD_PLANE_NUM 6
#elif defined OCTAHEDRON
	#define SHARD_PLANE_NUM 8
#elif defined DODECAHEDRON
	#define SHARD_PLANE_NUM 12
#elif defined ICOSCAHEDRON
	#define SHARD_PLANE_NUM 20
#elif defined RHOMBICTRICONTAHEDRON
	#define SHARD_PLANE_NUM 30
#endif
const vec3 SHARD_SPECULAR_COLOR = vec3(1., .2, .2);
const vec3 SHARD_TRANSMISSIVE_COLOR = vec3(.2, 0., 1.);
const float SHARD_OPACITY = .875;
#define SAMPLES 2
#define CAMERA_DISTANCE 2.25

struct RSet2 {
    mat2 q;
    vec2 l;
    vec4 r;
};
    RSet2 rset2() {
        return RSet2(
            mat2(1.0),
            vec2(0.0),
            vec4(0.0)
        );
    }

const vec4 rs1 = vec4(
    3.3540e-1,   7.2490e-1,   6.9690e-1,   5.2740e-1
) * vec4(400, 20, 80, 100) + vec4(1200, 20, 120, 300);
const vec4 rs2 = vec4(
    2.7730e-2,   5.1710e-1,   3.3970e-1,   1.0300e-1
) * vec4(400, 20, 80, 100) + vec4(1200, 20, 120, 300);
const vec4 rs3 = vec4(
    4.4190e-1,   7.6040e-1,   7.7920e-1,   3.7080e-1
) * vec4(400, 20, 80, 100) + vec4(1200, 20, 120, 300);

float randSin(float x, vec4 ss) {return fract(sin(x * ss.x + ss.y) * ss.z + ss.z);}

#define RAND randSin

vec3 rvec3(float x, vec4 r1, vec4 r2, vec4 r3) {
    return vec3(RAND(x, r1), RAND(x, r2), RAND(x, r3));
}

float dot2(vec3 v) {
    return dot(v, v);
}

struct Ray3 {
    // base
    vec3 b;
    // direction
    vec3 d;
};
#define eray3 Ray3(O3,O3)
#define nray3(b,d) Ray3(b,normalize(d))

Ray3 cameraRay(
    vec3 camera,
    vec3 focus,
    vec3 up,
    vec2 fov,
    float screen,
    vec2 uv
) {
    vec3 z = normalize(camera - focus);
    vec3 x = normalize(cross(up, z));
    vec3 y = normalize(cross(z, x));

    mat3 raycaster = screen * mat3(
        fov.x,0.,0.,
        0.,fov.y,0.,
        0.,0.,1.
    ) * mat3(x, y, -z);

    Ray3 ray = nray3(
        camera,
        raycaster * vec3(uv, 1.)
    );

    return ray;
}

#define Sphere3 vec4
float SphereHit(Ray3 r, Sphere3 s) {
    vec3 D = s.xyz - r.b;
    if(length(D) <= s.w)
        return 1.;
    if(dot(r.d, D) < 0.)
        return 0.;
    float D2 = dot2(D), dD = dot(r.d, D);
    float h = D2 / dD * sqrt(1. - dD * dD / D2);
    return h <= s.w ? 1. : 0.;
}

mat4 SphereTrace(Ray3 r, Sphere3 s, mat3 o) {
    vec3 D = s.xyz - r.b;
    if(length(D) <= s.w || dot(r.d, D) < 0.)
        return mat4(0.);
    float D2 = dot2(D), dD = dot(r.d, D);
    float len = dD - sqrt(dD * dD - D2 + s.w * s.w);
    vec4 p = vec4(len * r.d + r.b, 1.);
    vec4 z = vec4(normalize(p.xyz - s.xyz), 0.);
    vec4 x = vec4(normalize(cross(o[1], z.xyz)), 0.);
    vec4 y = vec4(normalize(cross(z.xyz, x.xyz)), 0.);
    return mat4(x, y, z, p);
}

mat4 SphereInnerTrace(Ray3 r, Sphere3 s, mat3 o) {
    vec3 D = s.xyz - r.b;
    float D2 = dot2(D), dD = dot(r.d, D);
    float len = dD - sqrt(dD * dD - D2 + s.w * s.w);
    vec4 p = vec4(len * r.d + r.b, 1.);
    vec4 z = vec4(normalize(s.xyz - p.xyz), 0.);
    vec4 x = vec4(normalize(cross(o[1], z.xyz)), 0.);
    vec4 y = vec4(normalize(cross(z.xyz, x.xyz)), 0.);
    return mat4(x, y, z, p);
}

float ShardHit(Ray3 r, vec4 p[SHARD_PLANE_NUM]) {
    float x, xset = 0.;
    float n, nset = 0.;
    for(int i = 0; i < SHARD_PLANE_NUM; i++) {
        vec3 pn = p[i].xyz;
        if(dot(pn, r.d) < 0.) {
            vec3 pc = pn * p[i].w;
            vec3 D = r.b - pc;
            float l = -dot(pn, D) / dot(pn, r.d);
            if(l > x || xset < .5) {
                x = l;
                xset = 1.;
            }
        }
        else {
            vec3 pc = pn * p[i].w;
            vec3 D = r.b - pc;
            float l = -dot(pn, D) / dot(pn, r.d);
            if(l < n || nset < .5) {
                n = l;
                nset = 1.;
            }
        }
    }
    return x < n ? 1. : 0.;
}

mat4 ShardTrace(Ray3 r, vec4 p[SHARD_PLANE_NUM], mat3 o) {
    float d;
    vec3 n;
    int s = 0;
    for(int i = 0; i < SHARD_PLANE_NUM; i++) {
        vec3 pn = p[i].xyz;
        if(dot(pn, r.d) < 0.) {
            vec3 pc = pn * p[i].w;
            vec3 D = r.b - pc;
            float l = -dot(pn, D) / dot(pn, r.d);
            if(l > d || s == 0) {
                d = l;
                n = pn;
                s = 1;
            }
        }
    }
    vec4 t = vec4(d * r.d + r.b, 1.);
    vec4 z = vec4(n, 0.);
    vec4 x = vec4(normalize(cross(o[1], z.xyz)), 0.);
    vec4 y = vec4(normalize(cross(z.xyz, x.xyz)), 0.);
    return mat4(x, y, z, t);
}

mat4 ShardInnerTrace(Ray3 r, vec4 p[SHARD_PLANE_NUM], mat3 o) {
    float d;
    vec3 n;
    int s = 0;
    for(int i = 0; i < SHARD_PLANE_NUM; i++) {
        vec3 pn = p[i].xyz;
        if(dot(pn, r.d) > 0.) {
            vec3 pc = pn * p[i].w;
            vec3 D = r.b - pc;
            float l = -dot(pn, D) / dot(pn, r.d);
            if(l < d || s == 0) {
                d = l;
                n = -pn;
                s = 1;
            }
        }
    }
    vec4 t = vec4(d * r.d + r.b, 1.);
    vec4 z = vec4(n, 0.);
    vec4 x = vec4(normalize(cross(o[1], z.xyz)), 0.);
    vec4 y = vec4(normalize(cross(z.xyz, x.xyz)), 0.);
    return mat4(x, y, z, t);
}

void surface(
    const in vec3 i,
    mat4 s,
    float n1,
    float n2,
    out bool tr,
    out Ray3 r,
    out float cr,
    out Ray3 t,
    out float ct
) {
    vec3 p = s[3].xyz;
    vec3 n = s[2].xyz, nx = s[0].xyz, ny = s[1].xyz;
    float Rn = dot(i, n);
    vec3 rn = Rn * n;
    vec3 rp = i - rn;
    float Rp = length(rp);

    r.b = p;
    r.d = rp - rn;

    float cosi = -Rn, sini = sqrt(1. - cosi * cosi);
    float sint = sini * n1 / n2;
    if(sint > 1.) {
        tr = true;
        t.b = p;
        t.d = vec3(0.);
        cr = 1.;
        ct = 0.;
    }
    else {
        tr = false;
        t.b = p;
        float cost = sqrt(1. - sint * sint);
        t.d = sint * rp / Rp - cost * n;
        float cs = (n1*cosi-n2*cost)/(n1*cosi+n2*cost);
        float cp = (n1*cost-n2*cosi)/(n1*cost+n2*cosi);
        cr = .5 * (cs * cs + cp * cp);
        ct = 1. - cr;
    }
}

vec3 plane(vec3 a, vec3 b, vec3 c) {
    vec3 m = a + b + c;
    vec3 n = normalize(cross(b - a, c - a));
    if(dot(m, n) < 0.)
        n *= -1.;
    return n;
}

void genPlanes(out vec4 planes[SHARD_PLANE_NUM], float R, vec3 center) {
    #ifdef RANDOM
    float seeds[SHARD_PLANE_NUM];
    for(int i = 0; i < SHARD_PLANE_NUM; i++)
        seeds[i] = float(i * SHARD_SEED) * 3.14159265358;
    for(int i = 0; i < SHARD_PLANE_NUM; i++) {
        vec3 v = normalize(2. * rvec3(seeds[i], rs1, rs2, rs3) - 1.);
        vec3 p = R * v * SHARD_SHAPE;
        v = normalize(v * SHARD_SHAPE.yzx * SHARD_SHAPE.zxy);
        planes[i] = vec4(v, dot(p, v));
    }
    #elif defined TETRAHEDRON
    vec3 i = vec3(1,0,0), j = vec3(0,1,0), k = vec3(0,0,1);
    vec3 nnn = -i-j-k, nnp = -i-j+k, npn = -i+j-k, npp = -i+j+k;
    vec3 pnn =  i-j-k, pnp =  i-j+k, ppn =  i+j-k, ppp =  i+j+k;
    R /= sqrt(3.);
    planes[0] = vec4(plane(nnp,ppp,pnn), R);
    planes[1] = vec4(plane(nnp,ppp,npn), R);
    planes[2] = vec4(plane(pnn,npn,ppp), R);
    planes[3] = vec4(plane(pnn,npn,nnp), R);
    #elif defined CUBE
    //R /= sqrt(3.);
    planes[0] = vec4( 1.,  0.,  0., R);
    planes[1] = vec4( 0.,  1.,  0., R);
    planes[2] = vec4( 0.,  0.,  1., R);
    planes[3] = vec4(-1.,  0.,  0., R);
    planes[4] = vec4( 0., -1.,  0., R);
    planes[5] = vec4( 0.,  0., -1., R);
    #elif defined OCTAHEDRON
    vec3 i = vec3(1,0,0), j = vec3(0,1,0), k = vec3(0,0,1);
    vec3 nnn = -i-j-k, nnp = -i-j+k, npn = -i+j-k, npp = -i+j+k;
    vec3 pnn =  i-j-k, pnp =  i-j+k, ppn =  i+j-k, ppp =  i+j+k;
    planes[ 0] = vec4(normalize(nnn), R);
    planes[ 1] = vec4(normalize(nnp), R);
    planes[ 2] = vec4(normalize(npn), R);
    planes[ 3] = vec4(normalize(npp), R);
    planes[ 4] = vec4(normalize(pnn), R);
    planes[ 5] = vec4(normalize(pnp), R);
    planes[ 6] = vec4(normalize(ppn), R);
    planes[ 7] = vec4(normalize(ppp), R);
    #elif defined DODECAHEDRON
    float h = .5 * (sqrt(5.) - 1.);
    float a = 1. + h, b = 1. - h * h;
    vec3 i = vec3(1,0,0), j = vec3(0,1,0), k = vec3(0,0,1);
    vec3 nnn = -i-j-k, nnp = -i-j+k, npn = -i+j-k, npp = -i+j+k;
    vec3 pnn =  i-j-k, pnp =  i-j+k, ppn =  i+j-k, ppp =  i+j+k;
    float c = 1.;
    //R /= sqrt(3.);
    planes[ 0] = vec4(plane(npp,ppp,c*vec3( 0, a, b)), R);
    planes[ 1] = vec4(plane(pnp,nnp,c*vec3( 0,-a, b)), R);
    planes[ 2] = vec4(plane(nnn,pnn,c*vec3( 0,-a,-b)), R);
    planes[ 3] = vec4(plane(ppn,npn,c*vec3( 0, a,-b)), R);
    planes[ 4] = vec4(plane(pnp,ppp,c*vec3( b, 0, a)), R);
    planes[ 5] = vec4(plane(ppn,pnn,c*vec3( b, 0,-a)), R);
    planes[ 6] = vec4(plane(nnn,npn,c*vec3(-b, 0,-a)), R);
    planes[ 7] = vec4(plane(npp,nnp,c*vec3(-b, 0, a)), R);
    planes[ 8] = vec4(plane(ppn,ppp,c*vec3( a, b, 0)), R);
    planes[ 9] = vec4(plane(npp,npn,c*vec3(-a, b, 0)), R);
    planes[10] = vec4(plane(nnn,nnp,c*vec3(-a,-b, 0)), R);
    planes[11] = vec4(plane(pnp,pnn,c*vec3( a,-b, 0)), R);
	#elif defined ICOSCAHEDRON
    float h = .5 * (sqrt(5.) - 1.);
    float a = 1. + h, b = 1. - h * h;
    vec3 i = vec3(1,0,0), j = vec3(0,1,0), k = vec3(0,0,1);
    vec3 nnn = -i-j-k, nnp = -i-j+k, npn = -i+j-k, npp = -i+j+k;
    vec3 pnn =  i-j-k, pnp =  i-j+k, ppn =  i+j-k, ppp =  i+j+k;
    planes[ 0] = vec4(normalize(nnn), R);
    planes[ 1] = vec4(normalize(nnp), R);
    planes[ 2] = vec4(normalize(npn), R);
    planes[ 3] = vec4(normalize(npp), R);
    planes[ 4] = vec4(normalize(pnn), R);
    planes[ 5] = vec4(normalize(pnp), R);
    planes[ 6] = vec4(normalize(ppn), R);
    planes[ 7] = vec4(normalize(ppp), R);
    planes[ 8] = vec4(normalize(vec3( a, b, 0)), R);
    planes[ 9] = vec4(normalize(vec3( a,-b, 0)), R);
    planes[10] = vec4(normalize(vec3(-a, b, 0)), R);
    planes[11] = vec4(normalize(vec3(-a,-b, 0)), R);
    planes[12] = vec4(normalize(vec3( 0, a, b)), R);
    planes[13] = vec4(normalize(vec3( 0, a,-b)), R);
    planes[14] = vec4(normalize(vec3( 0,-a, b)), R);
    planes[15] = vec4(normalize(vec3( 0,-a,-b)), R);
    planes[16] = vec4(normalize(vec3( b, 0, a)), R);
    planes[17] = vec4(normalize(vec3(-b, 0, a)), R);
    planes[18] = vec4(normalize(vec3( b, 0,-a)), R);
    planes[19] = vec4(normalize(vec3(-b, 0,-a)), R);
    #elif defined RHOMBICTRICONTAHEDRON
    float h = .5 * (sqrt(5.) - 1.);
    float a = 1. + h, b = 1. - h * h;
    vec3 i = vec3(1,0,0), j = vec3(0,1,0), k = vec3(0,0,1);
    vec3 nnn = -i-j-k, nnp = -i-j+k, npn = -i+j-k, npp = -i+j+k;
    vec3 pnn =  i-j-k, pnp =  i-j+k, ppn =  i+j-k, ppp =  i+j+k;
    float c = 1.;
    //R /= sqrt(3.);
    planes[ 0] = vec4(normalize(vec3( a, b, 0)+ppp), R);
    planes[ 1] = vec4(normalize(vec3( a, b, 0)+ppn), R);
    planes[ 2] = vec4(normalize(vec3( a,-b, 0)+pnp), R);
    planes[ 3] = vec4(normalize(vec3( a,-b, 0)+pnn), R);
    planes[ 4] = vec4(normalize(vec3(-a, b, 0)+npp), R);
    planes[ 5] = vec4(normalize(vec3(-a, b, 0)+npn), R);
    planes[ 6] = vec4(normalize(vec3(-a,-b, 0)+nnp), R);
    planes[ 7] = vec4(normalize(vec3(-a,-b, 0)+nnn), R);
    planes[ 8] = vec4( i, R);
    planes[ 9] = vec4(-i, R);
    planes[10] = vec4(normalize(vec3( 0, a, b)+ppp), R);
    planes[11] = vec4(normalize(vec3( 0, a, b)+npp), R);
    planes[12] = vec4(normalize(vec3( 0, a,-b)+ppn), R);
    planes[13] = vec4(normalize(vec3( 0, a,-b)+npn), R);
    planes[14] = vec4(normalize(vec3( 0,-a, b)+pnp), R);
    planes[15] = vec4(normalize(vec3( 0,-a, b)+nnp), R);
    planes[16] = vec4(normalize(vec3( 0,-a,-b)+pnn), R);
    planes[17] = vec4(normalize(vec3( 0,-a,-b)+nnn), R);
    planes[18] = vec4( j, R);
    planes[19] = vec4(-j, R);
    planes[20] = vec4(normalize(vec3( b, 0, a)+ppp), R);
    planes[21] = vec4(normalize(vec3( b, 0, a)+pnp), R);
    planes[22] = vec4(normalize(vec3(-b, 0, a)+npp), R);
    planes[23] = vec4(normalize(vec3(-b, 0, a)+nnp), R);
    planes[24] = vec4(normalize(vec3( b, 0,-a)+ppn), R);
    planes[25] = vec4(normalize(vec3( b, 0,-a)+pnn), R);
    planes[26] = vec4(normalize(vec3(-b, 0,-a)+npn), R);
    planes[27] = vec4(normalize(vec3(-b, 0,-a)+nnn), R);
    planes[28] = vec4( k, R);
    planes[29] = vec4(-k, R);
    #endif
}

#define N1 1.0003
#define N2 1.3
void ray(vec4 planes[SHARD_PLANE_NUM], Ray3 r, float weight, inout vec4 pixel) {
    //#define SPHERE_FRONT
    #ifdef SPHERE_FRONT
    if(SphereHit(r, vec4(vec3(0.), 1.)) > 0.) {
        mat4 trace = SphereTrace(r, vec4(vec3(0.), 1.), mat3(1.));
    #else
    if(ShardHit(r, planes) > 0.) {
        mat4 trace = ShardTrace(r, planes, mat3(1.));
    #endif
        Ray3 rr, rt;
        float cr, ct, ca;
        bool tr;
        surface(
            r.d,
            trace,
            N1,
            N2,
            tr,
            rr,
            cr,
            rt,
            ct
        );
        float l = 0.;
        ca = ct;
        vec3 rcol = cr * texture(iChannel0, rr.d.xzy).rgb * SHARD_SPECULAR_COLOR;
        vec3 tcol;
        r = rt;
        vec3 col = rcol;
        #define BOUNCES 8
        if(!tr)
            for(int b = 0; b < BOUNCES && ca > .001; b++) {
    			#ifdef SPHERE_REAR
                trace = SphereInnerTrace(r, vec4(vec3(0.), 1.), mat3(1.));
    			#else
                trace = ShardInnerTrace(r, planes, mat3(1.));
    			#endif
                surface(
                    r.d,
                    trace,
                    N2,
                    N1,
                    tr,
                    rr,
                    cr,
                    rt,
                    ct
                );
                l += distance(r.b, rr.b);
                vec3 c = texture(iChannel0, rt.d.xzy).rgb;
                /*
c = mix(
c,
SHARD_TRANSMISSIVE_COLOR,
l * SHARD_OPACITY
);
*/
                vec3 C = c * mix(
                    vec3(1.),
                    SHARD_TRANSMISSIVE_COLOR,
                    SHARD_OPACITY
                );
                if(tr) {
                    if(BOUNCES - b > 1)
                        tcol = vec3(0.);
                    else
                        tcol = ca * C;
                }
                else
                    tcol = ca * ct * C;
                col += tcol;
                ca *= cr;
                r = rr;
            }
        pixel += weight * vec4(col, 1.);
    }
    else
        pixel += weight * texture(iChannel0, r.d.xzy);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Normalized pixel coordinates (from 0 to 1)
    vec2 uv = fragCoord/iResolution.xy;

    float tm = iTime * 0.25;

    float R = 1.;
    vec3 center = vec3(0.);

    vec4 planes[SHARD_PLANE_NUM];
    genPlanes(planes, R, center);

    vec3
        camera = vec3(
            CAMERA_DISTANCE * vec2(cos(tm), sin(tm)),
            CAMERA_DISTANCE * ALTITUDE
            //D * 2. * vec2(cos(tm), sin(tm)),
            //D * .5
        ),
        focus = vec3(0., 0., .25),
        up = vec3(0., 0., 1.);
    vec2 fov = vec2(1.);
    float screen = 1.;

    const int SAMPLES2 = SAMPLES * SAMPLES;
    vec4 pixel = vec4(0.);
    float weight = 1. / float(SAMPLES2);
    for(int si = 0; si < SAMPLES2; si++) {
        int sx = si / SAMPLES;
        int sy = si - sx * SAMPLES;
        vec2 shift = vec2(ivec2(sx, sy) + 1) / (float(SAMPLES) * 2.);
        shift /= iResolution.x;
        Ray3 r = cameraRay(
            camera,
            focus,
            up,
            fov,
            screen,
            ((uv + shift) * 2. - 1.) * vec2(1., iResolution.y / iResolution.x)
        );
        ray(planes, r, weight, pixel);
    }

    // Output to screen
    fragColor = pixel;
    //fragColor *= 0.;
}

void mainVR( out vec4 fragColor, in vec2 fragCoord, in vec3 fragRayOri, in vec3 fragRayDir )
{
    float R = 1.;
    vec3 center = vec3(0.);

    vec4 planes[SHARD_PLANE_NUM];
    genPlanes(planes, R, center);
    Ray3 r = Ray3(
        fragRayOri,
        fragRayDir
    );
    ray(planes, r, 1., fragColor);
}