precision highp float;
uniform vec2 resolution;
uniform vec2 mouse;
uniform float time;
uniform int frame;
out vec4 outColor;


#define PI 3.14159265359
#define ITER 256

//#define Use_Perlin
#define Use_Simplex

float gSceneTime;
float BPM = 120.0;
// ========= Hash ===========

// reference https://www.shadertoy.com/view/4djSRW
#define MOD3 vec3(443.8975,397.2973, 491.1871)

float hash12(vec2 p) {
	return fract(sin(dot(p, vec2(127.1, 311.7))) * 43758.5453);
}

vec2 hash22(vec2 p) {
	return fract(sin(vec2(dot(p, vec2(127.1, 311.7)), dot(p, vec2(269.5, 183.3)))) * 43758.5453);
}

float hash31(vec3 p3)
{
	p3  = fract(p3 * MOD3);
	p3 += dot(p3, p3.yzx + 19.19);
	return -1.0 + 2.0 * fract((p3.x + p3.y) * p3.z);
}

vec3 hash33(vec3 p3)
{
	p3 = fract(p3 * MOD3);
	p3 += dot(p3, p3.yxz+19.19);
	return -1.0 + 2.0 * fract(vec3((p3.x + p3.y)*p3.z, (p3.x+p3.z)*p3.y, (p3.y+p3.z)*p3.x));
}

// ========= Noise ===========

float noisePerline(vec3 p)
{
	vec3 pi = floor(p);
	vec3 pf = p - pi;
	
	vec3 w = pf * pf * (3.0 - 2.0 * pf);
	
	return 	mix(
				mix(
					mix(dot(pf - vec3(0, 0, 0), hash33(pi + vec3(0, 0, 0))), 
						dot(pf - vec3(1, 0, 0), hash33(pi + vec3(1, 0, 0))),
					   	w.x),
					mix(dot(pf - vec3(0, 0, 1), hash33(pi + vec3(0, 0, 1))), 
						dot(pf - vec3(1, 0, 1), hash33(pi + vec3(1, 0, 1))),
					   	w.x),
					w.z),
				mix(
					mix(dot(pf - vec3(0, 1, 0), hash33(pi + vec3(0, 1, 0))), 
						dot(pf - vec3(1, 1, 0), hash33(pi + vec3(1, 1, 0))),
					   	w.x),
				   	mix(dot(pf - vec3(0, 1, 1), hash33(pi + vec3(0, 1, 1))), 
						dot(pf - vec3(1, 1, 1), hash33(pi + vec3(1, 1, 1))),
					   	w.x),
					w.z),
				w.y);
}

float noiseSimplex(vec3 p)
{
	const float K1 = 0.333333333;
	const float K2 = 0.166666667;
	
	vec3 i = floor(p + (p.x + p.y + p.z) * K1);
	vec3 d0 = p - (i - (i.x + i.y + i.z) * K2);
	
	// thx nikita: https://www.shadertoy.com/view/XsX3zB
	vec3 e = step(vec3(0.0), d0 - d0.yzx);
	vec3 i1 = e * (1.0 - e.zxy);
	vec3 i2 = 1.0 - e.zxy * (1.0 - e);
	
	vec3 d1 = d0 - (i1 - 1.0 * K2);
	vec3 d2 = d0 - (i2 - 2.0 * K2);
	vec3 d3 = d0 - (1.0 - 3.0 * K2);
	
	vec4 h = max(0.6 - vec4(dot(d0, d0), dot(d1, d1), dot(d2, d2), dot(d3, d3)), 0.0);
	vec4 n = h * h * h * h * vec4(dot(d0, hash33(i)), dot(d1, hash33(i + i1)), dot(d2, hash33(i + i2)), dot(d3, hash33(i + 1.0)));
	
	return dot(vec4(31.316), n);
}

float noise(vec3 p) {
#ifdef Use_Perlin
	return noisePerline(p * 2.0);
#elif defined Use_Simplex
	return noiseSimplex(p);
#endif
	return 0.0;
}

// ======== Time Easing & Frequency ===========
// emphasisPower: 強調の度合い (2.0〜5.0程度で変化)
float getEasedTime(float curTime, float bpm, float emphasisPower) {
	float beatDuration = 60.0 / bpm;

	float curBeat = curTime / beatDuration;

	float startTimeOfBeat = floor(curBeat); 
	float timeInBeat = fract(curBeat);

	float p = emphasisPower;
	float easedTime = pow(timeInBeat, p) / (pow(timeInBeat, p) + pow(1.0 - timeInBeat, p));

	float easedBeat = startTimeOfBeat + easedTime;
	return easedBeat * beatDuration;
}

// beatPerCycle: アニメーションが1往復するのに必要な拍数 (1.0で1拍、4.0で4拍)
float getFrequencyFactor(float bpm, float beatPerCycle) {
	float BPS = bpm / 60.0;
	float factor = BPS * (2.0 * PI) / beatPerCycle;
	return factor;
}

// ========= SDF Primitives ===========
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));
}

float sdPlane( vec3 p, vec4 n )
{
  // n must be normalized
  return dot(p,n.xyz) + n.w;
}

float sdCausticsReciever(vec3 ro, vec3 ray)
{
	float t = 0.0;
	for (int i = 0; i < ITER; i++) {
		float res = sdPlane(ro + ray * t, vec4(0.0, 1.0, 0.0, 10.0));
		if( res < 0.00001 ) return t;
		t += res;
	}

	return -1.0;
}


// ========= Grid Subdivision & Traversal ===========
//ref: ob5vr https://www.shadertoy.com/view/dsjGDD
struct SubdivResult {
  vec2 size;
  vec2 cell;
  vec2 hash;
};
struct GridResult {
  SubdivResult subdiv;
  float d;
};
GridResult gGridResult;

/**
 * Return the grid cell information the given point belongs to.
 */
SubdivResult subdivision( vec2 p ) {
  SubdivResult result;
  
  result.size = vec2( 0.5 );
  
  for ( int i = 0; i < 5; i ++ ) {
    result.cell = ( floor( p / result.size ) + 0.5 ) * result.size;
    result.hash = hash22( result.cell );
    
    if ( i == 4 || result.hash.x < 0.45 ) {
      break; // no
    } else {
      result.size *= 0.5;
    }
  }
  
  return result;
}


GridResult gridTraversal( vec2 ro, vec2 rd ) {
  GridResult result;

  result.subdiv = subdivision( ro + rd * 1E-3 );

  vec2 src = -( ro - result.subdiv.cell ) / rd;
  vec2 dst = abs( 0.5 * result.subdiv.size / rd );
  vec2 bv = src + dst;
  result.d = min( bv.x, bv.y );
  
  return result;
}

float mapBox(vec3 p, SubdivResult subdiv)
{
	// Repeated
	p.xz -= subdiv.cell;
	float id = subdiv.hash.x;
	float power = 3.0;//(0.5/min(subdiv.size.x, 0.25)); // 3.2
	float easeTime = getEasedTime(gSceneTime, BPM, power);
	float freq = getFrequencyFactor(BPM, 1.) * subdiv.hash.x;
    
	float offset = hash12(subdiv.cell) * 0.06 + sin(easeTime * freq) * 0.06;// - (1./subdiv.size.x*0.2)* 0.01;

	float d = sdBox(p + vec3(0,offset,0), vec3(subdiv.size.x,2.0,subdiv.size.y));
	return d;
}

float map(vec3 p)
{
	return mapBox(p, gGridResult.subdiv);
}

vec2 scene(vec3 ro, vec3 ray)
{
	float t = 0.0;
	for (int i = 0; i < ITER; i++) {
		gGridResult = gridTraversal(ro.xz + ray.xz * t, ray.xz);
		float res = map(ro+ray*t);
		if( res < 0.00001 ) return vec2(t, res);

		t += min(gGridResult.d,res);
	}

	return vec2(-1.0);
}

// 水面との交差 変形前
vec2 sdWaterPlane(vec3 ro, vec3 ray)
{
	float t = 0.0;
	for (int i = 0; i < 64; i++) {
		float res = sdPlane(ro+ray*t, vec4(0.0, 1.0, 0.0, -1.0));
		if( res < 0.0001 ) return vec2(t, res);
		t += res;
	}
	return vec2(-1.0);
}

vec3 waterNormal(vec3 p, float e) {
	vec3 ee = vec3(e, 0., 0.);
	vec3 pp = p * 1.0;
		float h1 = noise(p + ee.xyy);
	float h2 = noise(p - ee.xyy);
	float h3 = noise(p + ee.yyx);
	float h4 = noise(p - ee.yyx);
	vec3 du = vec3(1., 0., h2 - h1);
	vec3 dv = vec3(0., 1., h4 - h3);
	return normalize(cross(du, dv)) * 0.5 + 0.5;
	//return vec3(h1, h2, h3);
}


float caustics(vec3 p, vec3 lp) {
	vec3 ray = normalize(p - lp);

	vec2 shadow = scene(lp, ray);
	float l = distance(lp + ray * shadow.x, p);
	if (l > 0.01) {
		return 0.0;
	}

	vec2 d = sdWaterPlane(lp, ray);
	vec3 waterSurface = lp + ray * d.x;

	vec3 refractRay = refract(ray, vec3(0., 1., 0.), 1.0/1.333);
	float beforeHit = sdCausticsReciever(waterSurface, refractRay);
	vec3 beforePos = waterSurface + refractRay * beforeHit;

	float easeTime = getEasedTime(gSceneTime, BPM * 0.5, 2.);
	vec3 noisePos = waterSurface + vec3(0.,easeTime * 2.0,0.);
	float height = noise(noisePos);
	vec3 deformedWaterSurface = waterSurface + vec3(0., height, 0.);

	refractRay = refract(ray, waterNormal(noisePos, 0.5), 1.0/1.333);
	float afterHit = sdCausticsReciever(deformedWaterSurface, refractRay);
	vec3 afterPos = deformedWaterSurface + refractRay * afterHit;

	float beforeArea = length(dFdx(beforePos)) * length(dFdy(beforePos));
	float afterArea = length(dFdx(afterPos)) * length(dFdy(afterPos));
	return max(beforeArea / afterArea, .001);
}

vec3 normal(vec3 pos, float e )
{
	vec3 eps = vec3(e,0.0,0.0);

	return normalize( vec3(
		   map(pos+eps.xyy) - map(pos-eps.xyy),
		   map(pos+eps.yxy) - map(pos-eps.yxy),
		   map(pos+eps.yyx) - map(pos-eps.yyx) ) );
}

mat3 createCamera(vec3 ro, vec3 ta, float cr )
{
	vec3 cw = normalize(ta - ro);
	vec3 cp = vec3(sin(cr), cos(cr),0.0);
	vec3 cu = normalize( cross(cw,cp) );
	vec3 cv = normalize( cross(cu,cw) );
	return mat3( cu, cv, cw );
}

// ref by xor: https://mini.gmshaders.com/p/tonemaps
vec3 Tonemap_ACES(vec3 x)
{
	const float a = 2.51;
	const float b = 0.03;
	const float c = 2.43;
	const float d = 0.59;
	const float e = 0.14;
	return (x * (a * x + b)) / (x * (c * x + d) + e);
}

vec3 tonemap(vec3 color) {
	color = Tonemap_ACES(color);
	return color;
}

void main()
{
	gSceneTime = mod(time, 120.0);
	// fragment position
	vec2 fp = (gl_FragCoord.xy * 2.0 - resolution.xy) / min(resolution.x, resolution.y);

	// camera
	vec3 ro = vec3((gSceneTime * 0.01) * 10.0, 1.6, (gSceneTime * 0.01) * 10.0) * 2.0;
	vec3 lookat = vec3(ro.x + 1. , ro.y - 5.0, ro.z + 10.);
	mat3 cm = createCamera(ro, lookat, 0.0);
	vec3 rd = cm * normalize(vec3(fp, 4.0));
	
	vec3 focalPoint = ro + rd * 5.0;
	ro.xy += 0.01 * hash22(fp) * max(0., (sin(gSceneTime) * 0.2));
	rd = normalize(focalPoint - ro);
	
	// marching loop
	vec2 res = scene(ro, rd);
	
	vec3 lightPos = ro + vec3(10.0, 3.0, 10.0); //light pos
	vec3 ambient = vec3(0.12,0.11,0.11);
	vec3 fog = vec3(0.2,0.2,0.3);

	float easeTime = getEasedTime(gSceneTime, BPM, 3.2);
	vec3 causticsColorScale = vec3(
	  sin(easeTime * 0.5) * 0.2 + 0.2,
	  cos(easeTime * 0.8) * 0.5 + .5 + 0.3,
	  1.
	 );

	// hit check
	if(res.y > -0.5) {
		vec3 p = ro + rd * res.x;
	   vec3 n = normal(p, 0.0001);
		vec3 v = -normalize(p - lightPos);

		float c = caustics(ro + rd * res.x, lightPos) * 10.;
		vec3 albedo = vec3(0.02,0.15,0.27);

		vec3 co = vec3(c) * causticsColorScale;
		float li = max(dot(v, n), 0.0);
		vec3 col = co * li * albedo + ambient;

		outColor = vec4(mix(col, fog, min(1., log(res.x) / 7.)), 1.0);
	}else{
		outColor = vec4(vec3(0.0), 1.0);
	}
	outColor = vec4(tonemap(outColor.rgb), 1.0);
}