#version 430
uniform float ttime;
uniform float DeltaT;
uniform float PosY;

//uniform image2D destTexPos;
//uniform image2D destTexCol;
float   hash(float v);
vec2    rhash(vec2 uv);

layout(binding = 0) writeonly uniform image2D destTexPos;
layout(binding = 1) writeonly uniform image2D destTexVel;
layout(binding = 2) writeonly uniform image2D destTexCol;

uniform float type;
uniform sampler2D PositionTexture;
uniform sampler2D ColorTexture;
uniform sampler3D CurlTexture;

uniform float effect1;
uniform float effect2;
uniform float effect3;
uniform float effect4;
uniform float effect5;
uniform float effect6;


uniform sampler2D OriginalPositionTexture;
uniform sampler2D OriginalColorTexture;
uniform sampler2D OriginalNormalTexture;
uniform sampler2D OriginalPositionTexture_Text;
uniform sampler2D OriginalColorTexture_Text;

uniform sampler2D Swirl_Texture;

uniform float k0;
uniform float k1;
uniform float k2;
uniform float k3;
uniform float k4;
uniform float k5;

uniform sampler1D Original1DTexture;
layout(local_size_x = 16, local_size_y = 16) in;

layout(std140, binding = 4) buffer Pos { vec4 Positions[]; };
layout(std140, binding = 5) buffer Vel { vec4 Velocities[]; };
layout(std140, binding = 6) buffer Col { vec4 Colors[]; };


const float InverseMaxInt = 1.0 / 4294967295.0;

float hash(float n)
{
	return fract(sin(n) * 43758.5453);
}

vec2 rhash(vec2 uv) {
	const mat2 t = mat2(.12121212, .13131313, -.13131313, .12121212);
	const vec2 s = vec2(1e4, 1e6);
	uv *= t;
	uv *= s;
	return  fract(fract(uv / s) * uv);
}



float voronoi(const in vec2 uv)
{
	vec2 p = floor(uv);
	vec2 f = fract(uv);
	float v = 0.;
	for (int j = -1; j <= 1; j++)
		for (int i = -1; i <= 1; i++)
		{
			vec2 b = vec2(i, j);
			vec2 r = b - f + rhash(p + b);
			v += 1. / pow(dot(r, r), 8.);
		}
	return pow(1. / v, 0.0625);
}




vec3 Rotate(vec3 Input, mat4 RotationM) {
	vec4 Inn = Input.xyzz;
	Inn.w = 1.0;
	Inn = RotationM * Inn;

	float fc = 1.0 - clamp(pow(length(Input) * 1.1, 1.0) + 0.0, 0.3, 1.0);
	fc = 0.4 + 0.4 * sin(1 * length(Input) + ttime);
	fc *= clamp((ttime - 100) * 0.5, 0.0, 1.0);
	fc = 0;
	Inn.xyz += (Input - Inn.xyz) * fc;

	return Inn.xyz;

}

vec2 rotate(vec2 v, float a) {
	float s = sin(a);
	float c = cos(a);
	mat2 m = mat2(c, s, -s, c);
	return m * v;
}



float randhash(uint seed, float b)
{
	uint i = (seed ^ 12345391u) * 2654435769u;
	i ^= (i << 6u) ^ (i >> 26u);
	i *= 2654435769u;
	i += (i << 5u) ^ (i >> 12u);
	return float(b * i) * InverseMaxInt;
}
vec3 randhash3(uint seed) {
	float x = randhash(seed, 1.0) - 0.5;
	float y = randhash(seed + 413, 1.0) - 0.5;
	float z = randhash(seed * 2 + 12, 1.0) - 0.5;
	return vec3(x, y, z);

}

vec3 GetCurl(vec3 Pos) {

	return texture(CurlTexture, Pos.xyz).xyz;
}



float Octave(vec2 UV) {
	float ff = 0.5;
	float fv = int(UV.x * ff * 20) / (ff * 20.0);
	float fvy = int(UV.y * ff * 40) / (ff * 40.0);
	fv = 0;
	fvy = 0;
	float fy = voronoi(13 * vec2(-fv + UV.x + ttime * 0.073, fvy + UV.y + ttime * 0.03 + ttime * 0.0));
	fy *= voronoi(11 * vec2(UV.x - ttime * 0.053, fvy + UV.y - ttime * 0.02 + ttime * 0.0));
	// fy+=sin(voronoi (6*vec2(UV.x+ttime*0.073,UV.y-ttime*0.02+ttime*0.0)));
	fy += (voronoi(10.4 * vec2(UV.x - ttime * 0.043, fv + UV.y - ttime * 0.023 + ttime * 0.0)));

	fy *= 1 * voronoi(4 * vec2(UV.x - ttime * 0.043, fvy + UV.y + ttime * 0.0 + ttime * 0.0));
	fy += 1.4 * voronoi(33.4 * vec2(fv + UV.x - ttime * 0.042, UV.y + fvy - ttime * 0.01 + ttime * 0.0)) * voronoi(6.4 * vec2(fv + UV.x - ttime * 0.033, UV.y + fvy - ttime * 0.03 + ttime * 0.0))
		;

	fy += 0.3 * sin(ttime + 15.0 * length(vec2(1.0, 0.2) * (UV - vec2(0.5))));
	//fy=0.5+0.5*sin(fy*2);

	// fy*=0.5;
	//fy=voronoi (vec2(UV.x,UV.y)*8.0)*0.0;
	return fy;



}



float Octave2(vec2 UV) {
	float ff = 0.5;
	float fv = int(UV.x * ff * 20) / (ff * 20.0);
	float fvy = int(UV.y * ff * 40) / (ff * 40.0);
	fv = 0;
	fvy = 0;
	UV.y *= 0.6;
	float fy = voronoi(13 * vec2(-fv + UV.x + ttime * 1.043, fvy + UV.y + ttime * 0.03 + ttime * 0.0));
	fy *= voronoi(11 * vec2(UV.x + ttime * 0.043, fvy + UV.y - ttime * 0.02 + ttime * 0.0));
	fy += voronoi(6 * vec2(UV.x + ttime * 0.143, UV.y - ttime * 0.02 + ttime * 0.0));
	fy += voronoi(10.4 * vec2(UV.x + ttime * 0.043, fv + UV.y - ttime * 0.02 + ttime * 0.0));
	fy += 1.4 * voronoi(6.4 * vec2(fv + UV.x + ttime * 0.053, UV.y + fvy - ttime * 0.01 + ttime * 0.0))
		* voronoi(6.4 * vec2(fv + UV.x + ttime * 0.033, UV.y + fvy - ttime * 0.03 + ttime * 0.0))
		;

	fy *= voronoi(4 * vec2(UV.x - ttime * 0.223, fvy + UV.y - ttime * 0.01 + ttime * 0.0));
	fy += 0.3 * sin(ttime + 15.0 * length(vec2(1.0, 0.0) * (UV - vec2(0.5))));
	//fy=0.5+0.5*sin(fy*2);
	return (fy);



}


float noise(in vec3 x)
{
	vec3 p = floor(x);
	vec3 f = fract(x);

	f = f * f * (3.0 - 2.0 * f);

	float n = p.x + p.y * 57.0 + 113.0 * p.z;

	float res = mix(mix(mix(hash(n + 0.0), hash(n + 1.0), f.x),
		mix(hash(n + 57.0), hash(n + 58.0), f.x), f.y),
		mix(mix(hash(n + 113.0), hash(n + 114.0), f.x),
			mix(hash(n + 170.0), hash(n + 171.0), f.x), f.y), f.z);
	return res;
}
float fbm(vec3 p)
{
	float f;
	f = 0.5000 * noise(p); p = p * 2.02;
	f += 0.2500 * noise(p); p = p * 2.03;
	f += 0.1250 * noise(p);
	return f;
}



float mandelbrot(in vec2 c)
{
#if 1
	{
		float c2 = dot(c, c);
		// skip computation inside M1 - https://iquilezles.org/articles/mset1bulb
		if (256.0 * c2 * c2 - 96.0 * c2 + 32.0 * c.x - 3.0 < 0.0) return 0.0;
		// skip computation inside M2 - https://iquilezles.org/articles/mset2bulb
		if (16.0 * (c2 + 2.0 * c.x + 1.0) - 1.0 < 0.0) return 0.0;
	}
#endif

	const float B = 256.0;
	float n = 0.0;
	vec2 z = vec2(0.0);
	for (int i = 0; i < 512; i++)
	{
		z = vec2(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c;
		if (dot(z, z) > (B * B)) break;
		n += 1.0;
	}

	if (n > 511.0) return 0.0;

	// ------------------------------------------------------
	// smooth interation count
	//float sl = l - log(log(length(z))/log(B))/log(2.0);

	// equivalent optimized smooth interation count
	float sn = n - log2(log2(dot(z, z))) + 4.0;

	float al = smoothstep(-0.1, 0.0, sin(0.5 * 6.2831 * ttime));
	return mix(n, sn, al);
}


vec3 cartesianToSpherical(vec3 pos) {
	float r = length(pos);                     // Radial distance
	float theta = atan(pos.z, pos.x);          // Azimuthal angle (angle in XZ-plane)
	float phi = acos(pos.y / r);               // Polar angle (angle from Y-axis)

	return vec3(r, theta, phi);                // (r, theta, phi)
}


#define INVOCATION_SIZE (3184)
void main() {

	float TotalParticles = 1.0 / float(3184 * 3184);

	uint MyIndex = gl_GlobalInvocationID.x + gl_GlobalInvocationID.y * INVOCATION_SIZE;
	float fMyIndex = float(MyIndex);
	ivec2 storePos = ivec2(gl_GlobalInvocationID.xy);
	float fMyIndexX = float(gl_GlobalInvocationID.x) * 1.0 / 3184.0;
	float fMyIndexY = float(gl_GlobalInvocationID.y) * 1.0 / 3184.0;

	vec2 fMyIndexXY = vec2(fMyIndexX, fMyIndexY);



	vec3 OriginalPos = 0.5*randhash3(MyIndex * 2346);

	//OriginalPos *= 0;
	
	
	if (ttime > 19.5)
	//if (int (ttime)%2==0)
	{

		OriginalPos = 0.3 * normalize(OriginalPos);

		vec3 cart = cartesianToSpherical(OriginalPos);
		float fm = 
			pow(1.0 * clamp (sin(15.0 * cart.y),0.0,1.0) *
				clamp (sin(15.0 * cart.z),0.0,1.0), 7.0);
		OriginalPos *= 1+2.0*fm;


	}
	OriginalPos.xz = rotate(OriginalPos.xz, 0.3 * ttime); //

	float kk3 = k3;
	float flm = fMyIndex * TotalParticles;


	if (ttime < 66.0) kk3 = 0;
	kk3 = clamp((ttime - 61 + 5.0 * flm) * 0.3, 0.0, 1.0);
	kk3 -= clamp((ttime - 74 + 1.0 * flm) * 0.3, 0.0, 1.0);

	kk3 += clamp((ttime - 12 - 1.0 * flm) * 0.7, 0.0, 1.0);
	kk3 -= clamp((ttime - 16 - 1.0 * flm) * 0.3, 0.0, 1.0);

	kk3 += clamp((ttime - 118 - 3.0 * flm) * 0.017, 0.0, 1.0);
	kk3 -= clamp((ttime - 153 - 1.0 * flm) * 0.4, 0.0, 1.0);

	//kk3 = 1;
	//if (ttime > 14.8)
	if (kk3 > 0) {
		
		float fm = fMyIndex * TotalParticles;
		fm *= 1000;
//		fm += ttime * 1;
		OriginalPos = vec3(sin(fm), sin(fm * 1.3), sin(fm * 1.524))*0.3;
		OriginalPos = normalize(OriginalPos)*0.3;

		OriginalPos.xy = rotate(OriginalPos.xy, 0.1*ttime); //


	}




	//if (ttime<80)
	//if (int (MyIndex%5)==0)
	if (ttime<104+flm*2.0)
	if (ttime > 65) {

		OriginalPos= 0.5 * randhash3(MyIndex * 2346);
		OriginalPos = normalize(OriginalPos)*randhash3 (MyIndex*4451).x;
	};

	/*
	if (ttime>45) {

		float fm = fMyIndex * TotalParticles;
		fm *= 100;
		//		fm += ttime * 1;
		OriginalPos = vec3(sin(fm), sin(fm * 1.3), sin(fm * 1.524)) * 0.3;
		OriginalPos = normalize(OriginalPos) * 0.2;
		Positions[MyIndex].xyz += (OriginalPos.xyz - Positions[MyIndex].xyz) * 0.01;
		//OriginalPos.xy = rotate(OriginalPos.xy, 0.1 * ttime); //


	}
	*/

	if (ttime>150)
	//if (1.0 > 1.0 - fMyIndex * TotalParticles)
	{

		Positions[MyIndex].xz = rotate(Positions[MyIndex].xz, length(Velocities[MyIndex].xyz) * 1.0 * -10.1); //
		Positions[MyIndex].xy = rotate(Positions[MyIndex].yx, 1.0 * 1.0 * 0.1); //


	};



	if (ttime < 2.0+0.0*clamp (OriginalPos.y,0.0,100.0))
	//if (sin (ttime)<sin (OriginalPos.z))
	{
		Positions[MyIndex].xyz = OriginalPos.xyz + vec3(0, 0, 0.0) * ttime;
		//Positions[MyIndex].xyz = normalize(Positions[MyIndex].xyz);
	//	Positions[MyIndex].y = 0;


	}
	else
	{

		/*
		if (k2 > 0.0)
		{
		vec3 cart = cartesianToSpherical(Positions[MyIndex].xyz);
		float fm = 2.0*sin(10.0*cart.y)* sin(10.0 * cart.z);
		fm = clamp(fm, 0.98, 1.02);

		//Positions[MyIndex].xyz += (Positions[MyIndex].xyz*1.5 - Positions[MyIndex].xyz) * fm;
		Positions[MyIndex].xyz = 0.8*normalize(Positions[MyIndex].xyz);
		}
		*/

		Positions[MyIndex].xyz += Velocities[MyIndex].xyz; // OriginalPos.xyz + vec3(0, 0, -0.01) * ttime;
		float theta = 0.1*pow (sin(OriginalPos.y*20.0+3.0*ttime),20.0);

		float twirl = 0.0;

		twirl = -0.02 * clamp((ttime - 86) * 0.1, 0.0, 1.0);
		twirl += 0.02 * clamp((ttime - 146) * 0.1, 0.0, 1.0);

				Positions[MyIndex].xz = rotate(Positions[MyIndex].xz, twirl +1.0*clamp(0.08 - 0.3 * length(Positions[MyIndex].xyz), 0.0, 1.0));
		
				//Positions[MyIndex].xy = rotate(Positions[MyIndex].xy, Positions[MyIndex].y*k2);

				float idf = 10.0;
				if (k2 > 0.0) idf = 10000.0; // +clamp(k2, 0.0, 1.0) * clamp(k2, 0.0, 1.0) * 2000;
				//idf = 10000.0;
		//Positions[MyIndex].xz = rotate(Positions[MyIndex].xz, -clamp(2.3 * length(Velocities[MyIndex].xyz), 0.0, 1.0));
				/*				Positions[MyIndex].xyz = vec3 (ivec3(idf*Positions[MyIndex].xyz))/idf;
				Positions[MyIndex].x += (-0.5 + 1.0 * hash(58.76*float(MyIndex))) * 0.01;
				Positions[MyIndex].y += (-0.5 + 1.0 * hash(13.2*float(MyIndex))) * 0.01;
				Positions[MyIndex].z += (-0.5 + 1.0 * hash(434.3*float(MyIndex))) * 0.01;
				*/
				//Positions[MyIndex].xyz *= 1.0 - clamp (k1,0.0,1.0)*GetCurl(Positions[MyIndex].xyz);

		//Positions[MyIndex].xyz += (OriginalPos - Positions[MyIndex].xyz) * theta; // 0.2 * voronoi(10.0 * Positions[MyIndex].xz); // clamp(sin(length(OriginalPos.xz) * 12.0 + 1.0 * ttime), 0.0, 1.0)* k0;
	//	Positions[MyIndex].y = clamp(Positions[MyIndex].y, 0.0, 1000.0);
	}

	//Positions[MyIndex].x = -2 + mod(Positions[MyIndex].x + 2.0, 4.0);
	//Positions[MyIndex].z = -2 + mod(Positions[MyIndex].z + 2.0, 4.0);


	Positions[MyIndex].xyz += (
		
		3.0*GetCurl(Positions[MyIndex].xyz*3.1)*vec3 (1,1,1) - Positions[MyIndex].xyz) *
		clamp (k1,0.0,1.0)*0.05;



/*
// ending
	Positions[MyIndex].xyz += (

		3.0 * GetCurl(Positions[MyIndex].xyz * 1.1) * vec3(1, 1, 1) - Positions[MyIndex].xyz) *
	clamp(k1, 0.0, 1.0) * 0.05;
	*/




	/*
	Positions[MyIndex].xyz += (

		vec3 (0,0,0) - Positions[MyIndex].xyz) *
		clamp(2.0*k3-length (Positions[MyIndex].xyz)*0.1, -1.0, 1.0) * 0.02;
		*/

	//Positions[MyIndex].xyz *= 0.999;
	Colors[MyIndex].xyz = OriginalPos.xyz; // , 1, 1);
	
/*	Velocities[MyIndex].xyz = 0.01 / (ttime * ttime) *
		GetCurl(vec3(-ttime * 0.0) + vec3(15.0) + 0.25 * Positions[MyIndex].xyz) * 1.0 + vec3(0.0, 0.0, 0.0);
		*/

	float kk0 = k0;
	kk0 += 10.0 * clamp((ttime - 65) * 0.2, 0.0, 1.0);
	kk0 -= 10.0 * clamp((ttime - 105) * 0.2, 0.0, 1.0);


	Velocities[MyIndex].xyz = 
	vec3(0.00, -0.0, 0.0) + vec3(1, 1, 1) *kk0* 0.01 *
		GetCurl(vec3(-ttime * 0.00) + 1.0 * 1 * Positions[MyIndex].xyz) * 1.0 + vec3(0.0, 0.0, 0.0);
	
	Velocities[MyIndex].xyz +=
		vec3(0.00, -0.0, 0.0) + vec3(1, 1, 1) * kk0 * 0.01 *
		GetCurl(vec3(-ttime * 0.00) + 1.0 * 0.452 * Positions[MyIndex].xyz) * 1.0 + vec3(0.0, 0.0, 0.0);
	
/*	Velocities[MyIndex].xyz +=
		vec3(0.00, -0.0, 0.0) + vec3(1, 1, 1) * k2 * 0.04 *
		GetCurl(vec3(-ttime * 0.00) + 1.0 * 300.352 * OriginalPos.xyz) * 1.0 + vec3(0.0, 0.0, 0.0);
		*/


	Velocities[MyIndex].xyz += vec3(0, 0.0015, 0);
	if (Positions[MyIndex].y > 1)
		Positions[MyIndex].y -= 2;
	

	if (kk3 > 1.0-fMyIndex * TotalParticles) {
		Positions[MyIndex].xyz +=

			(OriginalPos.xyz - Positions[MyIndex].xyz) * 0.04;

	};



//	Velocities[MyIndex] += (Velocities[1+MyIndex] - Velocities[0+MyIndex])*clamp (k2,0.0,1.0);

	/*
	Velocities[MyIndex].xyz +=
		vec3(0.00, -0.0, 0.0) + vec3(1, 1, 1) * k0 * 0.01 *
		GetCurl(vec3(-ttime * 0.02) + 1.0 * 0.032 * Positions[MyIndex].xyz) * 1.0 + vec3(0.0, 0.0, 0.0);
		*/
	if (MyIndex % 100 == 0) {
//		Velocities[MyIndex].x += 0.005;
	}
	

	//Velocities[MyIndex].xyz = vec3 (1.0)*0.01* voronoi(10.0*Positions[MyIndex].xz);

	

	
	
	/*
	if (effect1 > 0.0)
	//if (int (ttime)%2==0)
	{

		if (int(MyIndex) % 2 == 0) 
		{

			Positions[MyIndex].xyz = -Positions[MyIndex+1].xyz;
			 //= normalize(Positions[MyIndex].xyz) * 1.1;

		}

	};
	*/

//	Velocities[MyIndex].xyz += 0.001 *
	//	GetCurl(vec3(-ttime * 0.0) + vec3(15.0) + 0.55 * Positions[MyIndex].xyz) * 1.0 + vec3(0.0, 0.0, 0.0);
	//Velocities[MyIndex].x = 0;

	//Velocities[MyIndex].xyz += (target- Velocities[MyIndex].xyz ) * 0.1;

	//Colors[MyIndex].xyz = vec3 (1,0,1);
	//Positions[MyIndex].xyz = OriginalPos.xyz;
	imageStore(destTexPos, storePos, Positions[MyIndex]);
	imageStore(destTexVel, storePos, Velocities[MyIndex]);
	imageStore(destTexCol, storePos, Colors[MyIndex]);
}
