#define HW_PERFORMANCE 0
// The MIT License
// Copyright © 2013 Inigo Quilez
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
//
// The license is here only not because I want to (can one
// license pieces of math?), but because people get upset
// if I don't add one...

// A list of useful distance function to simple primitives. All
// these functions (except for ellipsoid) return an exact
// euclidean distance, meaning they produce a better SDF than
// what you'd get if you were constructing them from boolean
// operations (such as cutting an infinite cylinder with two planes).

// List of other 3D SDFs:
//    https://www.shadertoy.com/playlist/43cXRl
// and
//    https://iquilezles.org/articles/distfunctions

#if HW_PERFORMANCE==0
#define AA 1
#else
#define AA 1   // make this 2 or 3 for antialiasing
#endif

//------------------------------------------------------------------
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 sdPlane( vec3 p )
{
	return p.y;
}

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

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 sdBoxFrame( vec3 p, vec3 b, float e )
{
       p = abs(p  )-b;
  vec3 q = abs(p+e)-e;

  return min(min(
      length(max(vec3(p.x,q.y,q.z),0.0))+min(max(p.x,max(q.y,q.z)),0.0),
      length(max(vec3(q.x,p.y,q.z),0.0))+min(max(q.x,max(p.y,q.z)),0.0)),
      length(max(vec3(q.x,q.y,p.z),0.0))+min(max(q.x,max(q.y,p.z)),0.0));
}
float sdEllipsoid( in vec3 p, in vec3 r ) // approximated
{
    float k0 = length(p/r);
    float k1 = length(p/(r*r));
    return k0*(k0-1.0)/k1;
}

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

float sdCappedTorus(in vec3 p, in vec2 sc, in float ra, in float rb)
{
    p.x = abs(p.x);
    float k = (sc.y*p.x>sc.x*p.y) ? dot(p.xy,sc) : length(p.xy);
    return sqrt( dot(p,p) + ra*ra - 2.0*ra*k ) - rb;
}

float sdHexPrism( vec3 p, vec2 h )
{
    vec3 q = abs(p);

    const vec3 k = vec3(-0.8660254, 0.5, 0.57735);
    p = abs(p);
    p.xy -= 2.0*min(dot(k.xy, p.xy), 0.0)*k.xy;
    vec2 d = vec2(
       length(p.xy - vec2(clamp(p.x, -k.z*h.x, k.z*h.x), h.x))*sign(p.y - h.x),
       p.z-h.y );
    return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

float sdOctogonPrism( in vec3 p, in float r, float h )
{
  const vec3 k = vec3(-0.9238795325,   // sqrt(2+sqrt(2))/2 
                       0.3826834323,   // sqrt(2-sqrt(2))/2
                       0.4142135623 ); // sqrt(2)-1 
  // reflections
  p = abs(p);
  p.xy -= 2.0*min(dot(vec2( k.x,k.y),p.xy),0.0)*vec2( k.x,k.y);
  p.xy -= 2.0*min(dot(vec2(-k.x,k.y),p.xy),0.0)*vec2(-k.x,k.y);
  // polygon side
  p.xy -= vec2(clamp(p.x, -k.z*r, k.z*r), r);
  vec2 d = vec2( length(p.xy)*sign(p.y), p.z-h );
  return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

float sdCapsule( vec3 p, vec3 a, vec3 b, float r )
{
	vec3 pa = p-a, ba = b-a;
	float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
	return length( pa - ba*h ) - r;
}

float sdRoundCone( in vec3 p, in float r1, float r2, float h )
{
    vec2 q = vec2( length(p.xz), p.y );
    
    float b = (r1-r2)/h;
    float a = sqrt(1.0-b*b);
    float k = dot(q,vec2(-b,a));
    
    if( k < 0.0 ) return length(q) - r1;
    if( k > a*h ) return length(q-vec2(0.0,h)) - r2;
        
    return dot(q, vec2(a,b) ) - r1;
}

float sdRoundCone(vec3 p, vec3 a, vec3 b, float r1, float r2)
{
    // sampling independent computations (only depend on shape)
    vec3  ba = b - a;
    float l2 = dot(ba,ba);
    float rr = r1 - r2;
    float a2 = l2 - rr*rr;
    float il2 = 1.0/l2;
    
    // sampling dependant computations
    vec3 pa = p - a;
    float y = dot(pa,ba);
    float z = y - l2;
    float x2 = dot2( pa*l2 - ba*y );
    float y2 = y*y*l2;
    float z2 = z*z*l2;

    // single square root!
    float k = sign(rr)*rr*rr*x2;
    if( sign(z)*a2*z2 > k ) return  sqrt(x2 + z2)        *il2 - r2;
    if( sign(y)*a2*y2 < k ) return  sqrt(x2 + y2)        *il2 - r1;
                            return (sqrt(x2*a2*il2)+y*rr)*il2 - r1;
}

float sdTriPrism( vec3 p, vec2 h )
{
    const float k = sqrt(3.0);
    h.x *= 0.5*k;
    p.xy /= h.x;
    p.x = abs(p.x) - 1.0;
    p.y = p.y + 1.0/k;
    if( p.x+k*p.y>0.0 ) p.xy=vec2(p.x-k*p.y,-k*p.x-p.y)/2.0;
    p.x -= clamp( p.x, -2.0, 0.0 );
    float d1 = length(p.xy)*sign(-p.y)*h.x;
    float d2 = abs(p.z)-h.y;
    return length(max(vec2(d1,d2),0.0)) + min(max(d1,d2), 0.);
}

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

// arbitrary orientation
float sdCylinder(vec3 p, vec3 a, vec3 b, float r)
{
    vec3 pa = p - a;
    vec3 ba = b - a;
    float baba = dot(ba,ba);
    float paba = dot(pa,ba);

    float x = length(pa*baba-ba*paba) - r*baba;
    float y = abs(paba-baba*0.5)-baba*0.5;
    float x2 = x*x;
    float y2 = y*y*baba;
    float d = (max(x,y)<0.0)?-min(x2,y2):(((x>0.0)?x2:0.0)+((y>0.0)?y2:0.0));
    return sign(d)*sqrt(abs(d))/baba;
}

// vertical
float sdCone( in vec3 p, in vec2 c, float h )
{
    vec2 q = h*vec2(c.x,-c.y)/c.y;
    vec2 w = vec2( length(p.xz), p.y );
    
	vec2 a = w - q*clamp( dot(w,q)/dot(q,q), 0.0, 1.0 );
    vec2 b = w - q*vec2( clamp( w.x/q.x, 0.0, 1.0 ), 1.0 );
    float k = sign( q.y );
    float d = min(dot( a, a ),dot(b, b));
    float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y)  );
	return sqrt(d)*sign(s);
}

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 sdCappedCone(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) );
}

// c is the sin/cos of the desired cone angle
float sdSolidAngle(vec3 pos, vec2 c, float ra)
{
    vec2 p = vec2( length(pos.xz), pos.y );
    float l = length(p) - ra;
	float m = length(p - c*clamp(dot(p,c),0.0,ra) );
    return max(l,m*sign(c.y*p.x-c.x*p.y));
}

float sdOctahedron(vec3 p, float s)
{
    p = abs(p);
    float m = p.x + p.y + p.z - s;

    // exact distance
    #if 0
    vec3 o = min(3.0*p - m, 0.0);
    o = max(6.0*p - m*2.0 - o*3.0 + (o.x+o.y+o.z), 0.0);
    return length(p - s*o/(o.x+o.y+o.z));
    #endif
    
    // exact distance
    #if 1
 	vec3 q;
         if( 3.0*p.x < m ) q = p.xyz;
    else if( 3.0*p.y < m ) q = p.yzx;
    else if( 3.0*p.z < m ) q = p.zxy;
    else return m*0.57735027;
    float k = clamp(0.5*(q.z-q.y+s),0.0,s); 
    return length(vec3(q.x,q.y-s+k,q.z-k)); 
    #endif
    
    // bound, not exact
    #if 0
	return m*0.57735027;
    #endif
}

float sdPyramid( in vec3 p, in float h )
{
    float m2 = h*h + 0.25;
    
    // symmetry
    p.xz = abs(p.xz);
    p.xz = (p.z>p.x) ? p.zx : p.xz;
    p.xz -= 0.5;
	
    // project into face plane (2D)
    vec3 q = vec3( p.z, h*p.y - 0.5*p.x, h*p.x + 0.5*p.y);
   
    float s = max(-q.x,0.0);
    float t = clamp( (q.y-0.5*p.z)/(m2+0.25), 0.0, 1.0 );
    
    float a = m2*(q.x+s)*(q.x+s) + q.y*q.y;
	float b = m2*(q.x+0.5*t)*(q.x+0.5*t) + (q.y-m2*t)*(q.y-m2*t);
    
    float d2 = min(q.y,-q.x*m2-q.y*0.5) > 0.0 ? 0.0 : min(a,b);
    
    // recover 3D and scale, and add sign
    return sqrt( (d2+q.z*q.z)/m2 ) * sign(max(q.z,-p.y));;
}

// la,lb=semi axis, h=height, ra=corner
float sdRhombus(vec3 p, float la, float lb, float h, float ra)
{
    p = abs(p);
    vec2 b = vec2(la,lb);
    float f = clamp( (ndot(b,b-2.0*p.xz))/dot(b,b), -1.0, 1.0 );
	vec2 q = vec2(length(p.xz-0.5*b*vec2(1.0-f,1.0+f))*sign(p.x*b.y+p.z*b.x-b.x*b.y)-ra, p.y-h);
    return min(max(q.x,q.y),0.0) + length(max(q,0.0));
}

float sdHorseshoe( in vec3 p, in vec2 c, in float r, in float le, vec2 w )
{
    p.x = abs(p.x);
    float l = length(p.xy);
    p.xy = mat2(-c.x, c.y, 
              c.y, c.x)*p.xy;
    p.xy = vec2((p.y>0.0 || p.x>0.0)?p.x:l*sign(-c.x),
                (p.x>0.0)?p.y:l );
    p.xy = vec2(p.x,abs(p.y-r))-vec2(le,0.0);
    
    vec2 q = vec2(length(max(p.xy,0.0)) + min(0.0,max(p.x,p.y)),p.z);
    vec2 d = abs(q) - w;
    return min(max(d.x,d.y),0.0) + length(max(d,0.0));
}

float sdU( in vec3 p, in float r, in float le, vec2 w )
{
    p.x = (p.y>0.0) ? abs(p.x) : length(p.xy);
    p.x = abs(p.x-r);
    p.y = p.y - le;
    float k = max(p.x,p.y);
    vec2 q = vec2( (k<0.0) ? -k : length(max(p.xy,0.0)), abs(p.z) ) - w;
    return length(max(q,0.0)) + min(max(q.x,q.y),0.0);
}

//------------------------------------------------------------------

vec2 opU( vec2 d1, vec2 d2 )
{
	return (d1.x<d2.x) ? d1 : d2;
}

//------------------------------------------------------------------

#define ZERO (min(iFrame,0))

//------------------------------------------------------------------

vec2 map( in vec3 pos )
{
    vec2 res = vec2( pos.y, 0.0 );

    // bounding box
    if( sdBox( pos-vec3(-2.0,0.3,0.25),vec3(0.3,0.3,1.0) )<res.x )
    {
      res = opU( res, vec2( sdSphere(    pos-vec3(-2.0,0.25, 0.0), 0.25 ), 26.9 ) );
	  res = opU( res, vec2( sdRhombus(  (pos-vec3(-2.0,0.25, 1.0)).xzy, 0.15, 0.25, 0.04, 0.08 ),17.0 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(0.0,0.3,-1.0),vec3(0.35,0.3,2.5) )<res.x )
    {
	res = opU( res, vec2( sdCappedTorus((pos-vec3( 0.0,0.30, 1.0))*vec3(1,-1,1), vec2(0.866025,-0.5), 0.25, 0.05), 25.0) );
    res = opU( res, vec2( sdBoxFrame(    pos-vec3( 0.0,0.25, 0.0), vec3(0.3,0.25,0.2), 0.025 ), 16.9 ) );
	res = opU( res, vec2( sdCone(        pos-vec3( 0.0,0.45,-1.0), vec2(0.6,0.8),0.45 ), 55.0 ) );
    res = opU( res, vec2( sdCappedCone(  pos-vec3( 0.0,0.25,-2.0), 0.25, 0.25, 0.1 ), 13.67 ) );
    res = opU( res, vec2( sdSolidAngle(  pos-vec3( 0.0,0.00,-3.0), vec2(3,4)/5.0, 0.4 ), 49.13 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(1.0,0.3,-1.0),vec3(0.35,0.3,2.5) )<res.x )
    {
	res = opU( res, vec2( sdTorus(      (pos-vec3( 1.0,0.30, 1.0)).xzy, vec2(0.25,0.05) ), 7.1 ) );
    res = opU( res, vec2( sdBox(         pos-vec3( 1.0,0.25, 0.0), vec3(0.3,0.25,0.1) ), 3.0 ) );
    res = opU( res, vec2( sdCapsule(     pos-vec3( 1.0,0.00,-1.0),vec3(-0.1,0.1,-0.1), vec3(0.2,0.4,0.2), 0.1  ), 31.9 ) );
	res = opU( res, vec2( sdCylinder(    pos-vec3( 1.0,0.25,-2.0), vec2(0.15,0.25) ), 8.0 ) );
    res = opU( res, vec2( sdHexPrism(    pos-vec3( 1.0,0.2,-3.0), vec2(0.2,0.05) ), 18.4 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(-1.0,0.35,-1.0),vec3(0.35,0.35,2.5))<res.x )
    {
	res = opU( res, vec2( sdPyramid(    pos-vec3(-1.0,-0.6,-3.0), 1.0 ), 13.56 ) );
	res = opU( res, vec2( sdOctahedron( pos-vec3(-1.0,0.15,-2.0), 0.35 ), 23.56 ) );
    res = opU( res, vec2( sdTriPrism(   pos-vec3(-1.0,0.15,-1.0), vec2(0.3,0.05) ),43.5 ) );
    res = opU( res, vec2( sdEllipsoid(  pos-vec3(-1.0,0.25, 0.0), vec3(0.2, 0.25, 0.05) ), 43.17 ) );
    res = opU( res, vec2( sdHorseshoe(  pos-vec3(-1.0,0.25, 1.0), vec2(cos(1.3),sin(1.3)), 0.2, 0.3, vec2(0.03,0.08) ), 11.5 ) );
    }

    // bounding box
    if( sdBox( pos-vec3(2.0,0.3,-1.0),vec3(0.35,0.3,2.5) )<res.x )
    {
    res = opU( res, vec2( sdOctogonPrism(pos-vec3( 2.0,0.2,-3.0), 0.2, 0.05), 51.8 ) );
    res = opU( res, vec2( sdCylinder(    pos-vec3( 2.0,0.14,-2.0), vec3(0.1,-0.1,0.0), vec3(-0.2,0.35,0.1), 0.08), 31.2 ) );
	res = opU( res, vec2( sdCappedCone(  pos-vec3( 2.0,0.09,-1.0), vec3(0.1,0.0,0.0), vec3(-0.2,0.40,0.1), 0.15, 0.05), 46.1 ) );
    res = opU( res, vec2( sdRoundCone(   pos-vec3( 2.0,0.15, 0.0), vec3(0.1,0.0,0.0), vec3(-0.1,0.35,0.1), 0.15, 0.05), 51.7 ) );
    res = opU( res, vec2( sdRoundCone(   pos-vec3( 2.0,0.20, 1.0), 0.2, 0.1, 0.3 ), 37.0 ) );
    }
    
    return res;
}

// https://iquilezles.org/articles/boxfunctions
vec2 iBox( in vec3 ro, in vec3 rd, in vec3 rad ) 
{
    vec3 m = 1.0/rd;
    vec3 n = m*ro;
    vec3 k = abs(m)*rad;
    vec3 t1 = -n - k;
    vec3 t2 = -n + k;
	return vec2( max( max( t1.x, t1.y ), t1.z ),
	             min( min( t2.x, t2.y ), t2.z ) );
}

vec2 raycast( in vec3 ro, in vec3 rd )
{
    vec2 res = vec2(-1.0,-1.0);

    float tmin = 1.0;
    float tmax = 20.0;

    // raytrace floor plane
    float tp1 = (0.0-ro.y)/rd.y;
    if( tp1>0.0 )
    {
        tmax = min( tmax, tp1 );
        res = vec2( tp1, 1.0 );
    }
    //else return res;
    
    // raymarch primitives   
    vec2 tb = iBox( ro-vec3(0.0,0.4,-0.5), rd, vec3(2.5,0.41,3.0) );
    if( tb.x<tb.y && tb.y>0.0 && tb.x<tmax)
    {
        //return vec2(tb.x,2.0);
        tmin = max(tb.x,tmin);
        tmax = min(tb.y,tmax);

        float t = tmin;
        for( int i=0; i<70 && t<tmax; i++ )
        {
            vec2 h = map( ro+rd*t );
            if( abs(h.x)<(0.0001*t) )
            { 
                res = vec2(t,h.y); 
                break;
            }
            t += h.x;
        }
    }
    
    return res;
}

// https://iquilezles.org/articles/rmshadows
float calcSoftshadow( in vec3 ro, in vec3 rd, in float mint, in float tmax )
{
    // bounding volume
    float tp = (0.8-ro.y)/rd.y; if( tp>0.0 ) tmax = min( tmax, tp );

    float res = 1.0;
    float t = mint;
    for( int i=ZERO; i<24; i++ )
    {
		float h = map( ro + rd*t ).x;
        float s = clamp(8.0*h/t,0.0,1.0);
        res = min( res, s );
        t += clamp( h, 0.01, 0.2 );
        if( res<0.004 || t>tmax ) break;
    }
    res = clamp( res, 0.0, 1.0 );
    return res*res*(3.0-2.0*res);
}

// https://iquilezles.org/articles/normalsSDF
vec3 calcNormal( in vec3 pos )
{
#if 0
    vec2 e = vec2(1.0,-1.0)*0.5773*0.0005;
    return normalize( e.xyy*map( pos + e.xyy ).x + 
					  e.yyx*map( pos + e.yyx ).x + 
					  e.yxy*map( pos + e.yxy ).x + 
					  e.xxx*map( pos + e.xxx ).x );
#else
    // inspired by tdhooper and klems - a way 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).x;
      //if( n.x+n.y+n.z>100.0 ) break;
    }
    return normalize(n);
#endif    
}

// https://iquilezles.org/articles/nvscene2008/rwwtt.pdf
float calcAO( in vec3 pos, in vec3 nor )
{
	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 ).x;
        occ += (h-d)*sca;
        sca *= 0.95;
        if( occ>0.35 ) break;
    }
    return clamp( 1.0 - 3.0*occ, 0.0, 1.0 ) * (0.5+0.5*nor.y);
}

// https://iquilezles.org/articles/checkerfiltering
float checkersGradBox( in vec2 p, in vec2 dpdx, in vec2 dpdy )
{
    // filter kernel
    vec2 w = abs(dpdx)+abs(dpdy) + 0.001;
    // analytical integral (box filter)
    vec2 i = 2.0*(abs(fract((p-0.5*w)*0.5)-0.5)-abs(fract((p+0.5*w)*0.5)-0.5))/w;
    // xor pattern
    return 0.5 - 0.5*i.x*i.y;                  
}

vec3 render( in vec3 ro, in vec3 rd, in vec3 rdx, in vec3 rdy )
{ 
    // background
    vec3 col = vec3(0.7, 0.7, 0.9) - max(rd.y,0.0)*0.3;
    
    // raycast scene
    vec2 res = raycast(ro,rd);
    float t = res.x;
	float m = res.y;
    if( m>-0.5 )
    {
        vec3 pos = ro + t*rd;
        vec3 nor = (m<1.5) ? vec3(0.0,1.0,0.0) : calcNormal( pos );
        vec3 ref = reflect( rd, nor );
        
        // material        
        col = 0.2 + 0.2*sin( m*2.0 + vec3(0.0,1.0,2.0) );
        float ks = 1.0;
        
        if( m<1.5 )
        {
            // project pixel footprint into the plane
            vec3 dpdx = ro.y*(rd/rd.y-rdx/rdx.y);
            vec3 dpdy = ro.y*(rd/rd.y-rdy/rdy.y);

            float f = checkersGradBox( 3.0*pos.xz, 3.0*dpdx.xz, 3.0*dpdy.xz );
            col = 0.15 + f*vec3(0.05);
            ks = 0.4;
        }

        // lighting
        float occ = calcAO( pos, nor );
        
		vec3 lin = vec3(0.0);

        // sun
        {
            vec3  lig = normalize( vec3(-0.5, 0.4, -0.6) );
            vec3  hal = normalize( lig-rd );
            float dif = clamp( dot( nor, lig ), 0.0, 1.0 );
          //if( dif>0.0001 )
        	      dif *= calcSoftshadow( pos, lig, 0.02, 2.5 );
			float spe = pow( clamp( dot( nor, hal ), 0.0, 1.0 ),16.0);
                  spe *= dif;
                  spe *= 0.04+0.96*pow(clamp(1.0-dot(hal,lig),0.0,1.0),5.0);
                //spe *= 0.04+0.96*pow(clamp(1.0-sqrt(0.5*(1.0-dot(rd,lig))),0.0,1.0),5.0);
            lin += col*2.20*dif*vec3(1.30,1.00,0.70);
            lin +=     5.00*spe*vec3(1.30,1.00,0.70)*ks;
        }
        // sky
        {
            float dif = sqrt(clamp( 0.5+0.5*nor.y, 0.0, 1.0 ));
                  dif *= occ;
            float spe = smoothstep( -0.2, 0.2, ref.y );
                  spe *= dif;
                  spe *= 0.04+0.96*pow(clamp(1.0+dot(nor,rd),0.0,1.0), 5.0 );
          //if( spe>0.001 )
                  spe *= calcSoftshadow( pos, ref, 0.02, 2.5 );
            lin += col*0.60*dif*vec3(0.40,0.60,1.15);
            lin +=     2.00*spe*vec3(0.40,0.60,1.30)*ks;
        }
        // back
        {
        	float dif = clamp( dot( nor, normalize(vec3(0.5,0.0,0.6))), 0.0, 1.0 )*clamp( 1.0-pos.y,0.0,1.0);
                  dif *= occ;
        	lin += col*0.55*dif*vec3(0.25,0.25,0.25);
        }
        // sss
        {
            float dif = pow(clamp(1.0+dot(nor,rd),0.0,1.0),2.0);
                  dif *= occ;
        	lin += col*0.25*dif*vec3(1.00,1.00,1.00);
        }
        
		col = lin;

        col = mix( col, vec3(0.7,0.7,0.9), 1.0-exp( -0.0001*t*t*t ) );
    }

	return vec3( clamp(col,0.0,1.0) );
}

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 )
{
    vec2 mo = iMouse.xy/iResolution.xy;
	float time = 32.0 + iTime*1.5;

    // camera	
    vec3 ta = vec3( 0.25, -0.75, -0.25 );
    vec3 ro = ta + vec3( 6.0*cos(0.1*time + 7.0*mo.x), 4.5, 6.0*sin(0.1*time + 7.0*mo.x) );
    // camera-to-world transformation
    mat3 ca = setCamera( ro, ta, 0.0 );

    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;
#else    
        vec2 p = (2.0*fragCoord-iResolution.xy)/iResolution.y;
#endif

        // focal length
        const float fl = 2.5;
        
        // ray direction
        vec3 rd = ca * normalize( vec3(p,fl) );

         // ray differentials
        vec2 px = (2.0*(fragCoord+vec2(1.0,0.0))-iResolution.xy)/iResolution.y;
        vec2 py = (2.0*(fragCoord+vec2(0.0,1.0))-iResolution.xy)/iResolution.y;
        vec3 rdx = ca * normalize( vec3(px,fl) );
        vec3 rdy = ca * normalize( vec3(py,fl) );
        
        // render	
        vec3 col = render( ro, rd, rdx, rdy );

        // gain
        // col = col*3.0/(2.5+col);
        
		// gamma
        col = pow( col, vec3(0.4545) );

        tot += col;
#if AA>1
    }
    tot /= float(AA*AA);
#endif
    
    fragColor = vec4( tot, 1.0 );
}

#define AS_LIB 1

int get_shape() {
    return int(iTime) % 52;
}

#include "/forest/shader/implicit.glsl"

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = 2.*(fragCoord-iResolution.xy/2.)/iResolution.y; // contains [-1,1]^2
    vec3 col = vec3(0.);
    
    // Camera rays
    vec3 camPos = vec3(4.,0.,0.);
    vec3 camDir = - normalize(camPos);
    vec3 rayPos, rayDir;
    float zoom = 1.3; // 1.8*cos(iTime);
    // if (checkKey(KEY_E)) zoom = 0.5;
    float fov = 0.4*zoom;
    float fov_ortho = 1.5*zoom;
    #if perspective
        // perspective cam
        rayPos = camPos;
        rayDir = normalize(camDir + fov*vec3(0., uv.x, uv.y));
    #else
        // orthographic cam
        rayPos = camPos + fov_ortho*vec3(0., uv.x, uv.y);
        rayDir = camDir;
    #endif
    // for perspective background in orthographic mode
    vec3 cubemapDir = normalize(camDir + fov*vec3(0., uv.x, uv.y));
    
    // Mouse-controlled rotation
    vec2 mouse = initMouse + vec2(0.015625*sin(iTime*PI), 0.0); // initMouse; // iMouse.xy == vec2(0.,0.) ? initMouse : (iMouse.xy/iResolution.xy - 0.5);
    float yaw = clamp(- mouse.x * 2.*PI * 1., -PI,PI);
    float pitch = clamp( mouse.y * PI * 1.2, -PI*0.5, PI*0.5);
    // pitch and yaw rotations (column-wise matrices)
    mat3 rot = mat3(cos(yaw), sin(yaw), 0., -sin(yaw), cos(yaw), 0., 0., 0., 1.);
    rot = rot * mat3(cos(pitch), 0., -sin(pitch), 0., 1., 0., sin(pitch), 0., cos(pitch));
    // apply
    camPos = rot*camPos;
    camDir = rot*camDir;
    rayPos = rot*rayPos;
    rayDir = rot*rayDir;
    cubemapDir = rot*cubemapDir;
    
    //cubemapDir = vec3(cubemapDir.x, cubemapDir.z, cubemapDir.y);
    
    vec3 hitPoint = raycast(rayPos, rayDir);
    
    if (hitPoint == BINGO) { fragColor = vec4(BINGO,1.0); return; }
    //if (hitPoint == NOHIT) { fragColor = vec4(NOHIT,1.0); return; }
    //if (hitPoint == NOBOUNDHIT) { fragColor = vec4(NOBOUNDHIT,1.0); return; }
    //if (hitPoint == ESCAPEDBOUNDS) { fragColor = vec4(ESCAPEDBOUNDS,1.0); return; }
    //if (hitPoint == MAXDISTREACHED) { fragColor = vec4(MAXDISTREACHED,1.0); return; }
    //if (hitPoint == MAXITERREACHED) { fragColor = vec4(MAXITERREACHED,1.0); return; }
    
    if (hitPoint == NOBOUNDHIT  || hitPoint == NOHIT || hitPoint == ESCAPEDBOUNDS || hitPoint == MAXITERREACHED) {
        //fragColor = vec4(vec3(0.2),1.0); return;

        // make background transparent
        fragColor = vec4(0.0,0.0,0.0,0.0); return;
        
        col = with_background(cubemapDir);
        
        #if showBoundingCube
        // darken bounding cube
        if (hitPoint != NOBOUNDHIT) { col *= vec3(0.7); }
        #endif

        fragColor = vec4(col,1.0); return;
    }

    vec3 grad = gradf(hitPoint+1.1*EPS*(-rayDir));
    float s = -sign(dot(grad,rayDir));
    
    col = with_color_mode(grad, s, hitPoint, camPos);

    col = clamp(col, 0., 1.);

    col = with_surface_pattern(col, hitPoint);

    col = with_shading(col, grad, s, rayDir);
    
    col = clamp(col, 0., 1.);

    fragColor = vec4(col,1.0);
}

//********************************************************************
//
// Begin of automatic differentiation header
// Full code with additional functions (gradients, jacobians, ...) can be found
// at: https://github.com/sibaku/glsl-autodiff
//
//********************************************************************
#ifndef HESSNUM_3_H_
#define HESSNUM_3_H_

// This file contains methods to compute the gradient and hessian 
// of a scalar valued 3 dimensional function using automatic forward differentiation

//--------------------------------
// Types
//--------------------------------

// Data type to hold information about a scalar valued 3 dimensional function
// These should be created by the constH3 (for constants) and varH3 (for variables) helpers
struct HNum3
{
    // The current value
    float val;
    // The current gradient
    vec3 g;
    // The current hessian
    mat3 h;
};

//--------------------------------
// Prototypes
//--------------------------------

/**
* Creates a constant HNum3
* @param val The current value of the constant
*/
HNum3 constH3(in float val);
/**
* Creates a HNum3 corresponding to the variable with the given index
* @param val The current value of the variable
* @param index The variable's index
*/
HNum3 varH3(in float val, in int index);
/**
* Creates a HNum3 corresponding to the variable x (index = 0)
* @param val The current value of the variable
*/
HNum3 varH3x(in float val);
/**
* Creates a HNum3 corresponding to the variable y (index = 1)
* @param val The current value of the variable
*/
HNum3 varH3y(in float val);
/**
* Creates a HNum3 corresponding to the variable z (index = 2)
* @param val The current value of the variable
*/
HNum3 varH3z(in float val);
HNum3 add(in HNum3 a, in HNum3 b);
HNum3 add(in HNum3 a, in float b);
HNum3 add(in float a, in HNum3 b);
HNum3 sub(in HNum3 a, in HNum3 b);
HNum3 sub(in HNum3 a, in float b);
HNum3 sub(in float a, in HNum3 b);
HNum3 mult(in HNum3 a, in HNum3 b);
HNum3 mult(in HNum3 a, in float b);
HNum3 mult(in float a, in HNum3 b);
HNum3 neg(in HNum3 a);
HNum3 div(in HNum3 a, in HNum3 b);
HNum3 div(in HNum3 a, in float b);
HNum3 div(in float a, in HNum3 b);
HNum3 inv(in HNum3 a);
HNum3 a_pow(in HNum3 a, in HNum3 b);
HNum3 a_pow(in HNum3 a, in float b);
HNum3 a_pow(in float a, in HNum3 b);
HNum3 a_min(in HNum3 a, in HNum3 b);
HNum3 a_max(in HNum3 a, in HNum3 b);
HNum3 a_exp2(in HNum3 a);
HNum3 a_inversesqrt(in HNum3 a);
HNum3 a_atan(in HNum3 a);
HNum3 a_sqrt(in HNum3 a);
HNum3 a_sinh(in HNum3 a);
HNum3 a_ceil(in HNum3 a);
HNum3 a_tan(in HNum3 a);
HNum3 a_asinh(in HNum3 a);
HNum3 a_asin(in HNum3 a);
HNum3 a_acosh(in HNum3 a);
HNum3 a_abs(in HNum3 a);
HNum3 a_exp(in HNum3 a);
HNum3 a_cosh(in HNum3 a);
HNum3 a_floor(in HNum3 a);
HNum3 a_log(in HNum3 a);
HNum3 a_atanh(in HNum3 a);
HNum3 a_log2(in HNum3 a);
HNum3 a_acos(in HNum3 a);
HNum3 a_tanh(in HNum3 a);
HNum3 a_cos(in HNum3 a);
HNum3 a_sin(in HNum3 a);
HNum3 a_atan2(in HNum3 y, in HNum3 x);
HNum3 a_atan2(in HNum3 y, in float x);
HNum3 a_atan2(in float y, in HNum3 x);
HNum3 a_mix(in HNum3 a, in HNum3 b, in HNum3 t);
HNum3 a_mix(in HNum3 a, in HNum3 b, in float t);
HNum3 a_mix(in HNum3 a, in float b, in HNum3 t);
HNum3 a_mix(in HNum3 a, in float b, in float t);
HNum3 a_mix(in float a, in HNum3 b, in HNum3 t);
HNum3 a_mix(in float a, in HNum3 b, in float t);
HNum3 a_mix(in float a, in float b, in HNum3 t);

//--------------------------------
// Macros
//--------------------------------

#define HESSIAN3(f,x, y, z,result)  {     result = f(varH3x(x), varH3y(y), varH3z(z)); }

//--------------------------------
// Utilities prototypes
//--------------------------------

mat3 a_outerProduct(in vec3 a, in vec3 b);

//--------------------------------
// Implementation
//--------------------------------

HNum3 constH3(in float val)
{
    return HNum3(val, vec3(0.0), mat3(0.0));
}
//--------------------------------
HNum3 varH3(in float val, in int index)
{   
    vec3 g = vec3(0.0);
    g[index] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 varH3x(in float val)
{   
    vec3 g = vec3(0.0);
    g[0] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 varH3y(in float val)
{   
    vec3 g = vec3(0.0);
    g[1] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 varH3z(in float val)
{   
    vec3 g = vec3(0.0);
    g[2] = 1.0;
    return HNum3(val, g, mat3(0.0));
}
//--------------------------------
HNum3 add(in HNum3 a, in HNum3 b)
{
    return HNum3(a.val + b.val , a.g + b.g, a.h + b.h);
}
//--------------------------------
HNum3 add(in HNum3 a, in float b)
{
    return HNum3(a.val + b , a.g, a.h);
}
//--------------------------------
HNum3 add(in float a, in HNum3 b)
{
    return HNum3(a + b.val , b.g, b.h);
}
//--------------------------------
HNum3 sub(in HNum3 a, in HNum3 b)
{
    return HNum3(a.val - b.val , a.g - b.g, a.h - b.h);
}
//--------------------------------
HNum3 sub(in HNum3 a, in float b)
{
    return HNum3(a.val - b , a.g, a.h);
}
//--------------------------------
HNum3 sub(in float a, in HNum3 b)
{
    return HNum3(a - b.val , - b.g, - b.h);
}
//--------------------------------
HNum3 mult(in HNum3 a, in HNum3 b)
{
    return HNum3(a.val * b.val, 
        a.val*b.g + b.val*a.g, 
        a.val*b.h + b.val*a.h + a_outerProduct(b.g,a.g) + a_outerProduct(a.g,b.g)
    );
}
//--------------------------------
HNum3 mult(in HNum3 a, in float b)
{
    return HNum3(a.val * b, b*a.g, b*a.h);
}
//--------------------------------
HNum3 mult(in float a, in HNum3 b)
{
    return HNum3(a * b.val, a*b.g, a*b.h);
}
//--------------------------------
HNum3 neg(in HNum3 a)
{
    return mult(-1.0,a);
}
//--------------------------------
HNum3 div(in HNum3 a, in HNum3 b)
{
    float b1 = b.val;
    float b2 = b1*b1;
    float b3 = b2*b1;

    return HNum3(a.val / b.val , 
        (b.val*a.g - a.val*b.g)/b2, 
        2.0*a.val/b3*a_outerProduct(b.g,b.g) 
        - a.val/b2*b.h
        + a.h/b1 
        - a_outerProduct(b.g/b2, a.g)
        - a_outerProduct(a.g/b2, b.g)
    );
}
//--------------------------------
HNum3 div(in HNum3 a, in float b)
{
    return HNum3(a.val / b, a.g/b, a.h/b);
}
//--------------------------------
HNum3 div(in float a, in HNum3 b)
{
    float b1 = b.val;
    float b2 = b1*b1;
    float b3 = b2*b1;

    return HNum3(a / b.val, 
        -a*b.g/b2, 
        2.0*a/b3*a_outerProduct(b.g,b.g) - a/b2*b.h
    );
}
//--------------------------------
HNum3 inv(in HNum3 a)
{
    return div(1.0, a);
}
//--------------------------------
HNum3 a_pow(in HNum3 a, in HNum3 b)
{
    return a_exp(mult(b,a_log(a)));
}
//--------------------------------
HNum3 a_pow(in HNum3 a, in float b)
{
    // constant exponent -> make special case
    float v = pow(a.val, b); // value f(a(x))
    float da = b*pow(a.val,b-1.0); // first derivative f'(a(x))
    float dda = b*(b-1.0)*pow(a.val,b-2.0); // second derivative f''(a(x))
    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_pow(in float a, in HNum3 b)
{
    return a_exp(mult(b,log(a)));
}
//--------------------------------
HNum3 a_min(in HNum3 a, in HNum3 b)
{
    if(a.val < b.val)
    {
        return a;
    }
    return b;
}
//--------------------------------
HNum3 a_max(in HNum3 a, in HNum3 b)
{
    if(a.val > b.val)
    {
        return a;
    }
    return b;
}
//--------------------------------
HNum3 a_exp2(in HNum3 a)
{
    float v = exp2(a.val); // value f(a(x))
    float da = log(2.0)*exp2(a.val); // first derivative f'(a(x))
    float dda = log(2.0)*log(2.0)*exp2(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_inversesqrt(in HNum3 a)
{
    float v = inversesqrt(a.val); // value f(a(x))
    float da = -0.5/pow(sqrt(a.val),3.0); // first derivative f'(a(x))
    float dda = 0.75/pow(sqrt(a.val),5.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_atan(in HNum3 a)
{
    float v = atan(a.val); // value f(a(x))
    float da = 1.0/(1.0 + a.val * a.val); // first derivative f'(a(x))
    float dda = -2.0*a.val/pow(1.0 + a.val * a.val, 2.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_sqrt(in HNum3 a)
{
    float v = sqrt(a.val); // value f(a(x))
    float da = 0.5/sqrt(a.val); // first derivative f'(a(x))
    float dda = -0.25/pow(sqrt(a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_sinh(in HNum3 a)
{
    float v = sinh(a.val); // value f(a(x))
    float da = cosh(a.val); // first derivative f'(a(x))
    float dda = sinh(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_ceil(in HNum3 a)
{
    float v = ceil(a.val); // value f(a(x))
    float da = 0.0; // first derivative f'(a(x))
    float dda = 0.0; // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_tan(in HNum3 a)
{
    float v = tan(a.val); // value f(a(x))
    float da = 1.0 + pow(tan(a.val),2.0); // first derivative f'(a(x))
    float dda = 2.0*tan(a.val)*(1.0 + pow(tan(a.val),2.0)); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_asinh(in HNum3 a)
{
    float v = asinh(a.val); // value f(a(x))
    float da = 1.0/sqrt(1.0 + a.val * a.val); // first derivative f'(a(x))
    float dda = -a.val/pow(sqrt(1.0 + a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_asin(in HNum3 a)
{
    float v = asin(a.val); // value f(a(x))
    float da = 1.0/sqrt(1.0 - a.val * a.val); // first derivative f'(a(x))
    float dda = a.val/pow(sqrt(1.0 - a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_acosh(in HNum3 a)
{
    float v = acosh(a.val); // value f(a(x))
    float da = 1.0/sqrt(-1.0 + a.val * a.val); // first derivative f'(a(x))
    float dda = -a.val/pow(sqrt(-1.0 + a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_abs(in HNum3 a)
{
    float v = abs(a.val); // value f(a(x))
    float da = a.val < 0.0 ? -1.0 : 1.0; // first derivative f'(a(x))
    float dda = 0.0; // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_exp(in HNum3 a)
{
    float v = exp(a.val); // value f(a(x))
    float da = exp(a.val); // first derivative f'(a(x))
    float dda = exp(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_cosh(in HNum3 a)
{
    float v = cosh(a.val); // value f(a(x))
    float da = sinh(a.val); // first derivative f'(a(x))
    float dda = cosh(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_floor(in HNum3 a)
{
    float v = floor(a.val); // value f(a(x))
    float da = 0.0; // first derivative f'(a(x))
    float dda = 0.0; // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_log(in HNum3 a)
{
    float v = log(a.val); // value f(a(x))
    float da = 1.0/a.val; // first derivative f'(a(x))
    float dda = -1.0/(a.val * a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_atanh(in HNum3 a)
{
    float v = atanh(a.val); // value f(a(x))
    float da = 1.0/(1.0 - a.val * a.val); // first derivative f'(a(x))
    float dda = 2.0*a.val/pow(1.0 - a.val * a.val,2.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_log2(in HNum3 a)
{
    float v = log2(a.val); // value f(a(x))
    float da = 1.0/(a.val * log(2.0)); // first derivative f'(a(x))
    float dda = -1.0/(a.val * a.val * log(2.0)); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_acos(in HNum3 a)
{
    float v = acos(a.val); // value f(a(x))
    float da = -1.0/sqrt(1.0 - a.val * a.val); // first derivative f'(a(x))
    float dda = -a.val/pow(sqrt(1.0 - a.val * a.val),3.0); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_tanh(in HNum3 a)
{
    float v = tanh(a.val); // value f(a(x))
    float da = 1.0 - pow(tanh(a.val),2.0); // first derivative f'(a(x))
    float dda = -2.0*tanh(a.val)*(1.0 - pow(tanh(a.val),2.0)); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_cos(in HNum3 a)
{
    float v = cos(a.val); // value f(a(x))
    float da = -sin(a.val); // first derivative f'(a(x))
    float dda = -cos(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_sin(in HNum3 a)
{
    float v = sin(a.val); // value f(a(x))
    float da = cos(a.val); // first derivative f'(a(x))
    float dda = -sin(a.val); // second derivative f''(a(x))

    return HNum3(v , da * a.g,  da * a.h + dda * a_outerProduct(a.g,a.g));
}
//--------------------------------
HNum3 a_atan2(in HNum3 y, in HNum3 x)
{
    const float pi = 3.14159265; 
    // from https://en.wikipedia.org/wiki/Atan2
    if(x.val > 0.0)
    {
        HNum3 n = a_sqrt(add(mult(x,x),mult(y,y)));
        HNum3 inner = div(y, add(n,x));
        
        return mult(2.0,a_atan(inner));
        
    }else if(x.val <= 0.0 && abs(y.val) > 1E-6)
    {
        HNum3 n = a_sqrt(add(mult(x,x),mult(y,y)));
        HNum3 inner = div(sub(n,x),y);
         return mult(2.0,a_atan(inner));
    }else if(x.val < 0.0 && abs(y.val) <= 1E-6)
    {
        return constH3(pi);
    }
    // return 0 for undefined
    return constH3(0.0); 
}
//--------------------------------
HNum3 a_atan2(in HNum3 y, in float x)
{
    return a_atan2(y,constH3(x));
}
//--------------------------------
HNum3 a_atan2(in float y, in HNum3 x)
{
    return a_atan2(constH3(y),x);
}
//--------------------------------
HNum3 a_mix(in HNum3 a, in HNum3 b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in HNum3 a, in HNum3 b, in float t)
{
    return add(mult(a, 1.0 - t), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in HNum3 a, in float b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in HNum3 a, in float b, in float t)
{
    return add(mult(a, 1.0 - t), b*t);
}

//--------------------------------
HNum3 a_mix(in float a, in HNum3 b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in float a, in HNum3 b, in float t)
{
    return add(a * (1.0 - t), mult(b, t));
}

//--------------------------------
HNum3 a_mix(in float a, in float b, in HNum3 t)
{
    return add(mult(a, sub(1.0, t)), mult(b, t));
}


//--------------------------------
// Implementation prototypes
//--------------------------------

mat3 a_outerProduct(in vec3 a, in vec3 b)
{
    return mat3(a * b[0], a * b[1], a * b[2]);
}


#endif // HESSNUM_3_H_
//********************************************************************
//
// End automatic differentiation header
//
//********************************************************************

//********************************************************************
//
// Metric functions
//
//********************************************************************
#define PI 3.14159265359

// float pSphere(vec3 p, float s) {return length(p)-s;}
HNum3 fSphere(in HNum3 x, in HNum3 y, in HNum3 z,
                 float globalTime) {             
                 
   return sub(a_sqrt(add(add(mult(x,x), mult(y,y)), mult(z,z))), 5.0)  ;
}

// Predefined functions!
// Just comment in the one you want to see! (comment out the other ones)
// If none is defined, the default identity operation is used (f(x,y,z) = 0).
// You can comment out all predefined ones and implement your own in fCustom

// #define EXPONENTIAL
// #define SADDLE// 
// #define WAVES
// #define PARABOLA
// #define IDENTITY

HNum3 fExponential(in HNum3 x, in HNum3 y, in HNum3 z,
                      float globalTime);
HNum3 fSaddle(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime);
HNum3 fWaves(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime);
HNum3 fParabola(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime);
HNum3 fIdentity(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime);
                   
// Implement your custom height function here and comment out all the defines
// above! You can use all the mathematical operations defined in the automatic
// differentiation header above
HNum3 fCustom(in HNum3 x, in HNum3 y, in HNum3 z,
                 float globalTime) {
  return div(fExponential(x, y, z, globalTime), 1.0);
  // return a_max(fExponential(x, y, z, globalTime), fSphere(x, y, z, globalTime));
}

// Height function used for the metric. More information in compute_update
HNum3 f(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime) {

#if defined EXPONENTIAL
  return fExponential(x, y, z, globalTime);
#elif defined SADDLE
  return fSaddle(x, y, z, globalTime);
#elif defined WAVES
  return fWaves(x, y, z, globalTime);
#elif defined PARABOLA
  return fParabola(x, y, z, globalTime);
#elif defined IDENTITY
  return fIdentity(x, y, z, globalTime);
#else
  return fCustom(x, y, z, globalTime);
#endif
}

HNum3 hessF(vec3 p, float globalTime) {
  // compute hessian with addition time parameter

  vec3 uGrad = vec3(1., 0., 0.);
  HNum3 uHessian = HNum3(p.x, uGrad, mat3(0.));
  vec3 vGrad = vec3(0., 1., 0.);
  HNum3 vHessian = HNum3(p.y, vGrad, mat3(0.));
  vec3 wGrad = vec3(0., 0., 1.);
  HNum3 wHessian = HNum3(p.z, wGrad, mat3(0.));
  return f(uHessian, vHessian, wHessian, globalTime);
}
vec3 compute_update(vec3 p, vec3 v, float globalTime) {

  // explicit form of the christoffel symbols for the metric given by graph of
  // the function f(x_1,x_1,x_3) df/dx_i = f_i d^2f/(dx_j dx_j) = f_{i,j}
  // gamma_{i,j}^k = f_k*f_{i,j}/(1 + |grad(f)|^2)

  // This is the update given by the differential geodesic equation used for
  // numerical integration in local coordinates (x_i) for a curve c_i(t) = x_i
  // d^2 c_k/dt^2 = gamma_{i,j}^k * d c_i/dt * d c_j/dt

  // The local position's velocity is given by the ray's current direction

  // compute first and second order derivatives
  HNum3 r = hessF(p, globalTime);

  vec3 vp = vec3(0.0);

  for (int k = 0; k < 3; k++) {
    for (int i = 0; i < 3; i++) {
      for (int j = 0; j < 3; j++) {
        vp[k] += r.g[k] * r.h[j][i] * v[i] * v[j];
      }
    }
  }

  return -vp / (1.0 + dot(r.g, r.g));
}

//********************************************************************
//
// Data
//
//********************************************************************

const float deg2Rad = PI / 180.0;
const vec2 D_lastMouse = vec2(0., 0.);
vec4 D_thetaPhi;
const vec3 D_pos = vec3(0.0);

//********************************************************************
//
// Camera functions
//
//********************************************************************
// This assumes the pixel position px to be in [0,1],
// which can be done by (x+0.5)/w or (y+0.5)/h (or h-y +0.5 for screens
// with top left origin) to sample pixel centers
vec3 createRay(vec2 px, mat4 PInv, mat4 VInv);

mat4 translate(vec3 t);
mat4 translateInv(vec3 t);

mat4 scale(vec3 s);
mat4 scaleInv(vec3 s);

mat4 rightToLeft();
mat4 rightToLeftInv();

mat4 ortho(float l, float r, float b, float t, float n, float f);
mat4 orthoInv(float l, float r, float b, float t, float n, float f);

mat4 projection(float n, float f);
mat4 projectionInv(float n, float f);

mat4 perspective(float fov, float aspect, float n, float f);
mat4 perspectiveInv(float fov, float aspect, float n, float f);

mat4 lookAt(vec3 eye, vec3 center, vec3 up);
mat4 lookAtInv(vec3 eye, vec3 center, vec3 up);

//********************************************************************
//
// Implementation
//
//********************************************************************
HNum3 fExponential(in HNum3 x, in HNum3 y, in HNum3 z,
                      float globalTime) {
  HNum3 rx = sub(x, 0.0);
  HNum3 ry = sub(y, 1.0);
  HNum3 rz = sub(z, 2.0);

  HNum3 x2 = mult(rx, rx);
  HNum3 y2 = mult(ry, ry);
  HNum3 z2 = mult(rz, rz);

  HNum3 sum = neg(add(x2, add(y2, z2)));

  HNum3 f1 = a_exp(mult(sum, 1.0 / (abs(sin(globalTime / 5.0)) + 0.1)));

  rx = sub(x, -6.0);
  ry = sub(y, -0.25);
  rz = sub(z, 3.5);

  x2 = mult(rx, rx);
  y2 = mult(ry, ry);
  z2 = mult(rz, rz);

  sum = neg(add(x2, add(y2, z2)));
  sum = mult(sum, 0.1);
  HNum3 f2 = mult(a_exp(sum), abs(sin(globalTime / 5.0)) * 5.0);
  return add(f1, f2);
}
HNum3 fSaddle(in HNum3 x, in HNum3 y, in HNum3 z,
                 float globalTime) {
  return mult(0.25 * pow(sin(globalTime / 2.0), 2.0),
              sub(add(mult(x, x), mult(z, z)), mult(y, y)));
}
HNum3 fWaves(in HNum3 x, in HNum3 y, in HNum3 z, float globalTime) {
  float frequency = 2.0 * PI / (10.0 + 3.0 * sin(0.1 * globalTime));
  return a_cos(add(mult(x, frequency), mult(z, frequency)));
}
HNum3 fParabola(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime) {
  HNum3 rel = sub(x, 0.0);
  HNum3 sum = mult(rel, rel);

  rel = sub(y, 0.0);
  sum = add(sum, mult(rel, rel));

  rel = sub(z, -10.0);
  sum = add(sum, mult(rel, rel));
  return mult(sum, 0.1);
}
HNum3 fIdentity(in HNum3 x, in HNum3 y, in HNum3 z,
                   float globalTime) {
  return constH3(0.0);
}


//********************************************************************
vec3 createRay(vec2 px, mat4 PInv, mat4 VInv) {

  // convert pixel to NDS
  // [0,1] -> [-1,1]
  vec2 pxNDS = px * 2. - 1.;

  // choose an arbitrary point in the viewing volume
  // z = -1 equals a point on the near plane, i.e. the screen
  vec3 pointNDS = vec3(pxNDS, -1.);

  // as this is in homogenous space, add the last homogenous coordinate
  vec4 pointNDSH = vec4(pointNDS, 1.0);
  // transform by inverse projection to get the point in view space
  vec4 dirEye = PInv * pointNDSH;

  // since the camera is at the origin in view space by definition,
  // the current point is already the correct direction (dir(0,P) = P - 0 = P
  // as a direction, an infinite point, the homogenous component becomes 0
  // the scaling done by the w-division is not of interest, as the direction
  // in xyz will stay the same and we can just normalize it later
  dirEye.w = 0.;

  // compute world ray direction by multiplying the inverse view matrix
  vec3 dirWorld = (VInv * dirEye).xyz;

  // now normalize direction
  return normalize(dirWorld);
}

// matrix operations
mat4 translate(vec3 t) {
  return mat4(vec4(1., 0., 0., 0.), vec4(0., 1., 0., 0.), vec4(0., 0., 1., 0.),
              vec4(t, 1.));
}
mat4 translateInv(vec3 t) { return translate(-t); }

mat4 scale(vec3 s) {
  return mat4(vec4(s.x, 0., 0., 0.), vec4(0., s.y, 0., 0.),
              vec4(0., 0., s.z, 0.), vec4(0., 0., 0., 1.));
}
mat4 scaleInv(vec3 s) { return scale(1. / s); }

mat4 rightToLeft() {
  // 1 0 0  0
  // 0 1 0  0
  // 0 0 -1 0
  // 0 0 0  1
  return scale(vec3(1., 1., -1.));
}

mat4 rightToLeftInv() {
  // same matrix
  return rightToLeft();
}

mat4 ortho(float l, float r, float b, float t, float n, float f) {

  // translation and scale
  return scale(vec3(2. / (r - l), 2. / (t - b), 2. / (f - n))) *
         translate(vec3(-(l + r) / 2., -(t + b) / 2., -(f + n) / 2.));
}

mat4 orthoInv(float l, float r, float b, float t, float n, float f) {
  return translateInv(vec3(-(l + r) / 2., -(t + b) / 2., -(f + n) / 2.)) *
         scaleInv(vec3(2. / (r - l), 2. / (t - b), 2. / (f - n)));
}

mat4 projection(float n, float f) {
  // n 0 0 0	0
  // 0 n 0 0	0
  // 0 0 n+f	-fn
  // 0 0 1	0
  return mat4(vec4(n, 0., 0., 0.), vec4(0., n, 0., 0.), vec4(0., 0., n + f, 1.),
              vec4(0., 0., -f * n, 0.));
}

mat4 projectionInv(float n, float f) {
  // 1/n 	0 	0 		0
  // 0 	1/n	0 		0
  // 0	0	0 		1
  // 0	0	-1/fn	(f+n)/fn

  return mat4(vec4(1. / n, 0., 0., 0.), vec4(0., 1. / n, 0., 0.),
              vec4(0., 0., 0., -1. / (f * n)),
              vec4(0., 0., 1., (f + n) / (f * n)));
}

mat4 perspective(float fov, float aspect, float n, float f) {
  float l = tan(fov / 2.) * n;
  float b = l / aspect;

  return ortho(-l, l, -b, b, n, f) * projection(n, f) * rightToLeft();
}

mat4 perspectiveInv(float fov, float aspect, float n, float f) {
  float l = tan(fov / 2.) * n;
  float b = l / aspect;

  return rightToLeftInv() * projectionInv(n, f) * orthoInv(-l, l, -b, b, n, f);
}

mat4 lookAt(vec3 eye, vec3 center, vec3 up) {

  vec3 z = normalize(eye - center);
  vec3 x = normalize(cross(up, z));
  vec3 y = cross(z, x);

  mat4 v = mat4(vec4(x.x, y.x, z.x, 0.), vec4(x.y, y.y, z.y, 0.),
                vec4(x.z, y.z, z.z, 0.), vec4(0., 0., 0., 1.));

  return v * translate(-eye);
}

mat4 lookAtInv(vec3 eye, vec3 center, vec3 up) {
  vec3 z = normalize(eye - center);
  vec3 x = normalize(cross(up, z));
  vec3 y = cross(z, x);

  return translateInv(-eye) *
         mat4(vec4(x, 0.), vec4(y, 0.), vec4(z, 0.), vec4(0., 0., 0., 1.));
}

//********************************************************************
//
// INSTRUCTIONS:
//
// This shader implements basic distancefield raytracing of a regular grid on
// a 3D manifold given as the graph of a function f(x,y,z).
//
// Click and drag the mouse to look around. WASD for movement. You are moving
// through distorted space, so you might not move in a straight line, but
// instead move along the shortest paths in curved space.
// Press "R" to reset the view.
//
// WARNING: The space-bending and movement might be a bit nauseating. I tried to
// not make stuff move too fast, but wanted to give a short warning anyways.

// DISPLAY OPTIONS:
// Below you can find a few defines to control the look of the grid, background
// and general tracing
//
// FUNCTIONAL OPTIONS
// The actual function is defined in the common shader part under "Metric
// functions".  A few functions are predefined, but you can implement your own
// in "HessNum3 fCustom(in HessNum3 x, in HessNum3 y, in HessNum3 z, float
// globalTime)". Just comment in the one you want to see or comment all of them
// out to use the custom function. Additional information can be found there as
// well.
//
// Note: As calculating the metric from the function requires 1st and 2nd order
// derivatives and it would be cumbersome to calculate and implement them for
// each function, I use automatic differentiation, which is why you need to use
// the special types and operations defined in the automatic differentiation
// header below to implement your own f.
//
//
//********************************************************************

// Maximum number of trace steps
// Higher step number will let you look farther, especially with higher
// deformations
#define MAX_STEPS 100
// The maximum length traversed per step. Lower values produce higher quality
// paths through the curved space
#define STEP_LENGTH 0.2

// Grid parameters
#define GRID_SIZE 1.0
#define GRID_THICKNESS 0.0001

// Linear fog parameters
#define FOG_MAX_DIST 50.0
#define FOG_MIN_DIST 10.0
#define FOG_COLOR vec3(0.1)

// Surface epsilon
#define EPS 0.005

//********************************************************************
//
// Distance helpers and tracing
//
//********************************************************************
#define INF 3.402823466e+38

float sdf(in vec3 p);

vec3 g_eye;

// Tracing function
bool trace(inout vec3 p, in vec3 dir, out float t, out int steps) {

  steps = 0;
  t = 0.0;
  
  for (int n = 0; n < MAX_STEPS && t < FOG_MAX_DIST; n++) {

    float d = sdf(p);

    if (d < EPS) {
    
     if(length(p - g_eye) > 10.0 * GRID_SIZE) return false;
    
      return true;
    }


    // make step length dependent on distance field to not miss geometry, but
    // also use EPS as a min step
    float s = max(EPS, min(d, STEP_LENGTH));

    // "acceleration" of a particle traveling along the geodesic
    vec3 vp = compute_update(p, dir, iTime);

    // euler integration
    dir += vp * s;
    // position euler step
    p += dir * s;

    // maybe not fully correct, but try to preserve arc length by normalizing
    // dir
    dir = normalize(dir);

    // length traveled, used for fog
    t += s;

    steps++;
  }

  return false;
}

float sdfBox2d(vec2 p, vec2 b);
float sdfGrid(vec3 p, vec3 cellSize, vec3 gridHalfLength);
vec3 normal(in vec3 p);

//********************************************************************
//
// Image formation
//
//********************************************************************

vec3 shade(in vec3 P, in vec3 N, in vec3 color, in vec3 LPos);

float simpleColor(float coord, float gridSize, float repetition);

// vec4 loadValue(in ivec2 re) { return texelFetch(iChannel0, re, 0); }


void mainImage(out vec4 fragColor, in vec2 fragCoord) {

  vec2 uv = fragCoord.xy / iResolution.xy;

  float aspect = iResolution.x / iResolution.y;

  D_thetaPhi = vec4(90.0 * deg2Rad - clamp(2.0 * PI * iTime, 0., 2.0*PI) / 20., 0.0, 0.0, 0.0);

  vec3 rayDir;
  vec3 p;

  // {
    vec2 thetaPhi = D_thetaPhi.xy; // loadValue(D_thetaPhi).xy;
    float theta = thetaPhi.x;
    float phi = thetaPhi.y;

    vec3 eye = D_pos.xyz; // loadValue(D_pos).xyz;
    // vec3 eye = vec3(0.0);
    vec3 center =
        eye + vec3(sin(theta) * sin(phi), cos(theta), sin(theta) * cos(phi));
        
    g_eye = eye;
    // inverse projection and view matrices
    mat4 PInv = perspectiveInv(radians(90.), aspect, 0.1, 100.);
    mat4 VInv = lookAtInv(eye, center, vec3(0., 1., 0.));

    rayDir = createRay(uv, PInv, VInv);
    p = eye;
  // }

  float t = 0.0;
  int steps = 100;

  vec3 col = vec3(0.0);
  
  if (trace(p, rayDir, t, steps)) {
  
    vec3 N = normal(p);

    col =
        vec3(simpleColor(p.x, GRID_SIZE, 3.0), simpleColor(p.y, GRID_SIZE, 7.0),
             simpleColor(p.z, GRID_SIZE, 11.0));

    vec3 LPos = vec3(1.0, 2.0, 0.0);
    col = shade(p, N, col, LPos);

  } else {
    t = INF;
  }

  // simple linear fog
  float fogFactor = (FOG_MAX_DIST - t) / (FOG_MAX_DIST - FOG_MIN_DIST);
  fogFactor = clamp(fogFactor, 0.0, 1.0);
  col = mix(FOG_COLOR, col, fogFactor);

  fragColor = vec4(col, t == INF ? 0.0 : 1.0);
}

//***************

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

float sdfGrid(vec3 p, vec3 cellSize, vec3 gridHalfLength) {
  vec3 pg = floor(p / cellSize) * cellSize;

  vec3 center = pg + cellSize * 0.5;

  float d = INF;
  
  d = min(d, sdfBox2d(p.xz - center.xz, gridHalfLength.xz));
  d = min(d, sdfBox2d(p.xy - center.xy, gridHalfLength.xy));
  d = min(d, sdfBox2d(p.yz - center.yz, gridHalfLength.yz));
  
  return d;
}

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

float sdf(in vec3 p) {
  float d = sdfGrid(p, vec3(GRID_SIZE), vec3(GRID_THICKNESS));
  
  // d = max(d, pSphere(p - g_center, 5.0 * GRID_SIZE));
  
  return d;
}

vec3 normal(in vec3 p) {
  return normalize(
      vec3(sdf(p + vec3(EPS, 0., 0.)) - sdf(p - vec3(EPS, 0., 0.)),
           sdf(p + vec3(0., EPS, 0.)) - sdf(p - vec3(0., EPS, 0.)),
           sdf(p + vec3(0., 0., EPS)) - sdf(p - vec3(0., 0., EPS))));
}

//********************************************************************
//
// Shading and main
//
//********************************************************************
vec3 shade(in vec3 P, in vec3 N, in vec3 color, in vec3 LPos) {
  vec3 L = normalize(LPos - P);
  // simple two sided diffuse shading
  return vec3(0.1, 0.8, 0.1) * clamp(abs(dot(L, N)), 0.0, 1.0);
}

float simpleColor(float coord, float gridSize, float repetition) {

  float cr = (gridSize * 0.5 - coord) / gridSize;

  float v = pow(sin(PI * cr / 2.0), 2.0);

  return v;
}