#define R iResolution.xy
#define T ((iTime+90.)*2.1)
#define ZERO min(iFrame, 0)
#define NEAR 0.003
#define FAR 1600.0
#define MAX_STEPS 40
#define M_PI 3.1415
#define TAU (M_PI*2.0)

float luma(vec3 color) { return dot(color, vec3(0.299, 0.587, 0.114)); }
vec3 aces(vec3 x) {return x*(2.51*x + .03) / (x*(2.43*x + .59) + .14); }

float smod(float x, float y) {
    float k = abs(0.5-fract(x/y));
    return smoothstep(0.0, 1.0, k);
}

#define NOISE(p, seed) (textureLod(iChannel3, ((p+0.289128) + (seed * 156.0))/256.0, 0.1).rgb)

vec3 noise(in vec2 p, in float seed) {
    vec2 id = floor(p);
    vec2 lv = fract(p);
    lv = lv * lv * (3.0 - 2.0 * lv);
    return mix(
        mix(NOISE(id, seed), NOISE(id+vec2(1, 0), seed), lv.x),
        mix(NOISE(id+vec2(0, 1), seed), NOISE(id+vec2(1, 1), seed), lv.x),
        lv.y
    );
}

vec3 noise(in vec2 p, in float seed, in float freq, const in int octaves) {
    float div = 0.0;
    float amp = 1.0;
    vec3 n = vec3(0.0);
    
    for (int i = ZERO; i < octaves; i++) {
        n += amp * noise(p*freq, seed);
        div += amp;
        amp /= 2.0;
        freq *= 2.0;
    }
    
    return n / div;
}

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

#define NEG(x) ((x)*2.0-1.0)

float sdfWater(in vec3 p) {
    if (p.y > 10.23) return p.y;
    float t = T*0.15;
    float freq = 0.59;
    
    vec2 offset = vec2(cos(t*2.0), sin(t*2.0))*1.3;
    vec2 mp = vec2(smod(p.x, 200.0), smod(p.y, 200.0)) * 0.5;
    
    float h = 0.0;
    vec2 shift = vec2(sin(t), cos(t)) * 1.5;
    vec2 dir = vec2(1, 0);
    float k1 = (0.5+(0.5*sin(mp.x*shift.y*0.1)));
    float k2 = (0.5+(0.5*sin(mp.y*shift.x*0.1)));
    float k = (k1 + k2) * 0.5;
    dir = mix(dir, vec2(0, 1), k*0.25);
    
    dir *= rot(t);
    vec3 n0 = noise(p.xz + (shift * 100.), 0.00001231, 0.009, 3);
    shift *= (0.25+(0.75*dot(normalize(shift), dir)))*2.;
    shift += 0.15*offset;
    h += NEG(n0.x * n0.y + n0.z);
    shift = mix(shift, (0.06*offset)+dir*vec2(sin(t+h-n0.x), cos(t+h-n0.y)), 0.5*smoothstep(0.4, 0.7, n0.x));
    vec3 n1 = noise(p.zx + (shift * 50.0), 11.982715, 0.01, 4);
    h += NEG(n1.x * n1.y + n1.z);
    shift = mix(shift, dir*vec2(sin(t+h-n1.x), cos(t+h-n1.y)), 0.5*smoothstep(0.4, 0.7, n1.z));
    shift -= offset*0.33;
    vec3 n3 = noise(p.xz - (shift * 30.0), 55.555315, 0.06, 3);
    h += NEG(n3.x) * smoothstep(0.4, 0.7, (n1.y+(n0.z*0.5))*0.6);
    vec3 n4 = noise((offset*0.03)+p.zx + (shift * 30.0), 201.0928182, 0.08, 6);
    float b = 0.15;
    float ib = 1.0 - b;
    h += NEG(n4.x * n4.y + n4.z) * (b + (ib * smoothstep(0.4, 0.9, (n0.z+n3.y+n0.z)*0.333333)));
    
    h *= 0.3333333;
    
    h *= 30.;
    return ((p.y+h)/1.1)-4.;
}
float ground(in vec3 p) {
    float h = 0.0;
    
    
    vec3 n0 = noise(p.xz*1.9, 1.119281, 0.003, 6);
    h += n0.x*(10.0+(n0.y + n0.z));
    h += (h*h);
    
    vec3 n1 = noise(p.zx, 10.9281812, 0.02, 6);
    h += (n1.x * n1.y + n1.z)*30. * smoothstep(0.4, 0.7, n0.y);
    
    vec3 n2 = noise(p.xz, 166.03989182, 0.003, 3);
    
    float fg = n2.x * n2.x;
    float g = ((fg*fg) + (n0.z * n0.y + n0.z * n0.y)) * 40.0;
    
    h -= g*1.1111;
    h += smoothstep(0.015, 0.5, n2.y)*100.*n2.y*n2.y*max(0.0, 1.0-fg);

    
    return h/1.05;
}
float sdf(in vec3 p, bool water) {
    if (p.y < -200.) return p.y;
    if (water) return sdfWater(p);
    if (p.y > 16.2) return p.y;

    float h = ground(p);
    
    return ((p.y+h)/1.1)+((p.y+4.0)*0.4);
}

struct Data {
    vec3 p;
    vec3 n;
    float d;
};

bool march(in vec3 ro, in vec3 rd, inout Data data, bool water, in float far) {
    data.d = FAR;
    float d = 0.0;
    for (int i = ZERO; i < MAX_STEPS; i++) {
        vec3 p = ro + rd * d;
        float next = sdf(p, water);
        
        d += next;
        if (abs(next) <= (NEAR * (1.0 + (abs(d) * 2.5)))) break;
        if (d >= far) return false;
    }
    
    d = abs(d);
    vec3 p = ro + rd * d;
    vec2 e = vec2(water ? 0.09 : 0.05, 0.0);
    data.n = normalize(sdf(p, water) - vec3(
        sdf(p - e.xyy, water),
        sdf(p - e.yxy, water),
        sdf(p - e.yyx, water)
    ));
    data.d = d;
    data.p = p;
    return true;
}

vec3 getSky(in vec3 rd) {
    float dotup = max(0.0, dot(rd, vec3(0, 1, 0)));
    vec3 col = vec3(0.22, 0.66, 0.72);
    col = col * col;
    col = pow(col, vec3(1.0 + dotup*2.0));
    return col;
}

vec3 blit(in vec3 ro, in vec3 rd, in Data data, in vec3 L, in vec3 lcol, in vec3 alb, float spf) {
    vec3 col = vec3(0.0);
    vec3 N = data.n;
    vec3 ref = reflect(L, N);
    float VdotR = max(0.0, dot(rd, ref));
    float spec = pow(VdotR, 128.0) * spf;
    vec3 diffuse = alb / M_PI;
    float NdotL = max(0.0, dot(N, L));
    vec3 att = NdotL * NdotL * diffuse * lcol;
    col += (att + spec);
    
    
    vec3 ref2 = reflect(rd, N);
    

    vec3 env = getSky(ref2);
    
    col += env * col;
    return col;
}


vec3 tt(in sampler2D samp, in vec2 uv, in vec3 N, in float d) {
    vec3 n1 = noise(uv, 0.329812, 2.5, 6);
    vec3 n2 = noise(uv.yx, 66.329812, 0.5, 6);
    vec3 n3 = noise(uv.xy+0.38281, 231.281, 0.4, 6);
    vec3 n4 = noise(uv.yx+13.3333, 125.38281, 0.1, 3);
    vec3 c1 = vec3(0.5, 0.3, 0.15);
    vec3 c2 = vec3(0.8, 0.33, 0.2);
    vec3 c3 = vec3(0.75, 0.39, 0.25);
    vec3 col = mix(c1, c2, n1.x);
    col = mix(col, c3, smoothstep(0.4, 0.7, n3.x));
    col = mix(col, clamp((col/M_PI)+0.19, 0.0, 1.0), 0.5*smoothstep(0.4, 0.7, n4.y));
    float cc = 1.0-smoothstep(0.01, 0.1, abs(n3.z*2.0-1.0));
    
    col += 0.5*(c2+col*0.8)*cc*smoothstep(0.4, 0.7, n4.z);
    return col;
}
vec3 render(in vec3 ro, in vec3 rd, inout float depth) {
    depth = FAR;
    vec3 col = vec3(0.0);
    Data data1 = Data(vec3(0.0), vec3(0.0), FAR);
    Data data2 = Data(vec3(0.0), vec3(0.0), FAR);
    Data data3 = Data(vec3(0.0), vec3(0.0), FAR);
    vec3 L = normalize(vec3(1, 2, 3));
    vec3 lcol = vec3(0.8, 0.7, 0.65);

    vec3 tmp = vec3(0.0);
  
    // ground over water
    if (march(ro, rd, data3, false, FAR)) {
      vec3 alb = tt(iChannel2, data3.p.xz, data3.n, data3.d);
      tmp = blit(ro, rd, data3, L, lcol, alb, 0.01);
      depth = (data3.d / FAR);
    } else {
       depth = FAR;
       return getSky(rd);
    }
   
    // water
    if (march(ro, rd, data1, true, FAR)) {
        vec3 alb = vec3(0.21, 0.66, 0.75);
        if (data1.d-0.001 < data3.d) {
            col += blit(ro, rd, data1, L, lcol, 0.22222*((alb*alb*alb)+0.5), 1.0);
            depth  = data1.d / FAR;
        }
        
        if (data3.d > data1.d-0.001) {
            float ior = 1.33;
            vec3 rd2 = refract(rd, data1.n, 1.0/ior);
            vec3 ro2 = data1.p-data1.n*0.001;
            
            
            // ground under water
            if (march(ro2, rd2, data2, false, FAR)) {


                float dd = distance(data1.p, data2.p) / FAR;

                vec3 alb = tt(iChannel2, data2.p.xz, data2.n, data2.d);
                col +=  blit(ro2, rd2, data2, L, lcol, alb, 0.1) / (1.0 + dd*6.0+(data2.d/FAR));

                depth = max(depth,  data2.d / FAR);

                float kj = distance(data2.p, data3.p);
                // foam
                col += smoothstep(0.04, 0.06, depth)*smoothstep(3.0, -4.0, kj) * noise(data3.p.xz+vec2(sin(T), cos(T)), 0.23918, 2.0, 1).x;

            }
        } else {
            float dd = distance(data1.p, data3.p) / FAR;
            col += tmp /  (1.0 + (dd*6.0));
        }
    }
  
 
    return col;
}
vec3 nextRo(in vec3 ro, in float t) {
    ro.y += (35.) + 10.0*(0.5*(0.5+sin(t)));
    ro.z += t*25.;
    ro.x += cos((t-11.2989122)*0.15)*90.;
    return ro;
}
void mainImage( out vec4 O, in vec2 fc )
{   vec3 col = vec3(0.0);
    vec2 uv = (fc - 0.5 * R.xy)/R.y;
    
    vec3 ro = vec3(0, 0.0, -1.0);
    vec3 rd = normalize(vec3(uv.xy, mix(1.0, 0.6, 0.33*(0.5+(0.5*cos(0.33*(T+33.329281)))))));
    
    vec4 m = vec4((iMouse.xy - 0.5 * R.xy)/R.y, iMouse.zw);
    
    if (m.z > 0.001) {
      rd.yz *= rot(m.y * TAU);
      rd.xz *= rot(m.x * TAU);
    } else {
        rd.yz *= rot(radians(-(15. + (0.5*(0.5+cos(T+1.82818))))));
        rd.xz *= rot(0.5*sin((T+11.023321)*0.25));
    }
    
    ro = nextRo(ro, T);
    float g1 = ground(ro);
    float g2 = ground(nextRo(ro, T+0.015));
    float g3 = ground(nextRo(ro, T+0.05));
    float g4 = ground(nextRo(ro, T+0.08));
    
    float g = (g1 + g2 + g3 + g4) * 0.25;
    
    ro.y += 3.0*(1.0 - 0.85*smoothstep(0.0, 30.0, distance(vec3(0, ro.y, 0), vec3(0, g, 0))));
    
    float depth = FAR;
    col = render(ro, rd, depth);
    float b = 0.45;
    float ib = 1.0 - b;
    float dotup = max(0.0, dot(rd, vec3(0, 1, 0)));
    float d = smoothstep(0.1, 0.9, depth);
    col += (d*d)* max(0.0, 1.0 - smoothstep(0.0, 0.44, dotup));
    
    float l = luma(col);
  
    col += (l*l + (col*col));
    col += (col*col*col*col);
    
    col *= (1.0+col*4.);
    
    col /= (1.0+max(col-0.25, 0.0));
    
    col = aces(col);
    
    col = pow(col, vec3(1.0 / 2.2));
    O = vec4(col,1.0);
}