const vec3 sundir = normalize( vec3(-1.,.5,-1.) );
const vec3 suncol = vec3(1.,.8,.5);

mat2 m2 = mat2(0.8,  0.6, -0.6,  0.8);
mat3 m3 = mat3( 0.00,  0.80,  0.60,
              -0.80,  0.36, -0.48,
              -0.60, -0.48,  0.64 );

bool sharp=true;

mat2 r2D(float a)
{
	float si = sin(a);
	float co = cos(a);
	return mat2(si, co, -co, si);
}

float fbm(vec2 p)//ground
{
    float s=.3;
    mat2 m = r2D(1.);
    vec2 r= vec2(0.);
    for(int i=0;i<7;i++)r+=sin(p), p=m*(p*s+.2*r+.1),s*=1.6;
    return r.x+r.y;
}

float fbm1(vec2 p)//grass
{
    float s=2.;
    mat2 m = r2D(1.);
    vec2 r= vec2(0.);
    for(int i=0;i<8;i++)r+=(sin(p.yx+.5))/s, p=m*p*s+cos(r.yx*s), s*=1.02;
    return (r.x+r.y)/s;
}


float fbm2(vec3 p, float t)//clouds
{
    float s=1., r= -4.3;
    p += 4.*cos(.07*p)+5.*sin(.05*p.zyx+1.);
    p*=.4;
    p.xz *=.4;
    vec3 n = vec3(1);
    for(int i=0;i<8;i++)
        p = m3*p.zxy,
        n+=.8*cos(p*s+.05*t),
        r+=abs(dot(sin(p*s+n.zxy-1.)/s,vec3(1.))),
        s*=1.9;

    return -r;
}


//reef/waves combined field
vec3 map(in vec3 p) {

	float s=2.,e,f,o;
    vec3 q=p,r=p;
    vec3 n = vec3(0);
	vec2 l = vec2(2.);
	for(e=f=p.y;s<5e1;s*=1.3)
            p.xz*=m2,
            n.xy*=m2,
            q=p*s+n,
            r=p*s+n,
            r.x+=iTime*.5*s,
            e+=.25*abs(dot(sin(r.xz*.1)/s,.8*l)),
            f+=.2+.2*(dot(sin(q.xz*.5)/s,l)),
            n-=cos(q);
    e+=1.7;
    return vec3(min(e,f),f,e);
}

vec3 normalRocks(in vec3 p)
{
	const vec2 e = vec2(0.004, 0.0);
	return normalize(vec3(
		map(p + e.xyy).y - map(p - e.xyy).y,
        .008,
		map(p + e.yyx).y - map(p - e.yyx).y
		));
}

vec3 normalSea(in vec3 p)
{
	const vec2 e = vec2(0.005, 0.0);
	return normalize(vec3(
		map(p + e.xyy).z - map(p - e.xyy).z,
        .02,
		map(p + e.yyx).z - map(p - e.yyx).z
		));
}

const float sc =1.;
const float low = 7.;
const float high = 12.;
const float rt = 6e3;

float f(vec3 p)//clouds
{
    p *= sc;
    p.y += rt;
    p.y = length(p)-rt;
    p.z +=.2*iTime;
    float d =fbm2(p,iTime)/sc;
    return d+exp(2.*(low-p.y))+exp(p.y-high);
}

vec2 marchsky( vec3 ro, vec3 rd)
{
    int st = 64;
    if(rd.y<0.001)return vec2(0.);
    float t=low/(.07+rd.y),h=(high-low)/(.1+rd.y)/float(st),d,c=1.,dt=1.,k=1.5,e,a=0.;
    vec3 p;
    for(int i=0;i<st;++i)
    {
        p=ro+rd*t;
        d=f(p);
        e=sharp ? smoothstep(.4, .7, .6*d) : 1.-exp(-d*d);
        c*=e;
        t+=k*h*e;
        if(t>400.)return vec2(0.);
        a+=((f(p+sundir*dt)-d)/dt+.5)*(1.-e)*c;
        k*=1.015;
    }

    return vec2(a*.7,(1.-c)*1.2);
}

vec3 sky( in vec3 ro, in vec3 rd )//modified from IQ clouds
{
    // background sky
    float sun = clamp( dot(sundir,rd), 0.0, 1.0 );
    vec3 col = vec3(0.6,0.6,0.78) - rd.y*0.5*vec3(1.0,0.4,.05) ;
    col += 0.4*suncol*pow( sun, 8.0 );
    col *=.9;
    // clouds
    vec2 res = marchsky( ro, rd);
    float k = res.x, c=res.y;
    if(c>.0)
       col *= 1.-.1*c,
       col += (.2+k)*c*suncol,
       col += vec3(0.2,0.08,0.04)*pow( sun, 3.0 )*c;
    return col;
}


vec3 march(in vec3 ro, in vec3 rd)
{
	const float maxd = 50.0;
	const float precis = 0.001;
    float h = 0.0;
    float t = 0.0;
    float dt = .2;
	float res = -1.0;
    for(int i = 0; i < 128; i++)
    {
        if(h < precis*t || t > maxd) break;
	    h = map(ro + rd * t).x;
        t += h*dt;
        dt *= 1.015;
    }
    if(t < maxd) res = t;
    return vec3(res,map(ro + rd * t).yz);
}

bool keypress(int key) {
    return texelFetch(iChannel1, ivec2(key,2),0).x != 0.0;
}


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


	vec3 col = vec3(0.);
    vec3 li = sundir;
    vec2 uv = fragCoord/iResolution.xy -.5;
    uv.x *= iResolution.x/iResolution.y;
    vec2 m = iMouse.xy/iResolution.xy;


    float zoom=1.;
    sharp = true;
    if(keypress(88))sharp = !sharp;
    vec3 ro=1.*vec3(sin(3.*m.x),0.*m.y+2.,cos(3.*m.x));
    vec3 rd=normalize(vec3(uv.xy,zoom));
    vec3 target=vec3(0,2.,0);

    vec3 w=normalize(target-ro);
    vec3 u=normalize(cross(w,vec3(0,1,0)));
    vec3 v=normalize(cross(w,-u));

    rd=mat3(u,v,w)*rd;


    vec3 a = march(ro, rd);
    float t = a.x;
    float dh = a.z-a.y;

    vec3 skyCol = sky(ro,rd);
    col = skyCol;

    if(t > 0.)
    {

        vec3 pos = ro + t * rd;
        float k=map(pos).z*1.+1.3;
        vec3 nor = normalRocks(pos);
        float r = max(dot(nor, li),0.1)/2.5;
        col =r*vec3(k*k, k, .8)+.2*exp(-50.*dh*dh)+.05;
        if(dh<0.1){
        	vec3 nor = normalSea(pos+dh*rd);
        	nor = reflect(rd, nor);
            col +=vec3(0.9,.2,.05)*dh*.5;
        	col += pow(max(dot(li, nor), 0.0), 5.0)*vec3(.7);
            uv.y*=-1.;
            col*=.5;
        	col +=.6* sky(ro,nor);

        }
	    col = .1+col;

	}
    float maxd = 50.;
    col = mix(col, skyCol, smoothstep(.6, .99, min(t, maxd)/maxd));


   	fragColor = vec4(col, 1.0);
}