#define PI 3.14159265
#define EPSILON 0.00001

vec3 ambientColor = vec3(0.2, 0.2, 0.2);
vec3 diffuseColor = vec3(0.2, 0.6, 0.8);
vec3 specularColor = vec3(1.0, 1.0, 1.0);
vec3 lightDir = normalize(vec3(0.0, 4.0, 5.0));
bool mode2d;

vec3 mymod(vec3 p, vec3 m)
{
    return mod(p,m)-m/2.;
}

// https://commons.wikimedia.org/wiki/File:Dynamical_plane_with_Julia_set_for_c%3D0.35_with_binary_decomposition.png#cite_note-1

// compute color of pixel
float color(vec2 z, float time)
{
    const int iMax=100;
    vec2 c = vec2(0.35+sin(time/3.33)*0.7,0.0+sin(time/2.0));  // initial value
    float er2 = 10000.0;

    // iteration
    for ( int i = 0; i < iMax; i++)
    {
        z = vec2(z.x*z.x-z.y*z.y,2.0*z.x*z.y) +  c; // z= z^2+c
        if (dot(z,z)> er2)   // escape test
            // exterior
            if (z.x>0.0)
            {
                return (float(i)/1000.);
            } // upper part of the target set
       		else
                return (0.); //false; //lower part of the target set
    }
    return (0.); //false; //interior
}

vec3 palette( in float t, in vec3 a, in vec3 b, in vec3 c, in vec3 d )
{
    return a + b*cos( (2.0*PI)*(c*t+d) );
}

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

// Cylinder standing upright on the xz plane
float fCylinder(vec3 p, float r, float height)
{
	float d = length(p.xz) - r;
	d = max(d, abs(p.y) - height);
	return d;
}

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

float pModPolar(inout vec2 p, float repetitions) {
	float angle = 2.0*PI/repetitions;
	float a = atan(p.y, p.x) + angle/2.;
	float r = length(p);
	float c = floor(a/angle);
	a = mod(a,angle) - angle/2.;
	p = vec2(cos(a), sin(a))*r;
	// For an odd number of repetitions, fix cell index of the cell in -x direction
	// (cell index would be e.g. -5 and 5 in the two halves of the cell):
	if (abs(c) >= (repetitions/2.0)) c = abs(c);
	return c;
}
vec3 myFold(vec3 z)
{
    if(z.x+z.y<0.) z.xy = -z.yx; // fold 1
    if(z.x+z.z<0.) z.xz = -z.zx; // fold 2
    if(z.y+z.z<0.) z.zy = -z.yz; // fold 3

    return z;
}

vec3 twist(vec3 p, float tw)
{
    float c = cos(tw);
    float s = sin(tw);
    mat2  m = mat2(c,-s,s,c);
    vec3  q = vec3(m*p.xy,p.z);

    return q;
}

vec3 rotX(vec3 p, float a)
{
    mat3 rot = mat3(1.,0.,0.,
                    0., cos(a), -sin(a),
                   0., sin(a), cos(a));
    return p * rot;
}

vec3 rotZ(vec3 p, float a)
{
    mat3 rot = mat3(cos(a), 0., sin(a),
        			0., 1., 0.,
                   -sin(a), 0., cos(a));
    return p * rot;
}

vec3 rotY(vec3 p, float a)
{
    mat3 rot = mat3(cos(a), -sin(a), 0.,
                    sin(a), cos(a), 0.,
        			0., 0., 1.
                   );
    return p * rot;
}

// octohedron
float so(vec3 p,float s) {
     return dot(p,normalize(sign(p+1e-6)))-s;
}

//Chamfered box
float sbc(vec3 p,vec3 d,float c, float w) {
    p=max(abs(p)-d, w);
    return so(p,c);
}

float myshape(vec3 p,vec3 d,float c, float w) {
    p=max(abs(p)-d, w);
    return so(p, c);
}

float map3(vec3 p)
{
    float d;
    d =       myshape(rotZ(rotY(p  - vec3(-1.0, 0., 2.),iTime/2.),iTime/3.), vec3(0.4,  .2, 0.2 ), 0.25, 0. );
    d = min(d,myshape(rotX(rotZ(p  - vec3( 0.0, 0., 3.),iTime/4.),iTime/1.), vec3(0.15, .3, 0.15), 0.2 , 0.1));
    d = min(d,myshape(rotX(rotY(p  - vec3( 1.0, 0., 2.),iTime/3.),iTime/2.), vec3(0.15, .4, 0.15), 0.2 ,-0.10));

    return d;
}

float map4(vec3 p)
{
    return sdSphere(mymod(p, vec3(2.)), 0.5);
}

vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        map3(vec3(p.x + EPSILON, p.y, p.z) ),
        map3(vec3(p.x, p.y + EPSILON, p.z) ),
        map3(vec3(p.x, p.y, p.z  + EPSILON) )
    ));
}

vec4 ChamferedBoxes(vec2 uv)
{
    vec3 ro = vec3(0.);
    vec3 rd = normalize(vec3(2.*uv.x-1., 1.33*uv.y-0.66, 1.));
    vec3 color;
    vec4 fragColor = vec4(0.);
    float t = 0.;
    for(int i=0; i<130 && t<40.; ++i)
    {
        vec3 p = ro + t * rd;
        float d = map3(p);
        if(d<EPSILON)
        {
            vec3 norm = estimateNormal(p);

            float occ = 0.5 + 0.5 * norm.y;

            float amb = clamp(0.5 + 0.5 * norm.y, 0.0, 1.0);
            float dif = clamp(dot(lightDir, norm), 0.0, 1.0);

            vec3 h = normalize(-rd + lightDir);
            float spe = pow(clamp(dot(h, norm), 0.0, 1.0), 64.0);

            diffuseColor = vec3(map4(p+norm));
            color = amb * ambientColor
             + dif * (diffuseColor + spe * specularColor) ;

            fragColor = vec4(vec3(1.)-(color * occ),1.);

           break;
        }
        t += d;
    }

    return fragColor;
}

float map(in vec3 p, in float time)
{
	vec3 q = time>80.0?myFold(twist(p,time)):twist(p,time);
    float c = pModPolar(q.xy,8.0);
    /*
    if(c<-1.0)
        return 0.0;
    */
	float d = sdBox(mymod(q,vec3(1.0)),vec3(0.5*sin(time*3.33)+0.5,0.3,0.2));

	q.yz+=0.5*sin(-time*3.0);
    d = min(d,sdSphere(mymod(q,vec3(1.0)),0.1));

    if(!mode2d)
    {
		pModPolar(q.zy,4.0);
		d = min(d,fCylinder(mymod(q.yzx,vec3(1.0)),0.05,1.1));
    }
	return d;
}
/*
vec3 getNormal(vec3 p, float eps)
{
    vec3 n;
    n.y = map(p);
    n.x = map(vec3(p.x+eps,p.y,p.z)) - n.y;
    n.z = map(vec3(p.x,p.y,p.z+eps)) - n.y;
    n.y = eps;
    return normalize(n);
}
*/

vec4 trace(in vec3 rs, in vec3 rd, in float time)
{
    const float EPSILON_NRM = 0.1 / 1280.0;
    vec4 ret=vec4(0.0);
    float t = 0.0;
    const int maxSteps = 30;
    for(int i = 0; i < maxSteps ; ++i)
    {
        vec3 p = rs + normalize(rd) * t;
        float d = map(p, time);
        if(d < 0.000001 )
        {
            ret = vec4(palette(abs(d+1.0*p.z),vec3(0.5, 0.5, 0.5),vec3(0.5, 0.5, 0.5),vec3(1.0, 0.7, 0.4),vec3(0.00, 0.15, 0.20)),1.);
            break;
        }
        t += d*1.0;
    }

    return ret;
}

void  mainImage(out vec4 fragColor,in vec2 fragCoord)
{
    fragColor.a = 0.0;
    float myTime = iTime*1.0;

    //  camera
    vec3 dir = normalize(vec3(2.*gl_FragCoord.xy -iResolution.xy, 0.8*iResolution.y));
    vec3 org = vec3(0.,0.,myTime);

    /*if(myTime<20.0)
    {
        vec2 uv = fragCoord/iResolution.xy;
        fragColor = ChamferedBoxes(uv);
    }

    if(fragColor.a == 0.0)*/
    {
        mode2d = (myTime<30.0) || (fract(myTime/10.0)<0.5);

        fragColor = trace(org, twist(mode2d?dir:dir.xzy,myTime*0.5), myTime);
        float bindecomp = 0.0; // color(dir.xy, myTime);
        if((myTime>6.0) && bindecomp > 0.0)
        {
            fragColor.rgb = vec3(1.)-fragColor.rgb;
            fragColor.a = bindecomp;
        }
    }
}