// Fork of "TUT  Sphere Gears" by iq. https://www.shadertoy.com/view/tt2XzG
// 2022-04-03 09:00:38

// Created by curiouspers - 2022
// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

// Part 1: https://www.youtube.com/watch?v=sl9x19EnKng
// Part 2: https://www.youtube.com/watch?v=bdICU2uvOdU


#define AA 1  // Set AA to 2 if you have 3090ti and 720p display

#define TAU 6.283185


mat2 Rot(float a) {
    float s=sin(a), c=cos(a);
    return mat2(c, -s, s, c);
}  


// https://iquilezles.org/articles/distfunctions
float sdRoundedCylinder( vec3 p, float ra, float rb, float h )
{
  vec2 d = vec2( length(p.xz)-2.0*ra+rb, abs(p.y) - h );
  return min(max(d.x,d.y),0.0) + length(max(d,0.0)) - rb;
}

vec2 iSphere( in vec3 ro, in vec3 rd, in vec3 cen, in float rad )
{
    ro -= cen;
    float b = dot(rd, ro);
    float c = dot(ro,ro) - rad*rad;
    float h = b*b - c;
    if ( h<0.0 ) return vec2(-1.0);
    h = sqrt(h);
    
    return vec2( -b-h, -b+h );
}
             
             
//----------------------------------


vec2 rot( vec2 v )
{
    return vec2(v.x-v.y,v.y+v.x)*0.707107;
}

float remap(float value, float min1, float max1, float min2, float max2) {
  return min2 + (value - min1) * (max2 - min2) / (max1 - min1);
}

vec4 flower( in vec3 p, float time, float mat, float segments, vec3 rotation, float ra1, float rb, float h )
{
    
    //p = fract(p*0.02);
    vec4 d = vec4(1.);
    float t = time;
    
    //float segments = 8.;
    p.yz *= Rot(rotation.x);
    p.xz *= Rot(rotation.y);
    p.xy *= Rot(rotation.z);
    
    // please teach me how to do this with domain bending and not in cycles
    float leafsPeriod = TAU/32.;
    for (float j = 0.; j < 1.; j += 1./segments){
        t+=j*0.5;
        vec3 ppp = p;
        ppp.xz *= Rot(j*TAU);
        ppp.x -= 0.2;
        //float leaf = 0.;
        for (float i = leafsPeriod*1.; i < TAU/2.6; i += leafsPeriod){
            //leaf++;
            vec3 pp = ppp;
            //pp.yz *= Rot((i+1.) * .1);

            pp.xy *= Rot(i );
            //pp.xy *= Rot(0.1);
            pp.xy *= Rot(fract(t)*leafsPeriod);

            float ra = ra1;
            
            //ra = sin(t*0.)*ra*0.5 + 0.6*ra;
            //float tt = abs(fract(t*1.)*2.-1.);
            //ra += tt*ra*i;
            //ra += ra*i;
            
            pp.x += ra+0.1;
            //pp.y -= 0.03*i;
            pp.x += pow(abs(pp.z),1.5); // bend the leaf at center
            pp.xy *= Rot(0.9);

            //float halfWay = TAU/2.5/2.;
            //float cur = i / halfWay;
            //float inc = (i < halfWay) ? ra*0.25 : -ra*0.25;
            float inc = ra*1.25 ;
            ra = ra + mix(i*inc, (i+leafsPeriod)*inc, fract(t));
            //ra = ra<0.08 ? ra : mix(ra, 0.04*(i+leafsPeriod), fract(t)); // try to shrinj last leafs
            //ra = ra<0.08 ? 0. : ra;
            
            //ra = ra + mix(leaf*inc, (leaf+1.)*inc, fract(t)*leafsPeriod);
            pp.x += mix(i*inc*2., (i+leafsPeriod)*inc*2., fract(t));
            
            //pp.y += (1.-i)*0.1;

            //pp.y += fract(t);
            //pp.x += t;
            pp.x *= 0.45;
            
            //float r1 = 0.003;
            //float r2 = 0.008;
            //float rb = r1;//mix(r1, r2, fract(t) * i / TAU);
            //float h = 0.001;
            d = min(d, vec4(sdRoundedCylinder(pp, ra, rb, h), vec3(1.)));
        }
    }
    d.w = mat;
    return d;
}
vec4 map( in vec3 p, float time )
{
    time *= .75;
    vec3 p1 = p+vec3(0.18, 0.19, -0.8);
    vec4 d = flower(p1, time, 1., 8., vec3(0.0,0.2,0.), 0.025, 0.003, 0.001);
    
    p1 = p+vec3(0.98,0.35,0.50);
    vec4 t = flower(p1, 1.-time, 2., 6., vec3(-0.32,0.12,0.2), 0.04, 0.009, 0.001);
    if( t.x<d.x ) d=t;
    
    //p *= 0.7;
    p1 = p+vec3(-0.83,0.25,0.1);
    t = flower(p1, 1.-time, 3., 9., vec3(-0.2,-0.10,-0.3), 0.025, 0.01, 0.000);
    if( t.x<d.x ) d=t;
    
	return d;
}

#define ZERO min(iFrame,0)

// https://iquilezles.org/articles/normalsSDF
vec3 calcNormal( in vec3 pos, in float time )
{
#if 0
    vec2 e = vec2(1.0,-1.0)*0.5773;
    const float eps = 0.00025;
    return normalize( e.xyy*map( pos + e.xyy*eps, time ).x + 
					  e.yyx*map( pos + e.yyx*eps, time ).x + 
					  e.yxy*map( pos + e.yxy*eps, time ).x + 
					  e.xxx*map( pos + e.xxx*eps, time ).x );
#else
    // klems's trick to prevent the compiler from inlining map() 4 times
    vec3 n = vec3(0.0);
    for( int i=ZERO; i<4; i++ )
    {
        vec3 e = 0.5773*(2.0*vec3((((i+3)>>1)&1),((i>>1)&1),(i&1))-1.0);
        n += e*map(pos+0.0005*e,time).x;
    }
    return normalize(n);
#endif    
}

float calcAO( in vec3 pos, in vec3 nor, in float time )
{
	float occ = 0.0;
    float sca = 1.0;
    for( int i=ZERO; i<5; i++ )
    {
        float h = 0.01 + 0.12*float(i)/4.0;
        float d = map( pos+h*nor, time ).x;
        occ += (h-d)*sca;
        sca *= 0.95;
    }
    return clamp( 1.0 - 3.0*occ, 0.0, 1.0 );
}


vec4 intersect( in vec3 ro, in vec3 rd, in float time )
{
    vec4 res = vec4(-1.0);
    
    vec2 vol = iSphere( ro, rd, vec3(0., 0., 0.), 1.9);
    if ( vol.y > 0.0)
    {
        // raymarch
        float t = max(0.0,vol.x);
        for( int i=0; i<128 && t<vol.y; i++ )
        {
            vec4 h = map(ro+t*rd,time);
            if( h.x<0.001 ) { res=vec4(t,h.yzw); break; }
            t += h.x;
        }
    }
    
    return res;
}

mat3 setCamera( in vec3 ro, in vec3 ta, float cr )
{
	vec3 cw = normalize(ta-ro);
	vec3 cp = vec3(sin(cr), cos(cr),0.0);
	vec3 cu = normalize( cross(cw,cp) );
	vec3 cv =          ( cross(cu,cw) );
    return mat3( cu, cv, cw );
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec3 tot = vec3(0.0);
    
    #if AA>1
    for( int m=ZERO; m<AA; m++ )
    for( int n=ZERO; n<AA; n++ )
    {
        // pixel coordinates
        vec2 o = vec2(float(m),float(n)) / float(AA) - 0.5;
        vec2 p = (2.0*(fragCoord+o)-iResolution.xy)/iResolution.y;
        float d = 0.5+0.5*sin(fragCoord.x*111.0)*sin(fragCoord.y*331.543);
        float time = iTime - 1.0*(1.0/48.0)*(float(m*AA+n)+d)/float(AA*AA);
        #else    
        vec2 p = (2.0*fragCoord-iResolution.xy)/iResolution.y;
        float time = iTime;
        #endif
        time *= 0.5;
        
	    // camera	
        //float an = 6.2831*time/40.0;
        //vec3 ta = vec3( 0.0+cos(time*2.12)*0.03, cos(time*2.23)*0.05, cos(time*3.8)*0.03 );
        //vec3 ro = ta + vec3( 1.3*cos(an), 0.5, 1.2*sin(an) );
        //vec3 ta = vec3( 0.0, 0.0, 0.0 );
        //vec3 ro = ta + vec3( 0.00, 0.8, 1.2 );
        vec3 ta = vec3(0.);//vec3( 0.0+cos(time*2.12)*0.005, cos(time*2.23)*0.005, cos(time*3.8)*0.005 );
        // good camera angle
        vec3 ro = ta + vec3( 0.1, 0.8, 1.2 );
        //vec3 ro = ta + vec3( 0.5, 1.3, 1. );
        
        
        // camera-to-world transformation
        mat3 ca = setCamera( ro, ta, 0.0 );
        
        // ray direction
        float fl = 2.0;
        vec3 rd = ca * normalize( vec3(p,fl) );
        
        // background
        vec3 col = vec3(.92+rd.y)*0.8;
        
        // raymarch geometry
        vec4 tuvw = intersect( ro, rd, time );
        if( tuvw.x>0.0 )
        {
            // shading/lighting	
            vec3 pos = ro + tuvw.x*rd;
            vec3 nor = calcNormal(pos, time);
            
            #if 1
            vec3 te = 0.5*texture( iChannel0, tuvw.yz*2.0 ).xyz+
                      0.5*texture( iChannel0, tuvw.yw*1.0 ).xyz;
            vec3 mate = 0.92*te;
            #else
            vec3 mate = vec3(1.);
            #endif
            
            if (tuvw.w < 2.){
                //mate =  tuvw.xyz*1.3-vec3(0.3,0.0,0.0);//tuvw.xyz;//vec3(1.000,1.000,1.000); // vec3(0.996,0.498,0.282);
                // purple - dark blue
                mate *= mix (vec3(0.3,0.0,0.0),vec3(1.,1.,0.6), tuvw.xyz-vec3(0.8,0.9,0.)) * 4.;
            } else if (tuvw.w < 3.){
                // golden - green
                mate *= mix (vec3(0.3,0.03,0.0),vec3(1.,1.,0.0), tuvw.xyz-vec3(0.5,0.9,0.)) * 4.5 + .2;
            } else if (tuvw.w < 4.){
                // light blue
                //mate = mix (vec3(0.1,0.03,0.2),vec3(0.0,2.,1.0), tuvw.xyz-vec3(1.0,0.7,0.0)) * 1.9;
                // red velvet
                mate *= mix (vec3(1.000,1.000,.0),vec3(1.500,0.1,0.1), tuvw.xyz-vec3(0.0,0.,0.0)) * 4.;
            }
            float len = length(pos);
            // coloring center
            //mate = mix( mate, vec3(1.0,.4,.2), 1.0-smoothstep(0.12,0.123,len) );
            
            
            vec3 f0 = mate;
            float ks = 0.5;//clamp(0.5+1.5*te.x, 0.0, 1.0);

            float focc = 1.0;
            focc = 0.55+0.45*dot( nor, normalize(pos) );
            focc *= 0.1+0.9*clamp(length(pos)/0.535,0.0,1.0);
            
            float occ = calcAO( pos, nor, time ) * focc;
            
            #if 0
                float li1 = 0.5+0.5*nor.y;
                col = 10.0*li1*mate*occ;
            #else
                col = vec3(0.0);
            #endif
            
            // top light
            {
                float dif = 0.5+0.5*nor.y;//clamp(
                //dif *= occ;
                
                vec3 ref = reflect(rd, nor);
                
                vec3 spe = vec3(1.)*smoothstep(0.5*sqrt(sqrt(focc)),0.9,ref.y);
                spe *= occ;
                
                float fre = clamp(1.0+dot(rd,nor),0.0,1.0);
                spe *= f0 + (1.0-f0)*pow(fre,1.0);
                spe *= 6.0;
                
                //vec3 lcol = vec3(0.7, 0.8, 1.1);
                vec3 lcol = vec3(1.);
                
                col += 0.05*lcol*dif*occ;
                col += ks*lcol*spe*dif;
                
            }
           
            // bottom light
            {
                float dif = clamp( 0.5-0.5*nor.y, 0.0, 1.0 );
                col += mate*dif*occ*2.*vec3(3.);//*vec3(0.7, 0.8, 1.1);
            }
           
        }
        
        // vignette
        col *= 1.0 - 0.2*dot(p,p);
        
        // gamma        
	    tot += pow(col,vec3(0.35) );
    #if AA>1
    }
    tot /= float(AA*AA);
    #endif
    
    // S-curve
    tot = clamp(tot, 0.0, 1.0);
    //tot = tot*tot*(3.0-2.0*tot); // same as below
    tot = smoothstep(0.0,1.0,tot); // make it more contrasty
    
    // dither
    tot += (1.0/512.0)*sin(fragCoord.x*111.0)*sin(fragCoord.y*321.5423);
    
    // show banding
    //tot = floor(tot*255.)/255.;
    //tot = (abs(dFdy(tot))+abs(dFdx(tot)))*200.0;
    
    fragColor = vec4( tot, 1.0 );
}