6precision mediump float;    
uniform vec2 u_resolution;
uniform float u_time;
uniform sampler2D texture_sampler;
uniform sampler2D texts;

struct Ray {
    vec3 rd;
    vec3 dir;    
};

mat2 scale(vec2 scale){
    return mat2(1. / scale.x, 0.0, 0.0, 1./scale.y);
}

////////////////////////////////////////////////////////////////
//
//                           HG_SDF
//
//     GLSL LIBRARY FOR BUILDING SIGNED DISTANCE BOUNDS
//
//     version 2021-07-28
//
//     Check https://mercury.sexy/hg_sdf for updates
//     and usage examples. Send feedback to spheretracing@mercury.sexy.
//
//     Brought to you by MERCURY https://mercury.sexy/
//
//
//
// Released dual-licensed under
//   Creative Commons Attribution-NonCommercial (CC BY-NC)
// or
//   MIT License
// at your choice.
//
// SPDX-License-Identifier: MIT OR CC-BY-NC-4.0
//
// /////


////////////////////////////////////////////////////////////////
//
//             HELPER FUNCTIONS/MACROS
//
////////////////////////////////////////////////////////////////

const float PI = 3.14159265;
const float  TAU = (2.*PI);
const float PHI = sqrt(5.)*0.5 + 0.5;

// Sign function that doesn't return 0
float sgn(float x) {
	return (x < 0. )? -1. : 1.;
}

vec2 sgn(vec2 v) {
	return vec2((v.x<0.)?-1.:1., (v.y<0.)?-1.:1.);
}

float square (float x) {
	return x*x;
}

vec2 square (vec2 x) {
	return x*x;
}

vec3 square (vec3 x) {
	return x*x;
}

float lengthSqr(vec3 x) {
	return dot(x, x);
}


// Maximum/minumum elements of a vector
float vmax(vec2 v) {
	return max(v.x, v.y);
}

float vmax(vec3 v) {
	return max(max(v.x, v.y), v.z);
}

float vmax(vec4 v) {
	return max(max(v.x, v.y), max(v.z, v.w));
}

float vmin(vec2 v) {
	return min(v.x, v.y);
}

float vmin(vec3 v) {
	return min(min(v.x, v.y), v.z);
}

float vmin(vec4 v) {
	return min(min(v.x, v.y), min(v.z, v.w));
}

////////////////////////////////////////////////////////////////
//
//             PRIMITIVE DISTANCE FUNCTIONS
//
////////////////////////////////////////////////////////////////
//
// Conventions:
//
// Everything that is a distance function is called fSomething.
// The first argument is always a point in 2 or 3-space called <p>.
// Unless otherwise noted, (if the object has an intrinsic "up"
// side or direction) the y axis is "up" and the object is
// centered at the origin.
//
////////////////////////////////////////////////////////////////

float fSphere(vec3 p, float r) {
	return length(p) - r;
}

// Plane with normal n (n is normalized) at some distance from the origin
float fPlane(vec3 p, vec3 n, float distanceFromOrigin) {
	return dot(p, n) + distanceFromOrigin;
}

// Cheap Box: distance to corners is overestimated
float fBoxCheap(vec3 p, vec3 b) { //cheap box
	return vmax(abs(p) - b);
}

// Box: correct distance to corners
float fBox(vec3 p, vec3 b) {
	vec3 d = abs(p) - b;
	return length(max(d, vec3(0))) + vmax(min(d, vec3(0)));
}

// Same as above, but in two dimensions (an endless box)
float fBox2Cheap(vec2 p, vec2 b) {
	return vmax(abs(p)-b);
}

float fBox2(vec2 p, vec2 b) {
	vec2 d = abs(p) - b;
	return length(max(d, vec2(0.))) + vmax(min(d, vec2(0.)));
}


// Endless "corner"
float fCorner (vec2 p) {
	return length(max(p, vec2(0.))) + vmax(min(p, vec2(0.)));
}

// Cylinder standing upright on the xz plane
float fCylinder(vec3 p, float r, float height) {
	float d = length(p.xz) - r;
	d = max(d, abs(p.y) - height);
	return d;
}

// Capsule: A Cylinder with round caps on both sides
float fCapsule(vec3 p, float r, float c) {
	return mix(length(p.xz) - r, length(vec3(p.x, abs(p.y) - c, p.z)) - r, step(c, abs(p.y)));
}

// Distance to line segment between <a> and <b>, used for fCapsule() version 2below
float fLineSegment(vec3 p, vec3 a, vec3 b) {
	vec3 ab = b - a;
	float t = clamp( dot(p - a, ab) / dot(ab, ab), 0., 1. );
	return length((ab*t + a) - p);
}

// Capsule version 2: between two end points <a> and <b> with radius r 
float fCapsule(vec3 p, vec3 a, vec3 b, float r) {
	return fLineSegment(p, a, b) - r;
}

// Torus in the XZ-plane
float fTorus(vec3 p, float smallRadius, float largeRadius) {
	return length(vec2(length(p.xz) - largeRadius, p.y)) - smallRadius;
}

// A circle line. Can also be used to make a torus by subtracting the smaller radius of the torus.
float fCircle(vec3 p, float r) {
	float l = length(p.xz) - r;
	return length(vec2(p.y, l));
}

// A circular disc with no thickness (i.e. a cylinder with no height).
// Subtract some value to make a flat disc with rounded edge.
float fDisc(vec3 p, float r) {
	float l = length(p.xz) - r;
	return l < 0. ? abs(p.y) : length(vec2(p.y, l));
}

// Hexagonal prism, circumcircle variant
float fHexagonCircumcircle(vec3 p, vec2 h) {
	vec3 q = abs(p);
	return max(q.y - h.y, max(q.x*sqrt(3.)*0.5 + q.z*0.5, q.z) - h.x);
	//this is mathematically equivalent to this line, but less efficient:
	//return max(q.y - h.y, max(dot(vec2(cos(PI/3), sin(PI/3)), q.zx), q.z) - h.x);
}

// Hexagonal prism, incircle variant
float fHexagonIncircle(vec3 p, vec2 h) {
	return fHexagonCircumcircle(p, vec2(h.x*sqrt(3.)*0.5, h.y));
}

// Cone with correct distances to tip and base circle. Y is up, 0 is in the middle of the base.
float fCone(vec3 p, float radius, float height) {
	vec2 q = vec2(length(p.xz), p.y);
	vec2 tip = q - vec2(0, height);
	vec2 mantleDir = normalize(vec2(height, radius));
	float mantle = dot(tip, mantleDir);
	float d = max(mantle, -q.y);
	float projected = dot(tip, vec2(mantleDir.y, -mantleDir.x));
	
	// distance to tip
	if ((q.y > height) && (projected < 0.)) {
		d = max(d, length(tip));
	}
	
	// distance to base ring
	if ((q.x > radius) && (projected > length(vec2(height, radius)))) {
		d = max(d, length(q - vec2(radius, 0)));
	}
	return d;
}

////////////////////////////////////////////////////////////////
//
//                DOMAIN MANIPULATION OPERATORS
//
////////////////////////////////////////////////////////////////
//
// Conventions:
//
// Everything that modifies the domain is named pSomething.
//
// Many operate only on a subset of the three dimensions. For those,
// you must choose the dimensions that you want manipulated
// by supplying e.g. <p.x> or <p.zx>
//
// <inout p> is always the first argument and modified in place.
//
// Many of the operators partition space into cells. An identifier
// or cell index is returned, if possible. This return value is
// intended to be optionally used e.g. as a random seed to change
// parameters of the distance functions inside the cells.
//
// Unless stated otherwise, for cell index 0, <p> is unchanged and cells
// are centered on the origin so objects don't have to be moved to fit.
//
//
////////////////////////////////////////////////////////////////

// Rotate around a coordinate axis (i.e. in a plane perpendicular to that axis) by angle <a>.
// Read like this: R(p.xz, a) rotates "x towards z".
// This is fast if <a> is a compile-time constant and slower (but still practical) if not.
void pR(inout vec2 p, float a) {
	p = cos(a)*p + sin(a)*vec2(p.y, -p.x);
}

// Shortcut for 45-degrees rotation
void pR45(inout vec2 p) {
	p = (p + vec2(p.y, -p.x))*sqrt(0.5);
}

// Repeat space along one axis. Use like this to repeat along the x axis:
// <float cell = pMod1(p.x,5);> - using the return value is optional.
float pMod1(inout float p, float size) {
	float halfsize = size*0.5;
	float c = floor((p + halfsize)/size);
	p = mod(p + halfsize, size) - halfsize;
	return c;
}

// Same, but mirror every second cell so they match at the boundaries
float pModMirror1(inout float p, float size) {
	float halfsize = size*0.5;
	float c = floor((p + halfsize)/size);
	p = mod(p + halfsize,size) - halfsize;
	p *= mod(c, 2.0)*2. - 1.;
	return c;
}

// Repeat the domain only in positive direction. Everything in the negative half-space is unchanged.
float pModSingle1(inout float p, float size) {
	float halfsize = size*0.5;
	float c = floor((p + halfsize)/size);
	if (p >= 0.)
		p = mod(p + halfsize, size) - halfsize;
	return c;
}

// Repeat only a few times: from indices <start> to <stop> (similar to above, but more flexible)
float pModInterval1(inout float p, float size, float start, float stop) {
	float halfsize = size*0.5;
	float c = floor((p + halfsize)/size);
	p = mod(p+halfsize, size) - halfsize;
	if (c > stop) { //yes, this might not be the best thing numerically.
		p += size*(c - stop);
		c = stop;
	}
	if (c <start) {
		p += size*(c - start);
		c = start;
	}
	return c;
}


// Repeat around the origin by a fixed angle.
// For easier use, num of repetitions is use to specify the angle.
float pModPolar(inout vec2 p, float repetitions) {
	float angle = 2.*PI/repetitions;
	float a = atan(p.y, p.x) + angle/2.;
	float r = length(p);
	float c = floor(a/angle);
	a = mod(a,angle) - angle/2.;
	p = vec2(cos(a), sin(a))*r;
	// For an odd number of repetitions, fix cell index of the cell in -x direction
	// (cell index would be e.g. -5 and 5 in the two halves of the cell):
	if (abs(c) >= (repetitions/2.)) c = abs(c);
	return c;
}

// Repeat in two dimensions
vec2 pMod2(inout vec2 p, vec2 size) {
	vec2 c = floor((p + size*0.5)/size);
	p = mod(p + size*0.5,size) - size*0.5;
	return c;
}

// Same, but mirror every second cell so all boundaries match
vec2 pModMirror2(inout vec2 p, vec2 size) {
	vec2 halfsize = size*0.5;
	vec2 c = floor((p + halfsize)/size);
	p = mod(p + halfsize, size) - halfsize;
	p *= mod(c,vec2(2.))*2. - vec2(1);
	return c;
}

// Same, but mirror every second cell at the diagonal as well
vec2 pModGrid2(inout vec2 p, vec2 size) {
	vec2 c = floor((p + size*0.5)/size);
	p = mod(p + size*0.5, size) - size*0.5;
	p *= mod(c,vec2(2.))*2. - vec2(1.);
	p -= size/2.;
	if (p.x > p.y) p.xy = p.yx;
	return floor(c/2.);
}

// Repeat in three dimensions
vec3 pMod3(inout vec3 p, vec3 size) {
	vec3 c = floor((p + size*0.5)/size);
	p = mod(p + size*0.5, size) - size*0.5;
	return c;
}

// Mirror at an axis-aligned plane which is at a specified distance <dist> from the origin.
float pMirror (inout float p, float dist) {
	float s = sgn(p);
	p = abs(p)-dist;
	return s;
}

// Mirror in both dimensions and at the diagonal, yielding one eighth of the space.
// translate by dist before mirroring.
vec2 pMirrorOctant (inout vec2 p, vec2 dist) {
	vec2 s = sgn(p);
	pMirror(p.x, dist.x);
	pMirror(p.y, dist.y);
	if (p.y > p.x)
		p.xy = p.yx;
	return s;
}

// Reflect space at a plane
float pReflect(inout vec3 p, vec3 planeNormal, float offset) {
	float t = dot(p, planeNormal)+offset;
	if (t < 0.) {
		p = p - (2.*t)*planeNormal;
	}
	return sgn(t);
}


////////////////////////////////////////////////////////////////
//
//             OBJECT COMBINATION OPERATORS
//
////////////////////////////////////////////////////////////////
//
// We usually need the following boolean operators to combine two objects:
// Union: OR(a,b)
// Intersection: AND(a,b)
// Difference: AND(a,!b)
// (a and b being the distances to the objects).
//
// The trivial implementations are min(a,b) for union, max(a,b) for intersection
// and max(a,-b) for difference. To combine objects in more interesting ways to
// produce rounded edges, chamfers, stairs, etc. instead of plain sharp edges we
// can use combination operators. It is common to use some kind of "smooth minimum"
// instead of min(), but we don't like that because it does not preserve Lipschitz
// continuity in many cases.
//
// Naming convention: since they return a distance, they are called fOpSomething.
// The different flavours usually implement all the boolean operators above
// and are called fOpUnionRound, fOpIntersectionRound, etc.
//
// The basic idea: Assume the object surfaces intersect at a right angle. The two
// distances <a> and <b> constitute a new local two-dimensional coordinate system
// with the actual intersection as the origin. In this coordinate system, we can
// evaluate any 2D distance function we want in order to shape the edge.
//
// The operators below are just those that we found useful or interesting and should
// be seen as examples. There are infinitely more possible operators.
//
// They are designed to actually produce correct distances or distance bounds, unlike
// popular "smooth minimum" operators, on the condition that the gradients of the two
// SDFs are at right angles. When they are off by more than 30 degrees or so, the
// Lipschitz condition will no longer hold (i.e. you might get artifacts). The worst
// case is parallel surfaces that are close to each other.
//
// Most have a float argument <r> to specify the radius of the feature they represent.
// This should be much smaller than the object size.
//
// Some of them have checks like "if ((-a < r) && (-b < r))" that restrict
// their influence (and computation cost) to a certain area. You might
// want to lift that restriction or enforce it. We have left it as comments
// in some cases.
//
// usage example:
//
// float fTwoBoxes(vec3 p) {
//   float box0 = fBox(p, vec3(1));
//   float box1 = fBox(p-vec3(1), vec3(1));
//   return fOpUnionChamfer(box0, box1, 0.2);
// }
//
////////////////////////////////////////////////////////////////


// The "Chamfer" flavour makes a 45-degree chamfered edge (the diagonal of a square of size <r>):
float fOpUnionChamfer(float a, float b, float r) {
	return min(min(a, b), (a - r + b)*sqrt(0.5));
}

// Intersection has to deal with what is normally the inside of the resulting object
// when using union, which we normally don't care about too much. Thus, intersection
// implementations sometimes differ from union implementations.
float fOpIntersectionChamfer(float a, float b, float r) {
	return max(max(a, b), (a + r + b)*sqrt(0.5));
}

// Difference can be built from Intersection or Union:
float fOpDifferenceChamfer (float a, float b, float r) {
	return fOpIntersectionChamfer(a, -b, r);
}

// The "Round" variant uses a quarter-circle to join the two objects smoothly:
float fOpUnionRound(float a, float b, float r) {
	vec2 u = max(vec2(r - a,r - b), vec2(0));
	return max(r, min (a, b)) - length(u);
}

float fOpIntersectionRound(float a, float b, float r) {
	vec2 u = max(vec2(r + a,r + b), vec2(0));
	return min(-r, max (a, b)) + length(u);
}

float fOpDifferenceRound (float a, float b, float r) {
	return fOpIntersectionRound(a, -b, r);
}


// The "Columns" flavour makes n-1 circular columns at a 45 degree angle:
float fOpUnionColumns(float a, float b, float r, float n) {
	if ((a < r) && (b < r)) {
		vec2 p = vec2(a, b);
		float columnradius = r*sqrt(2.)/((n-1.)*2.+sqrt(2.));
		pR45(p);
		p.x -= sqrt(2.)/2.*r;
		p.x += columnradius*sqrt(2.);
		if (mod(n,2.) == 1.) {
			p.y += columnradius;
		}
		// At this point, we have turned 45 degrees and moved at a point on the
		// diagonal that we want to place the columns on.
		// Now, repeat the domain along this direction and place a circle.
		pMod1(p.y, columnradius*2.);
		float result = length(p) - columnradius;
		result = min(result, p.x);
		result = min(result, a);
		return min(result, b);
	} else {
		return min(a, b);
	}
}

float fOpDifferenceColumns(float a, float b, float r, float n) {
	a = -a;
	float m = min(a, b);
	//avoid the expensive computation where not needed (produces discontinuity though)
	if ((a < r) && (b < r)) {
		vec2 p = vec2(a, b);
		float columnradius = r*sqrt(2.)/n/2.0;
		columnradius = r*sqrt(2.)/((n-1.)*2.+sqrt(2.));

		pR45(p);
		p.y += columnradius;
		p.x -= sqrt(2.)/2.*r;
		p.x += -columnradius*sqrt(2.)/2.;

		if (mod(n,2.) == 1.) {
			p.y += columnradius;
		}
		pMod1(p.y,columnradius*2.);

		float result = -length(p) + columnradius;
		result = max(result, p.x);
		result = min(result, a);
		return -min(result, b);
	} else {
		return -m;
	}
}

float fOpIntersectionColumns(float a, float b, float r, float n) {
	return fOpDifferenceColumns(a,-b,r, n);
}

// The "Stairs" flavour produces n-1 steps of a staircase:
// much less stupid version by paniq
float fOpUnionStairs(float a, float b, float r, float n) {
	float s = r/n;
	float u = b-r;
	return min(min(a,b), 0.5 * (u + a + abs ((mod (u - a + s, 2. * s)) - s)));
}

// We can just call Union since stairs are symmetric.
float fOpIntersectionStairs(float a, float b, float r, float n) {
	return -fOpUnionStairs(-a, -b, r, n);
}

float fOpDifferenceStairs(float a, float b, float r, float n) {
	return -fOpUnionStairs(-a, b, r, n);
}


// Similar to fOpUnionRound, but more lipschitz-y at acute angles
// (and less so at 90 degrees). Useful when fudging around too much
// by MediaMolecule, from Alex Evans' siggraph slides
float fOpUnionSoft(float a, float b, float r) {
	float e = max(r - abs(a - b), 0.);
	return min(a, b) - e*e*0.25/r;
}


// produces a cylindical pipe that runs along the intersection.
// No objects remain, only the pipe. This is not a boolean operator.
float fOpPipe(float a, float b, float r) {
	return length(vec2(a, b)) - r;
}

// first object gets a v-shaped engraving where it intersect the second
float fOpEngrave(float a, float b, float r) {
	return max(a, (a + r - abs(b))*sqrt(0.5));
}

// first object gets a capenter-style groove cut out
float fOpGroove(float a, float b, float ra, float rb) {
	return max(a, min(a + ra, rb - abs(b)));
}

// first object gets a capenter-style tongue attached
float fOpTongue(float a, float b, float ra, float rb) {
	return min(a, max(a - ra, abs(b) - rb));
}

//#endSection End of library

// https://stackoverflow.com/questions/4200224/random-noise-functions-for-glsl
// golden_noise
float noise(in vec2 xy, in float seed){
       return fract(tan(distance(xy*PHI, xy)*seed)*xy.x);
}

vec3 rnd23(vec2 p)
{
	vec3 p3 = fract(p.xyx * vec3(.1031, .1030, .0973));
    p3 += dot(p3, p3.yxz+33.33);
    return fract((p3.xxy+p3.yzz)*p3.zyx);
}

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

float opExtrusion( in vec3 p, in float sdf, in float h )
{
    vec2 w = vec2( sdf, abs(p.z) - h);
  	return min(max(w.x,w.y),0.0) + length(max(w,0.0));
}

float sdCog2d(vec2 pos) {        
    float r = length(pos)*2.;
    float a = atan(pos.y,pos.x);	    
    float f = 1. - smoothstep(-0.2, .8, sin(a * 12.))*0.14;
    f =  smoothstep(f,f + 2.,r);       
    return f;
}

float sdCog(vec3 pos, float angle) {        
    pos.xy *= Rot(angle);    
    float d1 = opExtrusion(pos, sdCog2d(pos.xy), 0.05);    
    float d2 = fCapsule(pos, vec3(0., 0.0, 0.), vec3(0., 0., 1.), 0.2);    
    return 0.8 * fOpDifferenceRound(d1,d2,0.05)-0.003;}


float sdHex(vec3 pos, float i, float angle) {                   
    vec3 po = pos;
    
    po.xz *= Rot(angle);
    po.yz *= Rot(angle);    
    pR(po.yz, PI/2.);

    float d1 = fHexagonCircumcircle(po, vec2(0.5+i, .1));
    float d2 = fHexagonCircumcircle(po, vec2(0.2+i, .1));    
    return fOpDifferenceRound(d1,d2,0.1);

}

// Scene
vec2 mapScene(in vec3 p) {
    float mat = 0.;
    float d = 1e10;

    //float dGround = p.y + 2.5;
    //d = min(d, dGround);

    vec3 po = p;    
    //po.y += sin(u_time);
    // pMod3(po, vec3(3.));    
    po.xy *= scale(vec2(1.3, 1.3));

    const float num = 6.;
    for (float i = 1.; i <= num; i++) {        
        // pos.z += i*.1;        
        float a = sdHex(po,i*0.35, u_time + abs( 2. + 0.4 * sin(u_time)) * i*3.1415/num);
        d = min(d,a);
        if (d == a) mat = 1. + mod(i,3.);
    }

    /*float c2 = sdCog(p, u_time);
    d = min(d, c2);
    if ( d == c2) mat = 4.;
*/
  //  float c3 = fBox(p+vec3(0.5, .87, 0.), vec3(1., 1.,1.));    
  //  d = min(d, c3);
   
 

    // if ( d == c1) mat = 1.;
 
    //if ( d == c3) mat = 3.;

    return vec2(d, mat); 
}

vec3 castRay(vec3 ro, vec3 rd, inout vec3 pos) {
	float t = 0.0;
    float mat = 0.;
    float hit = 0.;    
    for(int i=0; i < 150; i++) {
    	pos = ro + rd * t;
        vec2 res = mapScene(pos);
        t += res.x;
        mat = res.y;
        if (t > 80.) break;
        if (res.x < abs(0.001*t) ) {
            hit = 1.;
            break;
        }        
    }
    if (t > 80.) t = -1.0;
    return vec3(t, mat, hit);
}


vec3 castReflectedRay(vec3 ro, vec3 rd, vec3 pos) {
	float t = 0.0;
    float mat = 0.;
    float hit = 0.;    
    for(int i=0; i < 50; i++) {
    	pos = ro + rd * t;         
        vec2 res = mapScene(pos);
        t += res.x;
        mat = res.y;                        
        if (t > 40.) break;
        if (res.x < abs(0.001*t) ) {
            hit = 1.;
            break;
        }
        
    }
    if (t > 40.) t = -1.0;
    return vec3(t, mat, hit);
}

float softshadow( in vec3 ro, in vec3 rd, float mint, float maxt, float w )
{
    float res = 1.0;
    float t = mint;
    for( int i=0; i<40; i++ )
    {
        if (t > maxt) break;
        float h = mapScene(ro + t*rd).x;
        res = min( res, h/(w*t) );
        t += clamp(h, 0.005, 0.50);
        if( res < -1.0 || t>maxt ) break;

    }
    res = max(res,-1.0);
    return 0.25*(1.0+res)*(1.0+res)*(2.0-res);    
}

float castShadow(vec3 ro, vec3 rd) {
    float res = 1.0;
    float t = 0.001;
     for(int i = 0; i < 40; i++) {    	         
        float h = mapScene(ro + t* rd).x;
        res = min(res, 10.0*h/t);
        if (abs(h) < (0.001*t) ) break;
        t += h;
        if (t > 20.) break;        
    }
    return clamp(res,0., 1.);
}

vec3 calcNormal(vec3 pos) {
    vec2 e = vec2(.001, 0.);
    vec3 n = vec3(  mapScene(pos+e.xyy).x - mapScene(pos-e.xyy).x,
                    mapScene(pos+e.yxy).x - mapScene(pos-e.yxy).x,
                    mapScene(pos+e.yyx).x - mapScene(pos-e.yyx).x
                ); 
    return normalize(n);
}

vec3 fresnel( vec3 F0, vec3 h, vec3 l ) {
	return F0 + ( 1.0 - F0 ) * pow( clamp( 1.0 - dot( h, l ), 0.0, 1.0 ), 5.0 );
}

// https://suricrasia.online/blog/shader-functions/
vec3 eRot(vec3 p, vec3 ax, float ro) {
    return mix(dot(ax,p)*ax, p, cos(ro)) + sin(ro)*cross(ax,p);
}



vec3 addPointLight(vec3 light_pos, vec3 light_color, float shininess, vec3 v, vec3 dir, vec3 n, float occ) {
    vec3 Ks = vec3( .5454 );
    vec3 Kd = vec3( 1. );
    vec3 ref = reflect( dir, n  );       
	vec3 vl = normalize( v );	
	vec3 diffuse  = Kd * vec3( max( 0.0, dot( vl, n ) ) );
	vec3 specular = vec3( max( 0.0, dot( vl, ref ) ) );		
    vec3 F = fresnel( Ks, normalize( vl - dir ), vl )*occ;
	float shadow = softshadow(v+n*0.01,light_pos, .01, 30., 18.);  
    //float shadow = castShadow(v + n*0.02, light_pos);
    specular = pow( specular, vec3( shininess ) )*occ;    
	return light_color * mix( diffuse, specular, F ) * shadow;

}

float getAmbientOcc(vec3 p, vec3 n) {
    float occ = 0.;
    float weight = 1.;
    for (int i = 0; i < 8; i++) {
        float len = 0.01 + 0.02 * float(i*i);
        float dist = mapScene(p+n*len).x;
        occ += (len - dist) * weight;
        weight *=0.85;
    }
    return 1.0 - clamp(0.6 * occ, 0., 1.);
}

vec3 shading(vec3 v, vec3 n, vec3 dir, float material) {
	float shininess = 1.;
	float occ = getAmbientOcc(v,n);
    vec3 outMaterial = vec3(0.1529, 0.1529, 0.1529);

    if (material == 0.) {
        outMaterial = vec3(0.2863, 0.1059, 0.2431);
        shininess = 1.5;
    } else if (material == 1.) {
        outMaterial = vec3(0.3294, 0.0941, 0.6);
        shininess = 0.6;
    } else if (material == 2.) {
        outMaterial = vec3(0.5804, 0.9647, 1.0);
        shininess = 1.;
    } else if (material == 3.) {
        outMaterial = vec3(0.0, 0.0, 0.0);
        shininess = 100.;
    } else if (material == 4.) {
        outMaterial = vec3(0.9961, 1.0, 0.9922);
        shininess = .3;
    } 
    
    vec3 lights = vec3(0.);     
    lights += addPointLight(vec3( 0., -40., -1. ),vec3(0.56, 0.44, 0.18)*2., shininess, v, dir, n, occ);
    lights += addPointLight(vec3( -20.,4., 10. ),vec3(0.04, 0.2, 0.71)*2., shininess, v, dir, n, occ);
   //  lights += addPointLight(vec3( 2., -1.0, -1.0 ),vec3(0.12, 0.51, 0.63)*2., shininess, v,dir,n,occ );
    
    
    vec3 lightDir = vec3(0. , 4., 2.);
    float sun_dif = clamp(dot(n, lightDir), 0., 1.);
    float shadow = softshadow(v+n*0.01,lightDir, .01, 30., 18.);
    lights += vec3(0.6627, 0.7098, 0.8863) * sun_dif*shadow*occ; //* shadow; //* mix( vec3(sun_dif), specular, F )
        
    
     // final += texture( iChannel0, ref ).rgb * fresnel( Ks, n, -dir );
     // vec3 col = vec3(0.4)* ref.x;

    float ind = clamp( dot( n, normalize(lightDir*vec3(-1.0,.0,-1.0)) ), 0.0, 1.0 );
    lights += vec3(0.1333, 0.1333, 0.1216) * ind *occ; 

    return outMaterial * max(vec3(0.), lights);
}

vec3 postProcess(vec3 col) {             
    // float random = noise(gl_FragCoord.xy, 0.01+u_time);
    // float random2 = noise(gl_FragCoord.xy, .2+u_time);    
    //col += 0.075*clamp(vec3(0.5*random, 0.5*random2, 0.5*random), 0.02, 1.); // dither

    // Normalized pixel coordinates (from 0 to 1)
    vec2 screenCoord = gl_FragCoord.xy/u_resolution.xy;

    // Vignette
    float radius = 0.8;
    float d = smoothstep(radius, radius-0.4, length(screenCoord-vec2(0.5)));
    col = mix(col, col * d, .9);
    
    // Contrast
    float constrast = .5;
    col = mix(col, smoothstep(0.0, 1.0, col), constrast);
    

    // Colour mapping
    col *= vec3(1.0, 1.0, 1.0);

    col = pow( col, vec3(1.0/2.2) ); // gamma
    
    // fade in at the beginning
    //col*=vec3(clamp((u_time-1.8)*0.5,0., 1.));
    
    // fade out at the end
    // col*=vec3(clamp((120.-u_time)*.35, 0., 1.));
    
    return col;
}

vec3 getCameraRayDir(vec2 uv, vec3 camPos, vec3 camTarget)
{
    // Calculate camera's "orthonormal basis", i.e. its transform matrix components
    vec3 camForward = normalize(camTarget-camPos );
    vec3 camRight = normalize(cross(vec3(.0, 1.0, 0.0), camForward));
    vec3 camUp = normalize(cross(camForward, camRight));     
    float fov = 0.7;
    vec3 vDir = normalize(uv.x * camRight + uv.y * camUp + camForward * fov); 
    return vDir;
}


vec3 getCameraRayDir2(vec2 uv, vec3 camPos, vec3 lookAt, float zoom){
    vec3 f = normalize(lookAt - camPos);
    vec3 r = cross(vec3(0.0,1.0,0.0),f);
    vec3 u = cross(f,r);    
    vec3 c=camPos+f*zoom;
    vec3 i=c+uv.x*r+uv.y*u;
    return normalize(i-camPos);    
}

vec3 getCameraFov(vec2 uv, vec3 camPos, vec3 camTarget) {
    vec3 camForward = normalize(camTarget-camPos);
    vec3 camRight =  normalize(cross(vec3(0.0, 1.0, 0.0), camForward));
    vec3 camUp = normalize(cross(camForward,camRight));
	float fov = 1.7;
	// Depth of field 
    float dof = .25;
	vec2 h = vec2( noise(gl_FragCoord.xy, .13), noise(gl_FragCoord.xy, .4));
    //vec3 h= rnd23(gl_FragCoord.xy);
	vec3 voff = sqrt(h.x)*(camRight*sin(h.y*6.283)+camUp*cos(h.y*6.283))*dof;
	// camTarget -=voff;
    float focusdistance = 150.2;
    return normalize(uv.x * camRight + uv.y * camUp + fov * camForward + voff * fov/focusdistance);
}

vec3 applyFog(vec3 col, float t, vec3 rd, vec3  lightDir, float b ) {
    float fogAmount = 1.0 - exp(-t*b);
    float sunAmount = max( dot(rd, lightDir), 0.0 );
    vec3  fogColor  = mix( vec3(0.3686, 0.2431, 0.4392), // blue
                           vec3(0.4, 0.7294, 0.9216), // yellow
                           pow(sunAmount,8.0) );
    return mix( col, fogColor, fogAmount );
}

vec3 render(vec2 uv) {

    bool useDof = !true;
    //vec2 uv = (2.0 * gl_FragCoord.xy - u_resolution.xy) / u_resolution.y;        

    float angle = -2.4  +u_time*0.4;
    //angle = 0.;

    // camera
    vec3 camPos = vec3(0., 5., -3.);    
    vec3 camTarget = vec3(0., 0., 0.);     
    vec3 rayDir;
    
    if (useDof) {
        rayDir = getCameraFov(uv, camPos, camTarget);
    } else {    
        rayDir = getCameraRayDir2(uv, camPos, camTarget, 1.0);
    }    
    vec3 col =  vec3(0.051, 0.0667, 0.1529);
    //vec3  col = vec3(0.0314, 0.0118, 0.1255) + rayDir.y * 0.4;
    vec3 hitPos = vec3(0.);

    vec3 t = castRay(camPos, rayDir, hitPos);        
    vec3 rd = rayDir;

    if (t.z > 0.) {                
        vec3 nor = calcNormal(hitPos);
        col = shading(hitPos, nor, rayDir , t.y);      
        float fogAmount = 0.04;
        col = col*exp(-t.x*fogAmount) + applyFog(col, t.x, rd, vec3(0., .3, -1.), fogAmount) * (1.0-exp(-t.x*fogAmount));

        rayDir = normalize(reflect(rayDir, nor));
        vec3 rayOrigin = hitPos + (rayDir * 0.01);        
        vec3 t2 = castReflectedRay(rayOrigin, rayDir, hitPos);
    
        if (t2.z > 0.) {        
            hitPos =  rayOrigin + rayDir * t2.x;
            nor = calcNormal(hitPos);
            col += 0.1 * shading(hitPos, nor, rayDir , t2.y);
            
           /* rayDir = normalize(reflect(rayDir, nor));            
            rayOrigin = hitPos + (rayDir * 0.01);
            vec3 t3 = castReflectedRay(rayOrigin, rayDir, hitPos);
        
            if (t3.z > 0.) {        
                hitPos =  rayOrigin + rayDir * t3.x;
                nor = calcNormal(hitPos);
                col += 0.025 * shading(hitPos, nor, rayDir , t3.y);
            } */ 
        }        
    }
    // pixelColor*exp(-distance*b) + fogColor*(1.0-exp(-distance*b));

    return col;
}
vec2 getUV(vec2 offset) {
    vec2 uv =  2.0 *((gl_FragCoord.xy + offset*0.5)/u_resolution.xy - 0.5);
    uv.x *= u_resolution.x/u_resolution.y; // Correct for aspect ratio
    return uv;   
}

void main()
{
   
    vec3 finalColor = vec3(0.);
    const float AA_SIZE = 1.;
    float count = 0.0;
/*
    for (float aaY = 0.0; aaY < AA_SIZE; aaY++) {
        for (float aaX = 0.0; aaX < AA_SIZE; aaX++) {            
            finalColor += render(getUV(vec2(  aaX, aaY)));            
            count += 1.0;
        }
    }   
    finalColor /= count;        */
    finalColor += render(getUV(vec2(  0.,0.)));            
    
    finalColor = postProcess(finalColor);

    gl_FragColor = vec4(finalColor, 1.);

}