//Hyperbolic Dodecahedron + Icosahedron
// https://www.sciencedirect.com/topics/biochemistry-genetics-and-molecular-biology/virus-morphology

//forked from https://www.shadertoy.com/view/DlfXzB
 
#define AA         2
#define USE_NEWTON 0

#define CAVITY 0. //1.= with holes
//#define INVERT    //uncomment to see smooth polygon 

const float outerDodecahedronCircumradius = 2.;

// Sphere intersection function from IQ
vec2 sphIntersect( in vec3 ro, in vec3 rd, in vec3 ce, float ra )
{
    vec3 oc = ro - ce;
    float b = dot( oc, rd );
    vec3 qc = oc - b*rd;
    float h = ra*ra - dot( qc, qc );
    if( h<0.0 ) return vec2(-1.0); // no intersection
    h = sqrt( h );
    return vec2( -b-h, -b+h );
}

// Extended dual number
struct Dual
{
   float  f; // f
   vec2   d; // df/dt, d²f/dt²
};

const Dual w = Dual(1., vec2(0));
const float  phi = (1. + sqrt(5.)) / 2.;

Dual dMul(Dual a, Dual b)
{
    Dual res;
    res.f = a.f * b.f;
    res.d = vec2(a.d.x * b.f + a.f * b.d.x, a.f * b.d.y + a.d.y * b.f + 2. * a.d.x * b.d.x);
    return res;
}

Dual dSqr(Dual a)
{
    return dMul(a, a);
}

Dual dExp(Dual a)
{
    Dual res;
    res.f = exp(a.f);
    res.d = vec2(a.d.x*exp(a.f), (a.d.x*a.d.x+a.d.y) *exp(a.f));
    return res;
}


Dual dAdd(Dual a, Dual b)
{
   return Dual(a.f + b.f, a.d + b.d);
}

Dual dSub(Dual a, Dual b)
{
   return Dual(a.f - b.f, a.d - b.d);
}

Dual dConst(float x)
{
    return Dual(x, vec2(0));
}

Dual dLinear(float x)
{
    return Dual(0., vec2(x, 0.));
}

Dual dPow(Dual a, int n)
{
    Dual r = Dual(0., vec2(0));
    Dual b = a;
    for(int i = 0; i < 4; ++i)
    {
        if((n & (1 << i)) != 0)
            r = dAdd(r, b);
        b = dMul(b, a);
    }
    return r;
}

#define GENERATE_OVERLOADS(op) \
    Dual op(Dual a, Dual b, Dual c) { return op(op(a, b), c); } \
    Dual op(Dual a, Dual b, Dual c, Dual d) { return op(op(a, b, c), d); } \
    Dual op(Dual a, Dual b, Dual c, Dual d, Dual e) { return op(op(a, b, c, d), e); } \
    Dual op(Dual a, Dual b, Dual c, Dual d, Dual e, Dual f) { return op(op(a, b, c, d, e), f); }
    
GENERATE_OVERLOADS(dAdd)
GENERATE_OVERLOADS(dSub)
GENERATE_OVERLOADS(dMul)

// Hyperbolic dodecahedron
//e^((?^2 + 1) (x^2 + y^2 + z^2)) - e^((x ? - y)^2) - e^((x ? + y)^2) - e^((x - z ?)^2) - e^((x + z ?)^2) - e^((y ? - z)^2) - e^((y ? + z)^2) + 0
Dual dHyperDodecahedron(Dual x, Dual y, Dual z)
{
    float w = .98 +sin(iTime)*.02;
if(  sin(iTime*.3)>0.){
    return
#ifdef INVERT
        //SMOOTH DODECAHEDRON
        dAdd(dConst(-.022*CAVITY),
        dMul(dConst(.01), //invert sign and decrease step
            dSub(
                //e^((?^2 + 1) (x^2 + y^2 + z^2))  - (...)
                dExp(dMul(dConst(.5),dAdd(dAdd(dSqr(x),dSqr(y)),dSqr(z)))),
                dMul(dConst(.001),dAdd(
#else  
        //HYPERBOLIC DODECAHEDRON
        dAdd(dConst(-.063*CAVITY), //cavity
        dMul(dConst(-.01), //invert sign and decrease step
            dSub(
                //e^((?^2 + 1) (x^2 + y^2 + z^2))  - (...)
                dExp(dMul(dConst((phi*phi+1.)),dAdd(dAdd(dSqr(x),dSqr(y)),dSqr(z)))),
                dMul(dConst(w),dAdd(
#endif               
                    // e^((x ? + y)^2) +
                    dExp(dSqr(dAdd(dMul(dConst(phi),x),y))),
                    // e^((x ? - y)^2) +
                    dExp(dSqr(dAdd(dMul(dConst(-phi),x),y))),
                    // e^((y ? + x)^2) +                  
                    dExp(dSqr(dAdd(dMul(dConst(phi),y),z))),
                    // e^((y ? - x)^2) +
                    dExp(dSqr(dAdd(dMul(dConst(-phi),y),z))),                    
                    // e^((z ? + x)^2) +                   
                    dExp(dSqr(dAdd(dMul(dConst(phi),z),x))),
                    // e^((z ? - x)^2)+
                    dExp(dSqr(dAdd(dMul(dConst(-phi),z),x)))                                                         
                ))
            )            
        )); 
}else{       
    //Hyperbolic Icosahedron
    Dual dPhi2= dConst(phi*phi),dPhi2N= dConst(-phi*phi);   
    return
#ifdef INVERT
    dAdd(dConst(-.00315*CAVITY), 
    dMul(dConst(.002),dSub(
    //   exp( (phi^3 + phi^2 + 1 ) ( x^2 + y^2 + z^2 )) 
        dExp(dMul(dConst(.2),dAdd(dSqr(x),dSqr(y),dSqr(z)))),
        dMul(dConst(1e-9),dAdd(
#else
    dAdd(dConst(-.038*CAVITY), 
    dMul(dConst(-.002),dSub(
    //   exp( (phi^3 + phi^2 + 1 ) ( x^2 + y^2 + z^2 )) 
        dExp(dMul(dConst((phi*phi*phi+ phi*phi+1.)*1.),dAdd(dSqr(x),dSqr(y),dSqr(z)))),
        dMul(dConst(w*1.),dAdd(
#endif
            dAdd(
                // - exp( phi^2 ( x + y + z )^2)    
                dExp(dMul(dPhi2,dSqr(dAdd(x,y,z)))),
                // - exp( phi^2 ( -x + y + z )^2)                  
                dExp(dMul(dPhi2,dSqr(dSub(dAdd(x,y),z)))),
                // - exp( phi^2 (  x + y - z )^2)                
                dExp(dMul(dPhi2,dSqr(dSub(dAdd(y,z),x)))),
                // - exp( phi^2 ( x - y + z )^2 )             
                dExp(dMul(dPhi2,dSqr(dSub(dAdd(x,z),y))))
            ),
            dAdd( 
                  // - exp(( x + phi^2 y )^2)  
                  dExp(dSqr(dAdd(x,dMul(y,dPhi2)))),
                  // - exp(( y + phi^2 z )^2 ) 
                  dExp(dSqr(dAdd(y,dMul(z,dPhi2)))), 
                  // - exp(( z + phi^2 x  )^2) 
                  dExp(dSqr(dAdd(z,dMul(x,dPhi2)))),  
                  // - exp(( x - phi^2 y )^2 ) 
                  dExp(dSqr(dAdd(x,dMul(y,dPhi2N)))),  
                  // - exp(( y - phi^2 z )^2) 
                  dExp(dSqr(dAdd(y,dMul(z,dPhi2N)))),  
                  // - exp(( z - phi^2 x )^2)  
                  dExp(dSqr(dAdd(z,dMul(x,dPhi2N))))  
            )
        ))        
    )));
}
        
}

mat3 rotX(float a)
{
    return mat3(1., 0., 0., 0., cos(a), sin(a),  0., -sin(a), cos(a));
}

mat3 rotY(float a)
{
    return mat3(cos(a), 0., sin(a), 0., 1., 0., -sin(a), 0., cos(a));
}

vec3 traceRay(vec2 fragCoord)
{
    // Set up ray.
    vec3 ro = vec3(0., 0., outerDodecahedronCircumradius * 2.);
    vec3 rd = normalize(vec3((fragCoord.xy - iResolution.xy / 2.) / iResolution.y, -1.));

    mat3 m = mat3(1.);        
            
    if(iMouse.z > 0.)
        m = rotY((iMouse.x / iResolution.x * 2. - 1.) * 2.) *
            rotX((iMouse.y / iResolution.y * 2. - 1.) * 2.);
    else
        m = rotY(-sin(iTime / 2.)) * rotX(-sin(iTime / 3.));

    ro = m * ro;
    rd = m * rd;

    float t = 1e5;

    // Use ray-bounding-sphere intersection test to find initial estimate for root-finding.
    vec2 s = sphIntersect(ro, rd, vec3(0), outerDodecahedronCircumradius);

    bool hit = false;
    
    if(s.x > 0. && s.x < s.y)
    {
        t = s.x;
        
        for(int i = 0; i < 60; ++i)
        {
            vec3 rp = ro + rd * t;
            // Instead of the usual "divide by length of gradient", here the gradient along
            // the ray is computed directly.
            Dual res = dHyperDodecahedron(dAdd(dConst(rp.x), dLinear(rd.x)),
                                    dAdd(dConst(rp.y), dLinear(rd.y)),
                                    dAdd(dConst(rp.z), dLinear(rd.z)));
            if(t >= s.y)
                break;
                
            if(res.f > -2e-6)
            {
                hit = true;
                break;
            }
#if USE_NEWTON
            // Newton: x = x - f(x) / f'(x)
            t += min(.25, abs((res.f) / (res.d.x)));
#else
            // Halley: x = x - 2f(x)f'(x) / (2(f'(x)²) - f(x)f''(x))
            t += min(.25, abs(2. * res.f * res.d.x / (2. * res.d.x * res.d.x - res.f * res.d.y)));
#endif
        }
    }
    
    vec3 col = vec3(.06);

    if(hit)
    {
        vec3 rp = ro + rd * t;
        
        // Surface normal from gradient of scalar field.
        Dual res_x = dHyperDodecahedron(dAdd(dConst(rp.x), dLinear(1.)),
                                  dAdd(dConst(rp.y), dLinear(0.)),
                                  dAdd(dConst(rp.z), dLinear(0.)));
        Dual res_y = dHyperDodecahedron(dAdd(dConst(rp.x), dLinear(0.)),
                                  dAdd(dConst(rp.y), dLinear(1.)),
                                  dAdd(dConst(rp.z), dLinear(0.)));
        Dual res_z = dHyperDodecahedron(dAdd(dConst(rp.x), dLinear(0.)),
                                  dAdd(dConst(rp.y), dLinear(0.)),
                                  dAdd(dConst(rp.z), dLinear(1.)));
        vec3 n = normalize(vec3(res_x.d.x, res_y.d.x, res_z.d.x));

        // Base surface colour
        vec3 diff = .5+.5*cos(vec3(0,2,4) + length(rp)*10.);
        diff = mix(diff, diff * .75, mod(floor(rp.x * 4.) + floor(rp.y * 4.) + floor(rp.z * 4.), 2.));
        
        
        // Basic lighting
        
        float fr = mix(.01, 1., pow(1. - clamp(dot(rd, n), 0., 1.), 2.));
        vec3 r = reflect(-rd, n);

        vec3 l = normalize(vec3(1));
        col = diff * (dot(l, -n) * .5 + .5) * (1. - fr) * vec3(1., 1., .9);
        col += diff * mix(.25, 1., (dot(vec3(0,-1,0), -n) * .5 + .5) * (1. - fr)) * vec3(.2,.2,.4);
        col += smoothstep(.5, 1., max(0., dot(-r, l))) * fr * 1.1;
    }
    
    return clamp(col, 0., 1.);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec3 col = vec3(0);
    
    for(int y = 0; y < AA; ++y)
        for(int x = 0; x < AA; ++x)
        {
            col += traceRay(fragCoord + vec2(x, y) / float(AA));
        }

    col /= float(AA) * float(AA);

    fragColor = vec4(sqrt(col), 1.);
}