# version 440
//layout(binding = 0) buffer dcE { float E [  ]; };
//layout(binding = 1) buffer dcEnew { float Enew [  ]; };
//layout(binding = 2) buffer dcU { float U [  ]; };
//layout(binding = 3) buffer dcF { int F [  ]; };
//layout(rgba8, binding = 4) uniform image2D img;

layout(binding = 16) buffer dcE { float E [  ]; };
layout(binding = 17) buffer dcEnew { float Enew [  ]; };
layout(binding = 18) buffer dcU { float U [  ]; };
layout(binding = 19) buffer dcF { int F [  ]; };
//layout(binding = 20) buffer waveparticle { Particle p [ ]; };

layout(local_size_x = 10, local_size_y = 10, local_size_z = 10) in;

int per(int x, int nx)
{
	if (x < 0) x += nx;
	if (x >= nx) x -= nx;
	return x;
}

uniform int NX;
uniform int NY;
uniform int NZ;
uniform float mytime;
const float DT = 0.025;    // time dt
const float DX = 0.1;       // spatial dx
const float FREQ = 10.5;    // pulsating source frequency
const float EPS = 0.333;    // laplace operator
#define C_SRC 1
#define C_BND 2

// palette
vec4 color(float t, int k)
{
	const int siz = 12;
	float coltab[] = {
	0.5, 0.5, 0.5, 0.5, 0.65, 0.5, 1.0, 1.0, 1.0,  0.00, 0.33, 0.67,
	0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, 0.00, 0.10, 0.20 ,//!!
	0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0, 0.30, 0.20, 0.20 ,//good
	0.5, 0.5, 0.5, 0.5, 0.5, 0.5,  1.0, 1.0, 0.5,  0.80, 0.90, 0.30,
	0.5, 0.5, 0.5, 0.5, 0.5, 0.5,  1.0, 0.7, 0.4,  0.00, 0.15, 0.20,//gld
	0.5, 0.5, 0.5, 0.5, 0.5, 0.5,  2.0, 1.0, 0.0,  0.50, 0.20, 0.25, // nice bumpmap
	0.8, 0.5, 0.4, 0.2, 0.4, 0.2,  2.0, 1.0, 1.0,  0.00, 0.25, 0.25};
	int of = siz * k;   // offset
	vec4 col;
	col.r = coltab[0 + of] + coltab[3 + of] * cos(2 * 3.1416 * (coltab[6 + of] * t + coltab[9 + of]));
	col.g = coltab[1 + of] + coltab[4 + of] * cos(2 * 3.1416 * (coltab[7 + of] * t + coltab[10 + of]));
	col.b = coltab[2 + of] + coltab[5 + of] * cos(2 * 3.1416 * (coltab[8 + of] * t + coltab[11 + of]));
	col.a = 1;
	return col;
}

void main()
{
	int i, j, k, idx;
	i = int(gl_GlobalInvocationID.x);
	j = int(gl_GlobalInvocationID.y);
	k = int(gl_GlobalInvocationID.z);
	idx = i + NX * j + NX * NY * k;
	int ip, im, jp, jm, kp, km;
	if (i < NX - 1) ip = i + 1; else ip = 0; // warunki cykliczne
	if (i > 0) im = i - 1; else im = NX - 1;
	if (j < NY - 1) jp = j + 1; else jp = 0;
	if (j > 0) jm = j - 1; else jm = NY - 1;
	if (k < NZ - 1) kp = k + 1; else kp = 0;
	if (k > 0) km = k - 1; else km = NZ - 1;

	/*U[idx] = U[idx] +
			(
			(1 - EPS) * (E[ip+NX*j] + E[im+NX*j] + E[i+NX*jm] + E[i+NX*jp] - 4 * E[idx])
			+ EPS * (E[ip+NX*jp] + E[ip+NX*jm] + E[im+NX*jm] + E[im+NX*jp] - 4 * E[idx])
			) * (DT) / (DX * DX);*/
	//https://scicomp.stackexchange.com/questions/35920/3d-laplacian-operator

	U[idx] = U[idx] +
			(

			(E[ip + NX * j + NX * NY * k] + E[im + NX * j + NX * NY * k] + E[i + NX * jm + NX * NY * k] + E[i + NX * jp + NX * NY * k]
			+ E[i + NX * j + NX * NY * kp] + E[i + NX * j + NX * NY * km] - 6 * E[idx])
			) * (DT) / (DX * DX);

	// solve equation for E
	// dE/dt = U
	Enew[idx] = E[idx] + U[idx] * (DT);

	// start boundaries: pulsating sources

	float AMPLITUDE = 0.2;

	//if(mytime<1)
	if (F[idx] == C_SRC)
		Enew[idx] = sin(mytime * FREQ) * AMPLITUDE;

	// boundaries
	if (F[idx] == C_BND)
	{
		Enew[idx] = 0.0;
		U[idx] = 0.0f;
	}

	//Enew[idx] = 1.0 * (sqrt((i - NX / 2) * (i - NX / 2) + (j - NY / 2) * (j - NY / 2) + (k - NZ / 2) * (k - NZ / 2)));

	/*float b = Enew[idx] * 45.0;
	vec4 col1 = color(b, 1);
	vec4 col2 = color(b, 3);
	vec4 col = mix(col1, col2, abs(sin(mytime*0.2)));
	imageStore(img, ivec2(gl_GlobalInvocationID.xy), col);*/
}

