///////////////////////////////////////////////////////////////////////////////////
// NOTE: The music should start with the shader, if it does not click play then  //
// stop on the sound in channel2, then click reset time button "|<" then play |> //
// you should be able to get the shader and music to start together              //
///////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

#define TOO_FAR 100000000.0
#define EPSILON 0.0001
#define PI 3.14159

vec3 XAXIS = vec3(1.0, 0.0, 0.0);
vec3 YAXIS = vec3(0.0, 1.0, 0.0);
vec3 ZAXIS = vec3(0.0, 0.0, 1.0);

///////////////////////////////////////////////////////////////////////////////////

struct Material {
    vec3 colour;
    float diffuse;
    float specular;
};
    
struct Ray {
    vec3 pos;
    vec3 dir;
};

struct Light {
    vec3 dir;
    vec3 colour;
};
            
struct Result {
    vec3 pos;
    vec3 normal;
    float dist;
    float mindist;
    float travelled;
    Material mat;
};
    
struct SDFResult
{
	float dist;
    int matindex;
};

///////////////////////////////////////////////////////////////////////////////////
    
mat3 rotationmatrix(vec3 a)
{
    float cp=cos(a.x);
    float sp=sin(a.x);
    float cy=cos(a.y);
    float sy=sin(a.y);
    float cr=cos(a.z);
    float sr=sin(a.z);
    mat3 pitch = mat3(1, 0, 0, 0, cp, sp, 0, -sp, cp);
    mat3 yaw = mat3(cy, 0, -sy, 0, 1, 0, sy, 0, cy);
	mat3 roll = mat3(cr, sr, 0, -sr, cr, 0, 0, 0, 1);
    mat3 rotation = pitch*yaw*roll;    
    return rotation;
}

///////////////////////////////////////////////////////////////////////////////////
// SDF's & other spatial query functions

float dot2(in vec3 v ) {return dot(v,v);}
float dot2( in vec2 v ) { return dot(v,v); }

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

// 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;
}

// custom implementation of SD tri which skews space to bend the sail
float sdSail( vec3 p, vec3 a, vec3 b, vec3 c, float time )
{
    vec3 ba = b - a; 
    vec3 cb = c - b; 
    vec3 ac = a - c; 
    vec3 nor = normalize(cross( ba, ac ));
    vec3 acc = normalize(cross(ba, nor));
    float wid = dot(cb, acc);
    float w = dot(p-b, acc);
    float x = clamp((w/wid), 0.0, 1.0);
    float def = clamp(0.04405594 + 7.036014*x - 14.29915*x*x + 7.222222*x*x*x, 0.0, 1.0);
    //float def = sin(x*3.14159);
    def *= (0.6+0.2*cos((x+p.z)*3.14159*2.0 + time*3.0));
    p+=def*nor*0.07;
    float certainty = 1.0 - def*0.5;

    vec3 pa = p - a;
    vec3 pb = p - b;
    vec3 pc = p - c;
    
    if (sign(dot(cross(ba,nor),pa)) + sign(dot(cross(cb,nor),pb)) + sign(dot(cross(ac,nor),pc))<2.0)
    {
        return certainty*sqrt(min( min(
            dot2(ba*clamp(dot(ba,pa)/dot2(ba),0.0,1.0)-pa),
            dot2(cb*clamp(dot(cb,pb)/dot2(cb),0.0,1.0)-pb) ),
            dot2(ac*clamp(dot(ac,pc)/dot2(ac),0.0,1.0)-pc) ));
        
    }
    else
    {
    	return certainty*sqrt(dot(nor,pa)*dot(nor,pa)/dot2(nor));    
    }
}

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

vec2 opS( vec2 d1, vec2 d2 ) 
{ 
    return (-d1.x>d2.x) ? vec2(-d1.x,d1.y):d2; 
}

vec2 opI( vec2 d1, vec2 d2 ) 
{ 
    return (d1.x>d2.x) ? d1:d2; 
}

vec2 opsU( vec2 d1, vec2 d2, float k ) 
{
    float h = clamp( 0.5 + 0.5*(d2.x-d1.x)/k, 0.0, 1.0 );
    return vec2(mix( d2.x, d1.x, h ) - k*h*(1.0-h), d2.y); 
}

///////////////////////////////////////////////////////////////////////////////////

float hashfloat( uint n ) 
{
    // integer hash copied from Hugo Elias
	n = (n << 13U) ^ n;
    n = n * (n * n * 15731U + 789221U) + 1376312589U;
    return float( n & uvec3(0x7fffffffU))/float(0x7fffffff);
}

///////////////////////////////////////////////////////////////////////////////////

float planeIntersect( in vec3 ro, in vec3 rd, in vec4 p )
{
    return -(dot(ro,p.xyz)+p.w)/dot(rd,p.xyz);
}

vec2 sphereIntersect( in vec3 ro, in vec3 rd, in vec3 ce, float ra )
{
    vec3 oc = ro - ce;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - ra*ra;
    float h = b*b - c;
    if( h<0.0 ) return vec2(-1.0); // no intersection
    h = sqrt( h );
    return vec2( -b-h, -b+h );
}

///////////////////////////////////////////////////////////////////////////////////
    
float blerp(float x, float y0, float y1, float y2, float y3) {
	float a = y3 - y2 - y0 + y1;
	float b = y0 - y1 - a;
	float c = y2 - y0;
	float d = y1;
	return a * x * x * x + b * x * x + c * x + d;
}

float rand(vec2 co){
  return fract(sin(dot(co.xy ,vec2(12.9898,78.233))) * 43758.5453);
}

float perlin(float x, float h) {
	float a = floor(x);
	return blerp(mod(x, 1.0),
		rand(vec2(a-1.0, h)), rand(vec2(a-0.0, h)),
		rand(vec2(a+1.0, h)), rand(vec2(a+2.0, h)));
}

///////////////////////////////////////////////////////////////////////////////////

float sinbetween(float angle, float min, float max)
{
	return min + (1.0+sin(angle))*0.5*(max-min);
}

float perlinbetween(float a, float b, float min, float max)
{
	return min + (perlin(a,b))*(max-min);
}

///////////////////////////////////////////////////////////////////////////////////
// Non physical based atmospheric scattering made by robobo1221

const float pi = 3.14159265359;
const float invPi = 1.0 / pi;

const float zenithOffset = -0.02;
const float multiScatterPhase = 0.1;
const float density = 0.7;

const float anisotropicIntensity = 0.0; //Higher numbers result in more anisotropic scattering

const vec3 skyColor = vec3(0.39, 0.57, 1.0) * (1.0 + anisotropicIntensity); //Make sure one of the conponents is never 0.0

#define smooth(x) x*x*(3.0-2.0*x)
#define zenithDensity(x) density / pow(max(x - zenithOffset, 0.35e-2), 0.75)

vec3 getSkyAbsorption(vec3 x, float y){
	
	vec3 absorption = x * -y;
	     absorption = exp2(absorption) * 2.0;
	
	return absorption;
}

float getSunPoint(vec2 p, vec2 lp){
	return smoothstep(0.03, 0.026, distance(p, lp)) * 50.0;
}

float getRayleigMultiplier(vec2 p, vec2 lp){
	return 1.0 + pow(1.0 - clamp(distance(p, lp), 0.0, 1.0), 2.0) * pi * 0.5;
}

float getMie(vec2 p, vec2 lp){
	float disk = clamp(1.0 - pow(distance(p, lp), 0.1), 0.0, 1.0);
	
	return disk*disk*(3.0 - 2.0 * disk) * 2.0 * pi;
}

vec3 jodieReinhardTonemap(vec3 c){
    float l = dot(c, vec3(0.2126, 0.7152, 0.0722));
    vec3 tc = c / (c + 1.0);
    return mix(c / (l + 1.0), tc, tc);
}

vec3 getAtmosphericScattering(vec3 dir, vec3 sundir)
{
    // convert vec3's into the screen coordinates that robobo1221's code expected
    vec2 p;
    vec2 lp;
    sundir = sundir;
    sundir.y = -sundir.y;
    p.x = 0.0;
    p.y = atan(dir.y, length(dir.xz))/(pi*0.5);    
    lp.y = atan(sundir.y, length(sundir.xz))/(pi*0.5);
		
    vec3 dirflat = dir;
    dirflat.y=0.0;
    dirflat=normalize(dirflat);
    vec3 sundirflat = sundir;
    sundirflat.y=0.0;
    sundirflat=normalize(sundirflat);
    lp.x = 0.0 + acos(clamp(dot(-sundirflat, dirflat), -1.0, 1.0)) /(pi*0.5);
        
	float zenith = zenithDensity(p.y);
	float sunPointDistMult =  clamp(length(max(lp.y + multiScatterPhase - zenithOffset, 0.0)), 0.0, 1.0);
	
	float rayleighMult = getRayleigMultiplier(p, lp);
	
	vec3 absorption = getSkyAbsorption(skyColor, zenith);
    vec3 sunAbsorption = getSkyAbsorption(skyColor, zenithDensity(lp.y + multiScatterPhase));
	vec3 sky = skyColor * zenith * rayleighMult;
	vec3 sun = getSunPoint(p, lp) * absorption;
	vec3 mie = getMie(p, lp) * sunAbsorption;
	
	vec3 totalSky = mix(sky * absorption, sky / (sky + 0.5), sunPointDistMult);
         totalSky += sun + mie;
	     totalSky *= sunAbsorption * 0.5 + 0.5 * length(sunAbsorption);

	totalSky = jodieReinhardTonemap(totalSky);
	totalSky = pow(totalSky, vec3(2.2)); //Back to linear	
	return totalSky;
}

///////////////////////////////////////////////////////////////////////////////////
// Dave Hoskins's amazing tileable water caustic

#define TAU pi*2.0
#define CASUTIC_ITER 3

vec3 caustic(vec2 p, float time)
{
    p = mod(p*TAU, TAU)-250.0;
	vec2 i = vec2(p);
	float c = 1.0;
	float inten = .005;

	for (int n = 0; n < CASUTIC_ITER; n++) 
	{
		float t = time * (1.0 - (3.5 / float(n+1)));
		i = p + vec2(cos(t - i.x) + sin(t + i.y), sin(t - i.y) + cos(t + i.x));
		c += 1.0/length(vec2(p.x / (sin(i.x+t)/inten),p.y / (cos(i.y+t)/inten)));
	}
	c /= float(CASUTIC_ITER);
	c = 1.17-pow(c, 1.4);
	vec3 colour = vec3(pow(abs(c), 3.0));
    colour = clamp(colour + vec3(0.0, 0.35, 0.5), 0.0, 1.0);
    return colour;
}

///////////////////////////////////////////////////////////////////////////////////

mat4 brightnessMatrix( float brightness )
{
    return mat4( 1, 0, 0, 0,
                 0, 1, 0, 0,
                 0, 0, 1, 0,
                 brightness, brightness, brightness, 1 );
}

mat4 contrastMatrix( float contrast )
{
	float t = ( 1.0 - contrast ) / 2.0;
    
    return mat4( contrast, 0, 0, 0,
                 0, contrast, 0, 0,
                 0, 0, contrast, 0,
                 t, t, t, 1 );

}

mat4 saturationMatrix( float saturation )
{
    vec3 luminance = vec3( 0.3086, 0.6094, 0.0820 );
    
    float oneMinusSat = 1.0 - saturation;
    
    vec3 red = vec3( luminance.x * oneMinusSat );
    red+= vec3( saturation, 0, 0 );
    
    vec3 green = vec3( luminance.y * oneMinusSat );
    green += vec3( 0, saturation, 0 );
    
    vec3 blue = vec3( luminance.z * oneMinusSat );
    blue += vec3( 0, 0, saturation );
    
    return mat4( red,     0,
                 green,   0,
                 blue,    0,
                 0, 0, 0, 1 );
}
///////////////////////////////////////////////////////////////////////////////////
// The beauty of the sea ////////////////////////////////////////////////// peet //
///////////////////////////////////////////////////////////////////////////////////
// Put together from everything i've learned from all you awesome guys. ///////////
// Where there was an option between realism and aestehtics, I chose aesthetics ///
// Hope you like it :) ////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

#define DE_ITERATIONS 50		

///////////////////////////////////////////////////////////////////////////////////

#define WAVEHEIGHT	0.55

///////////////////////////////////////////////////////////////////////////////////

#define MSAA_X 1.0				
#define MSAA_Y 1.0				

///////////////////////////////////////////////////////////////////////////////////

Light g_light = Light(vec3(0.0, 0.0, 1.0), vec3(100.0, 100.0, 100.0)*0.25);
vec3 g_camera = vec3(0.0);

#define BASIC 0
#define ERROR 1

Material g_basic = Material(vec3(0.5, 0.5, 0.5), 0.5, 0.01);
Material g_error = Material(vec3(1.0, 0.0, 1.0), 1.0, 0.0);

Material g_mats[10];

const float g_movespeed = 3.0;

///////////////////////////////////////////////////////////////////////////////////

vec3 tritex(sampler2D tex, vec3 p, vec3 n)
{
 	return  (texture(tex,p.xy).rgb*n.z*n.z
            +texture(tex,p.zy).rgb*n.x*n.x
            +texture(tex,p.xz).rgb*n.y*n.y);
}

///////////////////////////////////////////////////////////////////////////////////
// IQ's texture based noise function (37,17) FTW!

float noise( in vec3 x )
{
    x *= 5.0;
    vec3 p = floor(x);
    vec3 f = fract(x);
	f = f*f*(3.0-2.0*f);
	vec2 uv = (p.xy+vec2(37.0,17.0)*p.z) + f.xy;
	vec2 rg = textureLod( iChannel0, (uv+ 0.5)/256.0, 0. ).yx;
	return mix( rg.x, rg.y, f.z );
}

///////////////////////////////////////////////////////////////////////////////////
// cheap environment backdrop - modified from bananaft's code

vec3 getfogcolour()
{    
    return getAtmosphericScattering(-g_light.dir, g_light.dir);
}

vec3 getlightcolour()
{    
    return getAtmosphericScattering(normalize(-g_light.dir - YAXIS*0.2), g_light.dir);
}

float getcloudtexture(vec3 p)
{   
    float c = 0.0;
    float weightsum = 0.0;
    float weight = 0.6;
    vec3 scale = vec3(0.0005, 0.0, 0.002);
    
    p+=vec3(0.1, 0.0, 0.2)*iTime*100.0;
    
    for (int i=0; i<4; i++)
    {
        weightsum+=weight;
        //c += texture(iChannel1, p.xz*scale.z*1.0).r*weight;
    	c += (noise(p*scale)*weight + noise(p*scale + YAXIS*1.0)*weight)*0.5;
        scale *= 1.9;
        weight *= 0.45;
    }
    c/=weightsum;
    return c;
}    
    
vec3 environment(vec3 origin, vec3 dir)
{
    float t = planeIntersect(origin, dir, vec4(0.0, -1.0, 0.0, 800.0));

    float c = getcloudtexture(origin+t*dir);
//    float delta = -(getcloudtexture(origin+(t+100.0)*dir) - c);
    float delta = getcloudtexture(origin+t*dir + 100.0*g_light.dir) - c;
    
    float alpha = clamp((1.0 - t/4000.0), 0.0, 1.0);
    if (t<0.0)
        alpha = 0.0;
    
    float cloudthreshold = 0.4;
    alpha *= clamp((c-cloudthreshold)/(0.75), 0.0, 1.0);
    
    vec3 col=vec3(c*40.0)*getlightcolour() + vec3(clamp(-g_light.dir.y, 0.0, 10.0))*80.0;
    //col += getlightcolour()*clamp((50.0+delta*8000.0), 0.0, 14000.5);
                                               
	vec3 sundir = g_light.dir; 
    vec3 env = getAtmosphericScattering(dir, sundir)*180.0;
    return mix(env, col, alpha);

}

///////////////////////////////////////////////////////////////////////////////////
// Fresnel

float fresnelR(vec3 d, vec3 n)
{
    float a = clamp(1.0-dot(n,-d), 0.0, 1.0);
    return clamp(exp((5.0*a)-5.0), 0.0, 1.0);
}
float fresnelT(vec3 d, vec3 n)
{
	return 1.0 - fresnelR(d,n);    
}

///////////////////////////////////////////////////////////////////////////////////
// DISTANCE ESTIMATION ////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

float anim(in vec3 pos)
{	// return the x offset at a given position(based on z)
    return sin(pos.z*3.0 + (iTime-pos.z*0.1)*8.0)*0.4;
}

vec2 sceneSDF( vec3 pos )
{
    vec2 res = vec2( TOO_FAR, ERROR );
    vec2 res2 = vec2( TOO_FAR, ERROR );
    
    pos.z -= iTime*g_movespeed;
        
    for (float i=0.0; i<=32.0; i++)
    {
        float rad = 0.06 + sin(3.14159*i/32.0)*0.02;
        vec3 p=vec3(0.0,-0.2,6.5 - i*3.0/32.0);
        p.x=anim(p);
        res = opsU(res, vec2(sdSphere((pos-p)*vec3(1.0, 0.4, 1.0), rad), 8.0), rad*1.75);
    }    
    return res;
}

///////////////////////////////////////////////////////////////////////////////////
// SDF system (Distance Estimate)

Result resultsDE(vec3 p)
{
    Result result = Result(vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0), TOO_FAR, TOO_FAR, 0.0, g_error);

    // normal generation
    result.normal = vec3(0.0);
    for( int i=min(iFrame,0); i<4; i++ )
    {
        vec3 e = 0.5773*(2.0*vec3((((i+3)>>1)&1), ((i>>1)&1), (i&1))-1.0);
        result.normal += e*sceneSDF(p+EPSILON*e).x;
    }
    result.normal = normalize(result.normal);        
    
    vec2 res = vec2(TOO_FAR, ERROR);    
    res = sceneSDF(p);
            
    result.mat=g_basic;
        
    result.pos = p;
    result.dist = res.x;
    return result;
}

///////////////////////////////////////////////////////////////////////////////////
// raymarch solid (Distance Estimate - based on value returned by sceneSDF())

Result raymarchDE(Ray ray, float distfactor, float maxtravel)
{
	float travelled=0.0;
    for (int i=0; i<DE_ITERATIONS; i++)
    {
    	vec2 res = sceneSDF(ray.pos);
                        
        ray.pos += res.x*ray.dir*distfactor; 
        travelled += res.x*distfactor;
        
        if (travelled>maxtravel)
            break;
    }     
    
    Result result = resultsDE(ray.pos);
    result.travelled=travelled;
    return result;
}

///////////////////////////////////////////////////////////////////////////////////
// HEIGHTFIELD ////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

float waterfield(vec3 p)
{	
    const float wavespeed = 1.0;
    float time = iTime;
    float cosa=cos(0.79);
    float sina=sin(0.79);
    mat2 rot = mat2(cosa, sina, -sina, cosa);
    
    vec2 phase = vec2(p.x, p.z);
    float w = 0.0;
    float mag = 1.0;
    float magtotal = 0.0;
    for (int i=0; i<3; i++)
    {
	    phase.x += time * 1.73 * (mag*wavespeed);
	    phase.y += time * 1.51 * (mag*wavespeed);
        w += (sin(phase.x) + cos(phase.y*0.3)) * mag;
        phase=rot*phase;
        magtotal += (mag*2.0);
        mag *= 0.675;	// should be approx 0.5
        phase *= 1.5;   // should be approx 2.0     
    }
    w/=magtotal;	// -1 to 1
    
    float y = w*WAVEHEIGHT*clamp(-g_light.dir.y, 0.0, 1.0)+0.0;
    
    vec3 snaketail = vec3(0.0, 0.0, g_movespeed*iTime + 6.0);
    vec3 totail = snaketail - p;
    if (totail.z > 0.0 && abs(totail.x)<totail.z)
    {
        float wake = sin(length(totail)*10.0-iTime*20.0)*0.1;
        wake *= 1.0 - (totail.x*totail.x)/(totail.z*totail.z);//(1.0-abs(totail.x)/totail.z)*(1.0-abs(totail.x)/totail.z);
        wake *= clamp(totail.z*0.25, 0.0, 1.0)*clamp((6.0-totail.z)*0.1, 0.0, 1.0);
    	y+=wake*clamp(-g_light.dir.y, 0.0, 1.0);    
    }
    
    return y;
}

///////////////////////////////////////////////////////////////////////////////////

float heightfield(vec3 p)
{	
    p*=0.035;    // scale down terrain frequency    
    float mag = 1.0;
    float magtotal = 0.0;
    float h = 0.0;
	for (int i=0; i<5; i++)
    {
        magtotal += mag;
        h += noise(vec3(p.x, 0.0, p.z))*mag;
        mag *= 0.5;
        p *= 2.03;
    }
    h = h/magtotal; //0 to 1
    h = (h*h*h*h-0.2)*32.0;
    h *= -sign(h);
    return h-noise(p*0.25)-(2.0-noise(p*0.1)*4.0);
}

///////////////////////////////////////////////////////////////////////////////////

float terrain(vec3 p, bool water)
{
    if (water)
		return waterfield(p) / clamp((dot(g_camera-p, g_camera-p)*0.01), 1.0, 100000.0);    
    else
		return heightfield(p) - 0.2;    
}

///////////////////////////////////////////////////////////////////////////////////
// raymarch heightfield (based on value returned by heightfield())

Result raymarchHF(Ray ray, float mindist, float maxdist, float stepsize, bool water)
{
    Result result;
    float dist = mindist;
    float h = 1.0;
    float dh = 1.0;
    float lastdh = -1.0;
    float lastdist = 0.0;
    float fracstep = 0.5;
    
    float maxy = water?WAVEHEIGHT:0.0;
    dist = max(dist, planeIntersect(ray.pos, ray.dir, vec4(0.0, 1.0, 0.0, -maxy)));
        
    for (float i=0.0; i<90.0; i++)
    {
        vec3 pos = ray.pos + ray.dir*dist;
        h = terrain(pos, water);
        dh = pos.y - h;
        if (dh<(0.001*i) || dist>maxdist)
        {
            break;    
        }
        lastdh = dh;
        lastdist = dist;
        dist+=stepsize;
        stepsize+=(0.0004*i);
    }
    
    if (dist<maxdist)
    {
        fracstep = lastdh/(lastdh-dh);    
        dist = mix(lastdist, dist, fracstep);
    }
    result.dist = dist;
    result.travelled = dist;
    result.pos = ray.pos + ray.dir*dist; 
    const float dn = 0.1;	
    result.normal.x = terrain(result.pos, water) - terrain(result.pos+XAXIS*dn, water);
    result.normal.y = dn;
    result.normal.z = terrain(result.pos, water) - terrain(result.pos+ZAXIS*dn, water);
    result.normal = normalize(result.normal);
    result.mat = g_basic;
    return result;
}

///////////////////////////////////////////////////////////////////////////////////
// Heightfield integrator

vec3 integrateHF( Ray inputray )
{
    const float exposure = 1e-2;
    const float gamma = 2.2;
    const float intensity = 100.0;
	const float maxdist = 400.0;
            
    float IOR = 1.4;
    if (inputray.pos.y < waterfield(inputray.pos))
        IOR = 1.0/IOR;
        
    vec3 colour = vec3(0.0, 0.0, 0.0);
    
    Ray ray=inputray;            
    Result result = raymarchHF(ray, 1.0, maxdist, 0.01, true);

    if (result.dist>=maxdist)
    {
        colour = environment(ray.pos, ray.dir);
	    colour.xyz = vec3(pow(colour * exposure, vec3(1.0 / gamma)));    
	    return colour;    
    }

    float fR=fresnelR(ray.dir, result.normal);
    float fT=1.0 - fR;
    vec3 watercolour = 2.0*vec3(0.0, 0.3, 0.2)*clamp(-g_light.dir.y, 0.0, 1.0);
    
    {
        vec3 refracted = refract(ray.dir, result.normal, 1.0/IOR);
        // underwater terrain intersection
        Ray groundray;
        groundray.pos = result.pos;
        groundray.dir = refracted;
        
	    Result snakeresult = raymarchDE(groundray, 1.0, 3.0);
        if (snakeresult.travelled < 3.0)
        {
            vec3 snakepos = groundray.pos + groundray.dir*(snakeresult.travelled);
            vec3 snakeposl = snakepos - vec3(0.0, 0.0, iTime*g_movespeed);
            vec3 snakecol = vec3(200.0);
            float f = mod(snakeposl.z, 0.15);
            if (f<(((snakepos.y+0.2)*0.2+0.0375)))
                snakecol = vec3(0.0);
            if (snakeposl.z>6.3)
            {
	            if (snakeposl.y>-0.05)
                    snakecol = vec3(0.0);
    			else
					snakecol = vec3(1.0, 0.8, 0.2)*100.0;
            }
            
            Result sr = resultsDE(snakepos);        
            snakecol*=0.25+vec3(tritex(iChannel1, snakeposl, sr.normal).r)*0.75;        
            snakecol = mix(snakecol, watercolour, clamp(snakeresult.travelled*snakeresult.travelled*0.01, 0.3, 1.0));        
            colour = snakecol*clamp(dot(sr.normal, -g_light.dir), 0.0, 1.0)*fT*g_light.colour;        
        }
        else
        {                    
            Result groundresult = raymarchHF(groundray, 0.1, 100.0, 0.1, false);
            vec3 ground = groundresult.pos;
            float grounddist = length(result.pos-ground);

            // texturing
            vec3 groundcol = texture(iChannel1, ground.xz*0.3).xyz;
            groundcol*=(groundcol+vec3(0.6));

            // caustics (Thanks Dave!)
            vec2 h = mod(groundresult.pos.xz, 10.0)*groundresult.pos.y*0.1;            
            vec3 cc = caustic(h, iTime);        
            groundcol*=0.5 + dot(cc,cc)*0.5;
            groundcol+=dot(cc,cc)*0.5*clamp(-g_light.dir.y, 0.0, 1.0);

            groundcol *= clamp(dot(groundresult.normal, -g_light.dir), 0.0, 1.0);
            groundcol = mix(groundcol, watercolour, clamp(grounddist*grounddist*0.01, 0.0, 1.0));        
            groundcol = groundcol*groundcol;
            colour = mix(colour, groundcol*40.0, fT)*g_light.colour;  
        }
    }
    
    Ray reflectray;
    reflectray.pos = result.pos + result.normal*0.02f;
    reflectray.dir = reflect(ray.dir, result.normal);

    vec3 ambient = environment(reflectray.pos, reflectray.dir);
    colour += fR*(ambient);                                     

    colour.xyz = mix(colour.xyz, environment(inputray.pos, inputray.dir), result.dist/maxdist);    
    colour.xyz = vec3(pow(colour * exposure, vec3(1.0 / gamma)));    
    return colour;    
}

///////////////////////////////////////////////////////////////////////////////////
// 'Camera'

Ray CreateRay(vec2 uv)
{
    // moveable cam
    vec3 p0, p1;
    p0=vec3(0.0, 2.0, 0.0);
    p1=vec3(0.0, 0.0, 3.0);

    vec2 vmouse;    
    if (iMouse.x<20.0)
    {
        // auto!
        const float introtime=17.3;
        const float beat = 4.02;
        if (iTime<introtime)
        {
        	vmouse = vec2(0.5+sin(iTime*0.3)*0.1, 0.25+sin(iTime*0.2)*0.15);
        }
        else
        {
        	float t = mod(iTime-introtime, beat);
            float s = floor((iTime-introtime) / beat);
            vec2 v1 = vec2(perlin(s, 1.0), 0.1+sqrt(perlin(s, 2.0))); 
            vec2 v2 = vec2(perlin(s, 3.0), 0.1+perlin(s, 4.0)*perlin(s, 4.0));
            
            if (mod(s,10.0)>=6.0 && mod(s,10.0)<=8.0)
            {
            	v1.y = 0.1+v1.y*0.25;    
            	v2.y = 0.1+v2.y*0.25;    
            }
            
            t=clamp(t/2.5, 0.0, 1.0);                   
            vmouse=mix(v1, v2, t*t*(3.0-2.0*t));            
        }
    }    
    else
    {
        vmouse = vec2(iMouse.x/iResolution.x, iMouse.y/iResolution.y);
    }    
	p0=vec3(0.0 + 16.0*vmouse.x - 8.0, 16.0*vmouse.y, 0.0);
    
    p0.z += iTime*g_movespeed;        
    p1.z += iTime*g_movespeed;    
                        
    Ray ray;
    float roll = 0.0;
    vec3 dir = (p1-p0) + vec3(0.0, 0.0, 4.0);
    dir = normalize(dir);
    vec3 up = vec3(dir.x*roll, 1.0, 0.0);
    up = normalize(up);
    vec3 right = cross(dir, up);
    right = normalize(right);
    up = cross(right, dir);
    up = normalize(up);
    
    ray.pos = vec3(0.0, 0.0, 0.0) + p0;
    ray.dir = dir*0.8 + up*uv.y + right*uv.x;
    ray.dir = normalize(ray.dir); 
    return ray;
}

///////////////////////////////////////////////////////////////////////////////////
// main loop, iterate over the pixels, doing MSAA

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{               
    fragColor = vec4(0.0, 0.0, 0.0, 1.0);
    float factor = 1.0/(MSAA_X*MSAA_Y);
             
    // materials
    g_mats[BASIC]=g_basic;
    g_mats[ERROR]=g_error;
    
    vec3 sundir = normalize(vec3(0.0, -0.87+sin((iTime*0.5+5.0)*0.3), -1.0));
    g_light.dir = sundir;
    g_light.colour = getlightcolour();
    
    for (float x=0.0; x<MSAA_X; x++)
    {
        for (float y=0.0; y<MSAA_Y; y++)
        {
            vec2 uv = fragCoord.xy / iResolution.xy * 2.0 - 1.0;
            uv.y *= iResolution.y / iResolution.x;

            uv.x += (1.0/(iResolution.x*MSAA_X))*x;
            uv.y += (1.0/(iResolution.y*MSAA_Y))*y;
            
    		Ray ray=CreateRay(uv);
		    g_camera = ray.pos;
                                    
#ifdef DITHER            
            float dither = hashfloat(uint(fragCoord.x+iResolution.x*fragCoord.y)+uint(iResolution.x*iResolution.y)*uint(iFrame));//Updated with iFrame dimension    
            ray.pos += ray.dir*dither*0.1;
#endif //DITHER            
            
            const float outofboundsdist = 3.0;
            const float maxvaliddist = 0.02;
            fragColor.xyz += integrateHF(ray)*factor;              
       	}        
    }    
    
    vec2 uv = fragCoord.xy / iResolution.xy * 2.0 - 1.0;
    uv.y *= iResolution.y / iResolution.x;    
    fragColor = brightnessMatrix( 0.05 ) * contrastMatrix( 1.5 ) * 
        saturationMatrix( 1.2 ) * fragColor;
    fragColor.xyz -= (uv.y*uv.y + uv.x*uv.x)*0.2;             
}

///////////////////////////////////////////////////////////////////////////////////