float almostIdentity(float x, float m, float n) {
	if (x > m) return x;
	float t = x / m;
	return (((2.0 * n - m) * t + (2.0 * m - 3.0 * n)) * t * t) + n;
}

float getWave(sampler2D tex, vec2 pos) {
    pos -= 0.5;
    ivec2 ip = ivec2(floor(pos));
    vec2 fp = fract(pos);

    float[4] samples;
    for (int i = 0; i < 4; ++i) {
        ivec2 o = ivec2(i % 2, i / 2);
        ivec2 p = (ip + o) % int(SIZE);
        samples[i] = texelFetch(tex, p, 0).x ;
    }

    return mix(
        mix(samples[0], samples[1], fp.x),
        mix(samples[2], samples[3], fp.x),
        fp.y
    );
}

float getWaveBicubic(sampler2D sampler, vec2 coord) {
	vec2 res = vec2(SIZE, SIZE);

	coord = coord - 0.5;

	vec2 f = fract(coord);
	coord -= f;

	vec2 ff = f * f;
	vec4 w0;
	vec4 w1;
	w0.xz = 1 - f; w0.xz *= w0.xz * w0.xz;
	w1.yw = ff * f;
	w1.xz = 3 * w1.yw + 4 - 6 * ff;
	w0.yw = 6 - w1.xz - w1.yw - w0.xz;

	vec4 s = w0 + w1;
	vec4 c = coord.xxyy + vec2(-0.5, 1.5).xyxy + w1 / s;
	//c /= res.xxyy;

	vec2 m = s.xz / (s.xz + s.yw);
	return mix(
		mix(getWave(sampler, c.yw), getWave(sampler, c.xw), m.x),
		mix(getWave(sampler, c.yz), getWave(sampler, c.xz), m.x),
		m.y);
}

float calculateWaves(in vec2 position) {
    position *= SIZE*rcp(SIDELENGTH);
	float adjustment = 0.1;
    float waves = getWaveBicubic(gaux4, position) - adjustment;
	return waves;
}

vec3 waterNormal(vec2 position) {
	cFloat dist = 4e-3;

	vec2 diffs;
    diffs.x = calculateWaves(position + vec2( dist, -dist));
    diffs.y = calculateWaves(position + vec2(-dist, dist));
    diffs  -= calculateWaves(position + vec2(-dist, -dist));

	vec3 wavenormal = vec3(-2.0 * dist, 4.0 * dist * dist, -2.0 * (dist * dist + dist));
	wavenormal.xz *= diffs;

    return wavenormal;
}