// created by florian berger (flockaroo) - 2017
// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

// moebius gears 2
// improvement of https://www.shadertoy.com/view/4d2fDc
// better performance, additional tangential-rotation, gears matching by offset instead of scaling

// feel free to play around with the defines below

#define PI2 6.2832
#define GEAR_NUM 9
#define TOOTH_NUM 15
#define GEAR_W .25
#define MOEB_R 2.
#define TSPEED 1.
#define RSPEED 1.

// iq's signed cylinder distance field
float cylinder( vec3 p, float r, float h )
{
  vec2 d = abs(vec2(length(p.xy),p.z)) - vec2(r,h);
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

float cylinderInf( vec3 p, float r )
{
  return length(p.xy)-r;
}

float gear(vec3 p)
{
    float d=10000.;
    float ang=atan(p.y,p.x);
    d=min(d,cylinder(p,1.+.1*sin(ang*float(TOOTH_NUM)),GEAR_W));
    d=max(d,-cylinderInf(p,.65));
    return d;
}

vec4 inverseQuat(vec4 q)
{
    //return vec4(-q.xyz,q.w)/length(q);
    // if already normalized this is enough
    return vec4(-q.xyz,q.w);
}

vec4 multQuat(vec4 a, vec4 b)
{
    return vec4(cross(a.xyz,b.xyz) + a.xyz*b.w + b.xyz*a.w, a.w*b.w - dot(a.xyz,b.xyz));
}

vec3 transformVecByQuat( vec3 v, vec4 q )
{
    return v + 2.0 * cross( q.xyz, cross( q.xyz, v ) + q.w*v );
}

vec4 axAng2Quat(vec3 ax, float ang)
{
    return vec4(normalize(ax),1)*sin(vec2(ang*.5)+vec2(0,PI2*.25)).xxxy;
}

// checking only nearest gears
float allGearsNear(vec3 p)
{
    float l=length(p.xy);
    float angp = -atan(p.y/l,p.x/l);
    if(angp<0.0) angp+=PI2;
    float d = 10000.0;
    float dang=PI2/(float(GEAR_NUM)-.5);
    float R=MOEB_R;
    for(float a=-1.;a<=1.;a++)
    {
      for(float b=0.;b<=1.;b++)
      {
        float ang=b*(PI2+dang*.5)+(a+floor(angp/dang))*dang;
        if(ang>PI2*2.) ang-=PI2*2.;
        if(ang<0.0) ang+=PI2*2.;
        float s=mod(floor(ang/dang+.01),2.)*2.-1.;
	    float r=R*tan(dang*.5);
        // radial offset of gears is different depending if lying or upright
        //R*(1+tan(dang*.5)*GEAR_W*.5) in upright case
        //R/cos(dang*.5) in lying case
        vec3 pos=mix(R*(1.+tan(dang*.5)*GEAR_W*.5),
                     R/cos(dang*.5),
                     -cos(ang*.5+iTime*TSPEED)*.5+.5
                    )*1.02*vec3(cos(vec2(ang,ang+PI2*.25)),0);
        vec4 quat=vec2(0,.707107).xyxy;
        quat=multQuat(quat,axAng2Quat(vec3(1,0,0),2.*sin(iTime*1.5*RSPEED)*s));
        quat=multQuat(quat,axAng2Quat(vec3(0,1,0),ang*.25+iTime*.5*TSPEED));
        quat=multQuat(quat,axAng2Quat(vec3(0,0,1),ang));
        d=min(d,gear(transformVecByQuat(p-pos,quat)/r)*r);
      }
    }
    return d;
}

float dist(vec3 p)
{
    return allGearsNear(p);
}

float march(inout vec3 pos, vec3 dir)
{
    float rval=0.0;
	float r = MOEB_R*tan(PI2/(float(GEAR_NUM)-.5)*.5);
    float h = r*sqrt(1.+4.*GEAR_W*GEAR_W);
    float R = MOEB_R+h*1.2;
    float d0=length(pos);

    // do some bounding check for better performance
    bool inside = false;
    // upper bounding plane
    inside = inside || length((pos-dir/dir.z*(pos.z-h)).xy)<sqrt(R*R-h*h);
    // lower bounding plane
    inside = inside || length((pos-dir/dir.z*(pos.z+h)).xy)<sqrt(R*R-h*h);
    // bounding sphere
    vec3 pn = pos-dir*dot(pos,dir);
    float d=length(pn);
    inside = inside || (d<R && abs((pn-dir*sqrt(R*R-d*d)).z)<h);
    if(!inside) return 0.0;

    float eps=.001;
    for(int i=0;i<80;i++)
    {
       	float d=dist(pos);
        if(d<eps) { rval=1.0; break; }
        if(d>d0+R) { rval=0.0; break; }
        pos+=dir*d*.6;
    }
    return rval;
}

vec3 getGrad(vec3 p)
{
    float eps=.01;
    vec3 d = vec3(eps,0,0);
    return vec3(
        dist(p+d.xyz)-dist(p-d.xyz),
        dist(p+d.zxy)-dist(p-d.zxy),
        dist(p+d.yzx)-dist(p-d.yzx)
        )/eps;
}

vec3 refl(vec3 dir)
{
	return texture(iChannel0,dir.yzx).zyx*1.2+.1;
}


void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
 	vec3 pos = vec3(0,0,8);
 	vec3 dir = normalize(vec3((fragCoord.xy-.5*iResolution.xy)/iResolution.x,-1.));
    vec2 ang=PI2*iMouse.xy/iResolution.xy*vec2(1,.7);
    if(iMouse.x==0. && iMouse.y==0.) ang = vec2(2.5,1.1);
    vec4 q = vec4(0,0,0,1);
    q = multQuat(q,axAng2Quat(vec3(0,0,1),ang.x));
    q = multQuat(q,axAng2Quat(vec3(1,0,0),ang.y));
    pos=transformVecByQuat(pos,q);
   	dir=transformVecByQuat(dir,q);
    float m = march(pos,dir);
    vec3 n = normalize(getGrad(pos));
    vec3 bg=mix(vec3(1),vec3(.3,.3,1),length(fragCoord.xy-.5*iResolution.xy)/iResolution.x);
    vec3 col=vec3(.9,.95,1);
    col*=refl(reflect(dir,n));
    col=.1+col*.9;
    bg = refl(dir);
    // calc ao by stepping along normal
    float ao=1.;
    ao*=.2+.8*dist(pos+n.xyz*.12)/.12;
    ao*=.4+.6*dist(pos+n.xyz*.25)/.25;
    ao*=.6+.4*dist(pos+n.xyz*.5)/.5;
    col*=clamp(ao,0.,1.);

    //ao=-dist(pos-n.xyz*.02)/.02;
    //col+=clamp(1.-ao,0.,1.);

    float vign = (1.-.3*length((fragCoord.xy/iResolution.xy)*2.-1.));
    fragColor=vec4(mix(col,bg,1.-m)*vign,1);
}