//built off of Adam's fresnel demo as skeleton code: https://www.shadertoy.com/view/Wdj3RK

const int RAY_STEPS = 256;
const float HALF_PI = 3.14159 * 0.5;
const float PI = 3.14159;
const float TWO_PI = 2.0 * 3.14159;

const int AO_SAMPLES = 256;
const float AO_DIST = 0.15;
const float FIVETAP_K = 2.0;
const float SHADOW_HARDNESS = 20.0;

const vec3 MAIN_LIGHT_DIR = normalize(vec3(1.0, 0.2, -0.4));
//------------------------------------------------------------------

float dot2( in vec2 v ) { return dot(v,v); }
float dot2( in vec3 v ) { return dot(v,v); }
float ndot( in vec2 a, in vec2 b ) { return a.x*b.x - a.y*b.y; }

float random1(vec2 p) {
  return fract(sin(dot(p, vec2(456.789, 20487145.123))) * 842478.5453);
}

float random1( vec3 p ) {
  return fract(sin(dot(p, vec3(127.1, 311.7, 191.999))) * 43758.5453);
}

vec2 random2(vec2 p) {
  return fract(sin(vec2(dot(p, vec2(127.1, 311.7)), dot(p, vec2(269.5, 183.3)))) * 43758.5453);
}

float mySmootherStep(float a, float b, float t) {
  t = t*t*t*(t*(t*6.0 - 15.0) + 10.0);
  return mix(a, b, t);
}

float interpNoise3D1(vec3 p) {
  vec3 pFract = fract(p);
  float llb = random1(floor(p));
  float lrb = random1(floor(p) + vec3(1.0,0.0,0.0));
  float ulb = random1(floor(p) + vec3(0.0,1.0,0.0));
  float urb = random1(floor(p) + vec3(1.0,1.0,0.0));

  float llf = random1(floor(p) + vec3(0.0,0.0,1.0));
  float lrf = random1(floor(p) + vec3(1.0,0.0,1.0));
  float ulf = random1(floor(p) + vec3(0.0,1.0,1.0));
  float urf = random1(floor(p) + vec3(1.0,1.0,1.0));

  float lerpXLB = mySmootherStep(llb, lrb, pFract.x);
  float lerpXHB = mySmootherStep(ulb, urb, pFract.x);
  float lerpXLF = mySmootherStep(llf, lrf, pFract.x);
  float lerpXHF = mySmootherStep(ulf, urf, pFract.x);

  float lerpYB = mySmootherStep(lerpXLB, lerpXHB, pFract.y);
  float lerpYF = mySmootherStep(lerpXLF, lerpXHF, pFract.y);

  return mySmootherStep(lerpYB, lerpYF, pFract.z);
}


float fbm(vec3 p, float octaves) {
  float amp = 0.5;
  float freq = 32.0;
  float sum = 0.0;
  float maxSum = 0.0;
  for(float i = 0.0; i < 10.0; ++i) {
    if(i == octaves)
    break;
    maxSum += amp;
    sum += interpNoise3D1(p * freq) * amp;
    amp *= 0.5;
    freq *= 2.0;
  }
  return sum / maxSum;
}

struct Intersection
{
    float t;
    vec3 color;
    vec3 p;
    int object;
};

vec3 rotateX(vec3 p, float amt) {
    return vec3(p.x, cos(amt) * p.y - sin(amt) * p.z, sin(p.y) + cos(p.z));
}
    
vec3 rotateY(vec3 p, float a) {
    return vec3(cos(a) * p.x + sin(a) * p.z, p.y, -sin(a) * p.x + cos(a) * p.z);
}

vec3 rotateZ(vec3 p, float amt) {
    return vec3(cos(amt) * p.x - sin(amt) * p.y, sin(amt) * p.x + cos(amt) * p.y, p.z);
}

void raycast(vec2 uv, out vec3 dir, out vec3 eye, out vec3 ref) {
    eye = rotateY(vec3(0.0, 0.0,10.0), HALF_PI / 2.0);
    //eye = vec3(0.0, 5.0, 5.0);
    ref = vec3(0.0, 0.0, 0.0);
    
    float len = tan(3.14159 * 0.125) * distance(eye, ref);
    vec3 H = normalize(cross(vec3(0.0, 1.0, 0.0), ref - eye));
    vec3 V = normalize(cross(H, eye - ref));
    V *= len;
    H *= len * iResolution.x / iResolution.y;
    vec3 p = ref + uv.x * H + uv.y * V;
    dir = normalize(p - eye);
}

float sphere(vec3 p, float r, vec3 c)
{
    return distance(p, c) - r;
}

float torus(vec3 p, vec2 t)
{
  vec2 q = vec2(length(p.xz)-t.x,p.y);
  return length(q) - t.y;
}

// http://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
float sdEllipsoid( in vec3 p, in vec3 c, in vec3 r )
{
    return (length( (p-c)/r ) - 1.0) * min(min(r.x,r.y),r.z);
}

// http://iquilezles.org/www/articles/smin/smin.htm
float smin( float a, float b, float k )
{
	float h = clamp( 0.5 + 0.5*(b-a)/k, 0.0, 1.0 );
	return mix( b, a, h ) - k*h*(1.0-h);
}

// http://iquilezles.org/www/articles/smin/smin.htm
float smax( float a, float b, float k )
{
	float h = clamp( 0.5 + 0.5*(b-a)/k, 0.0, 1.0 );
	return mix( a, b, h ) + k*h*(1.0-h);
}

// http://iquilezles.org/www/articles/smin/smin.htm
vec3 smax( vec3 a, vec3 b, float k )
{
	vec3 h = clamp( 0.5 + 0.5*(b-a)/k, 0.0, 1.0 );
	return mix( a, b, h ) + k*h*(1.0-h);
}

// Box with side lengths b
float box(vec3 p, vec3 b)
{
  return length(max(abs(p) - b, 0.0));
}

float infCone( vec3 p, vec2 c )
{
    // c is the sin/cos of the angle
    vec2 q = vec2( length(p.xz), -p.y );
    float d = length(q-c*max(dot(q,c), 0.0));
    return d * ((q.x*c.y-q.y*c.x<0.0)?-1.0:1.0);
}



float sdCone(vec3 p, vec3 a, vec3 b, float ra, float rb)
{
    float rba  = rb-ra;
    float baba = dot(b-a,b-a);
    float papa = dot(p-a,p-a);
    float paba = dot(p-a,b-a)/baba;

    float x = sqrt( papa - paba*paba*baba );

    float cax = max(0.0,x-((paba<0.5)?ra:rb));
    float cay = abs(paba-0.5)-0.5;

    float k = rba*rba + baba;
    float f = clamp( (rba*(x-ra)+paba*baba)/k, 0.0, 1.0 );

    float cbx = x-ra - f*rba;
    float cby = paba - f;
    
    float s = (cbx < 0.0 && cay < 0.0) ? -1.0 : 1.0;
    
    return s*sqrt( min(cax*cax + cay*cay*baba,
                       cbx*cbx + cby*cby*baba) );
}

float sdCappedCone( in vec3 p, in float h, in float r1, in float r2 )
{
    vec2 q = vec2( length(p.xz), p.y );
    
    vec2 k1 = vec2(r2,h);
    vec2 k2 = vec2(r2-r1,2.0*h);
    vec2 ca = vec2(q.x-min(q.x,(q.y < 0.0)?r1:r2), abs(q.y)-h);
    vec2 cb = q - k1 + k2*clamp( dot(k1-q,k2)/dot2(k2), 0.0, 1.0 );
    float s = (cb.x < 0.0 && ca.y < 0.0) ? -1.0 : 1.0;
    return s*sqrt( min(dot2(ca),dot2(cb)) );
}


float capsule( vec3 p, float h, float r )
{
  p.y -= clamp( p.y, 0.0, h );
  return length( p ) - r;
}

//helix spiral from comments of https://www.shadertoy.com/view/tlXSWM
float spiral(vec3 q)
{   
	float l = length(q.xz) - 2.2;
	float d = mod(atan(q.z,q.x)-q.y * 2.2, 4.0*HALF_PI) - 3.14159;
	return length(vec2(l, d)) - 1.15;
}


//failed attempt at taper function
float taperCap( in vec3 p ) 
{
    vec3 q = vec3(0.1*sqrt(p.y)*p.x, p.y, 0.1*sqrt(p.y)*p.z);
    return capsule(q, 5., 1.);
}

float sdCappedCylinder( vec3 p, float h, float r )
{
  vec2 d = abs(vec2(length(p.xz),p.y)) - vec2(h,r);
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

vec3 opCheapBend( in vec3 p )
{
    const float k = 0.6; // or some other amount
    float c = cos(k*p.x);
    float s = sin(k*p.x);
    mat2  m = mat2(c,-s,s,c);
    vec3  q = vec3(m*p.xy,p.z);
    return q;
}

//ultimately ended up using opTwist and opTwist2 (to have varying levels of curve/twist)
vec3 opTwist( in vec3 p )
{
    const float k = 1.7; // or some other amount
    float c = cos(k*p.y);
    float s = sin(k*p.y);
    mat2  m = mat2(c,-s,s,c);
    vec3  q = vec3(m*p.xz, p.y);
    q = vec3(q.x, p.y, q.y );
    //q = p;
    return q;
}

vec3 opTwist2( in vec3 p )
{
    const float k = 1.9; // or some other amount
    float c = cos(k*p.y);
    float s = sin(k*p.y);
    mat2  m = mat2(c,-s,s,c);
    vec3  q = vec3(m*p.xz, p.y);
    q = vec3(q.x, p.y, q.y );
    //q = p;
    return q;
}

vec3 axTwist( in vec3 p )
{
    const float k = 0.9; // or some other amount
    float c = cos(k*p.y);
    float s = sin(k*p.y);
    
    vec3 q = vec3(p.x*s + p.z*c - 0.2, k*p.y, p.z*s - p.x*c);
    return q;
}

float opSmoothUnion( float d1, float d2, float k ) 
{
    float h = clamp( 0.5 + 0.5*(d2-d1)/k, 0.0, 1.0 );
    return mix( d2, d1, h ) - k*h*(1.0-h); 
}


float coneTwist3( in vec3 p )
{
    //Translate(-X, 0, Z) * Rotate(Y, -45) * Twist(Y axis) * Rotate(Z, <45) * Basic Cone

    vec3 q = p + vec3(-0.4, -0.9, 2.5);
    q = rotateY(q, -PI * 0.14);
    q = opTwist(q);// - vec3(0., 1., 0.);
    q = rotateZ(q, PI * 0.17);
    //return sdCone(q, vec3(0.0,-1.5,0.0), vec3(0.0,1.1,0.0), 0.25, 0.0 );
    return infCone(q, vec2(0.173, 0.985)); 
    //return capsule(q, 2., 0.1);
    //return sdCappedCylinder(q, 0.6, 0.2);
}

float coneTwist4( in vec3 p )
{
    //Translate(-X, 0, Z) * Rotate(Y, -45) * Twist(Y axis) * Rotate(Z, <45) * Basic Cone

    //vec3 q = p + vec3(4., -6.6, 5.1);
    vec3 q = p + vec3(-0.4, 0.9, -0.2);
    
    //was experimenting - this wqs originally meant to be another twisty cone but i liked the look
    //of the resulting shape better
    q = rotateX(q, PI * 0.013); 
    
    q = rotateY(q, -PI * 0.95);
    q = opTwist2(q);// - vec3(0., 1., 0.);
    q = rotateZ(q, PI * 0.12);
    
    
    
    //return sdCone(q, vec3(0.0,-1.5,0.0), vec3(0.0,1.1,0.0), 0.25, 0.0 );
    //return infCone(q, vec2(sin(PI * 0.04), cos(PI * 0.04))); 
    return sdCone(q, vec3(0.0,-1.5,0.0), vec3(0.0,6.1,0.0), 0.75, 0.0 );
    //return capsule(q, 2., 0.1);
    //return sdCappedCylinder(q, 0.6, 0.2);
}

float torusTwist(in vec3 p) 
{
    vec3 q = p;
    //q = axTwist(q) - vec3(0., 1., 0.);
    q = opTwist(q);
    return torus(q, vec2(0.5,0.2));
}

//adapted from snail shader: https://www.shadertoy.com/view/ld3Gz2
float spiral2( in vec3 p ) 
{
    const float sc = 1.0/3.0;
    p -= vec3(0.05,0.12,-0.09);    

    p *= sc;

    vec3 q = mat3(-0.6333234236, -0.7332753384, 0.2474039592,
                   0.7738444477, -0.6034162289, 0.1924931824,
                   0.0081370606,  0.3133626215, 0.9495986813) * p;
    //vec3 q = p;
    //vec3 q = mat3(1., 0., 0.,
    //              0., 0., 1.,
    //              0., -1., 0.) * p;
             

    const float b = 0.1759;
    
    float r = length( q.xy );
    float t = atan( q.y, q.x );
 
    // https://swiftcoder.wordpress.com/2010/06/21/logarithmic-spiral-distance-field/
    float np = (log(   r)/b-t)/(TWO_PI);
    float nm = (log(0.11)/b-t)/(TWO_PI);
    float n = min(np,nm);
    
    float ni = floor( n );
    
    float r1 = exp( b * (t + TWO_PI*ni));
    float r2 = r1 * 3.019863;
    
    //-------

    float h1 = q.z + 1.5*r1 - 0.5; float d1 = sqrt((r1-r)*(r1-r)+h1*h1) - r1;
    float h2 = q.z + 1.5*r2 - 0.5; float d2 = sqrt((r2-r)*(r2-r)+h2*h2) - r2;
    
    float d, dx, dy;
    if( d1<d2 ) { d = d1; dx=r1-r; dy=h1; }
    else        { d = d2; dx=r2-r; dy=h2; }
    
    vec3 s = q;
    q = q - vec3(0.34,-0.1,0.03);
    q.xy = mat2(0.8,0.6,-0.6,0.8)*q.xy;
    d = smin( d, torus( q, vec2(0.28,0.05) ), 0.06);
    d = smax( d, -sdEllipsoid(q,vec3(0.0,0.0,0.0),vec3(0.24,0.36,0.24) ), 0.03 );
    d = smax( d, -sdEllipsoid(s,vec3(0.52,-0.0,0.0),vec3(0.42,0.23,0.5) ), 0.05 );
    
    return d/sc;
}




// Version that just returns the t value, for surface normal computation
float sceneMap3D(vec3 pos)
{
    float t = coneTwist3(pos - vec3(4., 0., 4.));
    t = min(t, coneTwist4(pos - vec3(0., 0., 4.)));
    //float t = spiral2(pos);
    //float t = torusTwist(pos - vec3(3., 0., 3.));
    t = min(t, sphere(pos - vec3(0.1, 0., 0.1), 0.5, vec3(4.5, 0.89, 3.0)));
    //t = min(t, sphere(pos, 2.0, vec3(-8.0, 0.0, 4.0)));
    t = min(t, box(pos + vec3(0.0, 3.0, 0.0), vec3(50.0, 1.0, 50.0)));
    //t = min(t, sdCone(pos + vec3(-8.0, 0.0, 4.0), vec2(0.5, 0.86), 5.0));
    //t = min(t, coneTwist(pos + vec3(0.0, 3.0, 0.0)));
    return t;
}

// Version that returns t value and the ID of the object hit
void sceneMap3D(vec3 pos, out float t, out int obj)
{
    t = coneTwist3(pos - vec3(4., 0., 4.));
    //t = spiral2(pos);
    //t = torusTwist(pos - vec3(3., 0., 3.));
    float t2;
    obj = 0; // 0 is first tentacle
    if((t2 =  coneTwist4(pos - vec3(0., 0., 4.))) < t) {
        t = t2;
        obj = 1; // 1 is second tentacle/green thingy
    }
    if((t2 = sphere(pos - vec3(0.1, 0., 0.1), 0.5, vec3(4.5, 0.89, 3.0))) < t) {
        t = t2;
        obj = 2; // 2 is sphere
    }

    if((t2 = box(pos + vec3(0.0, 3.0, 0.0), vec3(50.0, 1.0, 50.0))) < t) {
        t = t2;
        obj = 3; // 3 is floor
    }
}


float fiveTapAO(vec3 p, vec3 n, float k) {
    float aoSum = 0.0;
    for(float i = 0.0; i < 5.0; ++i) {
        float coeff = 1.0 / pow(2.0, i);
        aoSum += coeff * (i * AO_DIST - sceneMap3D(p + n * i * AO_DIST));
    }
    return 1.0 - k * aoSum;
}

float softShadow(vec3 dir, vec3 origin, float min_t, float k) {
    float res = 1.0;
    float t = min_t;
    for(int i = 0; i < RAY_STEPS; ++i) {
        float m = sceneMap3D(origin + t * dir);
        if(m < 0.0001) {
            return 0.0;
        }
        res = min(res, k * m / t);
        t += m;
    }
    return res;
}

float shadow(vec3 dir, vec3 origin, float min_t) { 
    return softShadow(dir, origin, min_t, SHADOW_HARDNESS);
}

void march(vec3 origin, vec3 dir, out float t, out int hitObj)
{
    t = 0.001;
    for(int i = 0; i < RAY_STEPS; ++i)
    {
        vec3 pos = origin + t * dir;
    	float m;
        sceneMap3D(pos, m, hitObj);
        if(m < 0.01)
        {
            return;
        }
        t += m;
    }
    t = -1.0;
    hitObj = -1;
}

vec3 computeNormal(vec3 pos)
{
    vec3 epsilon = vec3(0.0, 0.001, 0.0);
    return normalize( vec3( sceneMap3D(pos + epsilon.yxx) - sceneMap3D(pos - epsilon.yxx),
                            sceneMap3D(pos + epsilon.xyx) - sceneMap3D(pos - epsilon.xyx),
                            sceneMap3D(pos + epsilon.xxy) - sceneMap3D(pos - epsilon.xxy)));
}

//250, 200, 158
vec3 skyColor(vec3 dir) {
    float t = smoothstep(0.0, 1.0, dir.y);
    t = clamp(0.0, 1.0, t + fbm(dir / 2.0, 4.0) * 0.9);
    t = fbm((dir / 8.0), 6.0); //+ abs(cos(iTime *0.2));
    t = smoothstep(0.0, 1.0, t);
    vec3 cloudGrad = mix(vec3(180.0, 109.0, 84.0) / 255.0, vec3(173.0, 193.0, 198.0) / 255.0, t);
    vec3 cloudGrad2 = mix(vec3(250.0, 200.0, 158.0) / 255.0, vec3(173.0, 193.0, 198.0) / 255.0, t);
    vec3 dawnGradient = mix(vec3(253.0, 96.0, 81.0) / 255.0, vec3(255.0, 229.0, 119.0) / 255.0, t);
    vec3 duskGradient = mix(vec3(48.0, 24.0, 96.0) / 255.0, vec3(144.0, 96.0, 144.0) / 255.0, t);
    //t = max(0.0, dot(dir, MAIN_LIGHT_DIR*abs(2.0*sin(0.3*iTime * HALF_PI / 4.0))));
    //return mix(duskGradient, dawnGradient, t);
    //vec3 col1 = mix(duskGradient, dawnGradient, t*abs(cos(0.1*iTime * HALF_PI / 4.0)));
    //vec3 col2 = mix(cloudGrad, cloudGrad2, t);
    //return mix(cloudGrad, cloudGrad2, t);
    t = max(0.0, dot(dir, MAIN_LIGHT_DIR));
    vec3 col1 = mix(duskGradient, dawnGradient, t);
    vec3 col2 = mix(cloudGrad, cloudGrad2, t);
    return mix(col2, col1, t);
    //return max(col1 * abs(sin(iTime)), col2*abs(cos(iTime*20.0)));
}



//used specular reflection code from https://www.shadertoy.com/view/4l3GDM
vec3 computeMaterial(int hitObj, vec3 p, vec3 d, vec3 n) {
    switch(hitObj) {
        case 0:
        // first tentacle
        vec3 l = normalize(MAIN_LIGHT_DIR);
        vec3 v = normalize(d);
        float f = dot(l, n) > 0.0 ? 1.0 : 0.0;
        vec3 r = reflect(l, n);
        return vec3(1.,1.,1.) * 0.83 * f * pow(max(0.0, dot(r, v)), 2.0) + vec3(0.955,0.775,0.24)*0.95;
        break;
        case 1:
        // second tentacle
        //return vec3(0.67, 1.0, 0.67);
        vec3 l1 = normalize(MAIN_LIGHT_DIR);
        vec3 v1 = normalize(d);
        float f1 = dot(l1, n) > 0.0 ? 1.0 : 0.0;
        vec3 r1 = reflect(l1, n);
        return vec3(1.,1.,1.) * 0.95 * f1 * pow(max(0.0, dot(r1, v1)), 8.0) + vec3(0.67, 1.0, 0.67)*0.95;
        break;
        case 2:
        vec3 l2 = normalize(MAIN_LIGHT_DIR);
        vec3 v2 = normalize(d);
        float f2 = dot(l2, n) > 0.0 ? 1.0 : 0.0;
        vec3 r2 = reflect(l2, n);
        return vec3(1.,1.,1.) * 0.95 * f2 * pow(max(0.0, dot(r2, v2)), 84.0) + vec3(1.,1.,1.)*0.75;
        break;
        case 3:
        // Floor
        //float t = floor(mod((sin(p.x) + sin(p.z)) * 0.5, 2.0));
        //return mix(vec3(0.7, 0.4, 0.2), vec3(1.0), t);
        vec3 l3 = normalize(MAIN_LIGHT_DIR);
        vec3 v3 = normalize(d);
        float f3 = dot(l3, n) > 0.0 ? 1.0 : 0.0;
        vec3 r3 = reflect(l3, n);
        return vec3(1.,1.,1.) * 0.95 * f3 * pow(max(0.0, dot(r3, v3)), 4.0) + vec3(1.,1.,1.)*0.75;
        break;
        case -1:
        // Background
        break;
    }
    return vec3(1.0);
}

Intersection sdf3D(vec3 dir, vec3 eye)
{
    float t;
    int hitObj;
    march(eye, dir, t, hitObj);
    
    if(t == -1.0) {
        return Intersection(t, skyColor(dir), vec3(eye + 1000.0 * dir), -1);
    }

    vec3 isect = eye + t * dir;
    vec3 nor = computeNormal(isect);
    
    vec3 material = computeMaterial(hitObj, isect, dir, nor);

    float fresnel = 1.0 - max(0.0, dot(normalize(eye - isect), nor));
    fresnel = 0.75 + 0.75 * fresnel;
    
    vec3 lightDir = rotateY(normalize(MAIN_LIGHT_DIR), sin(iTime * 0.5));
    
    vec3 sdfColor = mix(material, skyColor(reflect(dir, nor)) * material, fresnel) * shadow(lightDir, isect, 0.1);
    
    
    return Intersection(t, sdfColor, isect, hitObj);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Normalized pixel coordinates (from 0 to 1)
    vec2 uv = fragCoord/iResolution.xy;
    // [-1, 1]
    vec2 uv2 = 2.0 * uv - vec2(1.0);
    
    // Stratified sampled 4x4 anti-aliasing
    Intersection aaIsects[16];
    vec3 dir, eye, ref;
    int idx = 0;
    for(float i = 0.0; i < 1.0; i += 0.25) {
        for(float j = 0.0; j < 1.0; j += 0.25) {
            raycast(uv2 + (vec2(i, j) + random2(vec2(i, j))) / iResolution.xy, dir, eye, ref);
            aaIsects[idx++] = sdf3D(dir, eye);
        }
    }
    vec3 avgColor = vec3(0.0);
    for(int i = 0; i < 16; ++i) {
        avgColor += aaIsects[i].color;
    }
    avgColor *= 0.0625;
    fragColor = vec4(avgColor, 1.0);
}