
vec4 grad4(float j, vec4 ip)
{
	const vec4 ones = vec4(1.0, 1.0, 1.0, -1.0);
	vec4 p,s;

	p.xyz = floor( fract (vec3(j) * ip.xyz) * 7.0) * ip.z - 1.0;
	p.w = 1.5 - dot(abs(p.xyz), ones.xyz);
	s = vec4(lessThan(p, vec4(0.0)));
	p.xyz = p.xyz + (s.xyz*2.0 - 1.0) * s.www; 

	return p;
}

#define F4 0.309016994374947451

vec4 simplexNoiseDerivatives (vec4 v)
{
	const vec4  C = vec4( 0.138196601125011,0.276393202250021,0.414589803375032,-0.447213595499958);

	vec4 i  = floor(v + dot(v, vec4(F4)) );
	vec4 x0 = v -   i + dot(i, C.xxxx);

	vec4 i0;
	vec3 isX = step( x0.yzw, x0.xxx );
	vec3 isYZ = step( x0.zww, x0.yyz );
	i0.x = isX.x + isX.y + isX.z;
	i0.yzw = 1.0 - isX;
	i0.y += isYZ.x + isYZ.y;
	i0.zw += 1.0 - isYZ.xy;
	i0.z += isYZ.z;
	i0.w += 1.0 - isYZ.z;

	vec4 i3 = clamp( i0, 0.0, 1.0 );
	vec4 i2 = clamp( i0-1.0, 0.0, 1.0 );
	vec4 i1 = clamp( i0-2.0, 0.0, 1.0 );

	vec4 x1 = x0 - i1 + C.xxxx;
	vec4 x2 = x0 - i2 + C.yyyy;
	vec4 x3 = x0 - i3 + C.zzzz;
	vec4 x4 = x0 + C.wwww;

	i = mod289(i); 
	float j0 = permute( permute( permute( permute(i.w) + i.z) + i.y) + i.x);
	vec4 j1 = permute( permute( permute( permute (   i.w + vec4(i1.w, i2.w, i3.w, 1.0 ))
												   + i.z + vec4(i1.z, i2.z, i3.z, 1.0 ))
												   + i.y + vec4(i1.y, i2.y, i3.y, 1.0 ))
												   + i.x + vec4(i1.x, i2.x, i3.x, 1.0 ));


	vec4 ip = vec4(1.0/294.0, 1.0/49.0, 1.0/7.0, 0.0) ;

	vec4 p0 = grad4(j0,   ip);
	vec4 p1 = grad4(j1.x, ip);
	vec4 p2 = grad4(j1.y, ip);
	vec4 p3 = grad4(j1.z, ip);
	vec4 p4 = grad4(j1.w, ip);

	vec4 norm = taylorInvSqrt(vec4(dot(p0,p0), dot(p1,p1), dot(p2, p2), dot(p3,p3)));
	p0 *= norm.x;
	p1 *= norm.y;
	p2 *= norm.z;
	p3 *= norm.w;
	p4 *= taylorInvSqrt(dot(p4,p4));

	vec3 values0 = vec3(dot(p0, x0), dot(p1, x1), dot(p2, x2)); //value of contributions from each corner at point
	vec2 values1 = vec2(dot(p3, x3), dot(p4, x4));

	vec3 m0 = max(0.5 - vec3(dot(x0,x0), dot(x1,x1), dot(x2,x2)), 0.0); //(0.5 - x^2) where x is the distance
	vec2 m1 = max(0.5 - vec2(dot(x3,x3), dot(x4,x4)), 0.0);

	vec3 temp0 = -6.0 * m0 * m0 * values0;
	vec2 temp1 = -6.0 * m1 * m1 * values1;

	vec3 mmm0 = m0 * m0 * m0;
	vec2 mmm1 = m1 * m1 * m1;

	float dx = temp0[0] * x0.x + temp0[1] * x1.x + temp0[2] * x2.x + temp1[0] * x3.x + temp1[1] * x4.x + mmm0[0] * p0.x + mmm0[1] * p1.x + mmm0[2] * p2.x + mmm1[0] * p3.x + mmm1[1] * p4.x;
	float dy = temp0[0] * x0.y + temp0[1] * x1.y + temp0[2] * x2.y + temp1[0] * x3.y + temp1[1] * x4.y + mmm0[0] * p0.y + mmm0[1] * p1.y + mmm0[2] * p2.y + mmm1[0] * p3.y + mmm1[1] * p4.y;
	float dz = temp0[0] * x0.z + temp0[1] * x1.z + temp0[2] * x2.z + temp1[0] * x3.z + temp1[1] * x4.z + mmm0[0] * p0.z + mmm0[1] * p1.z + mmm0[2] * p2.z + mmm1[0] * p3.z + mmm1[1] * p4.z;
	float dw = temp0[0] * x0.w + temp0[1] * x1.w + temp0[2] * x2.w + temp1[0] * x3.w + temp1[1] * x4.w + mmm0[0] * p0.w + mmm0[1] * p1.w + mmm0[2] * p2.w + mmm1[0] * p3.w + mmm1[1] * p4.w;

	return vec4(dx, dy, dz, dw) * 49.0;
}


// -------------------------------------------------------------------------------------
// Default _persistence should be 0.5
//
vec3 curlNoise( vec3 _p, float _time, int _octaves, float _persistence )
{
	vec4 xNoisePotentialDerivatives = vec4(0.0);
	vec4 yNoisePotentialDerivatives = vec4(0.0);
	vec4 zNoisePotentialDerivatives = vec4(0.0);
	
	float tmpPersistence = _persistence;
	
	for (int i = 0; i < _octaves; ++i)
	{
		float scale = (1.0 / 2.0) * pow(2.0, float(i));
		
		float noiseScale = pow(tmpPersistence, float(i));
		if (tmpPersistence == 0.0 && i == 0) //fix undefined behaviour
		{
			noiseScale = 1.0;
		}
		
		xNoisePotentialDerivatives += simplexNoiseDerivatives(vec4(  _p 								  * pow(2.0, float(i)), _time)) * noiseScale * scale;
		yNoisePotentialDerivatives += simplexNoiseDerivatives(vec4(( _p + vec3(123.4, 129845.6, -1239.1)) * pow(2.0, float(i)), _time)) * noiseScale * scale;
		zNoisePotentialDerivatives += simplexNoiseDerivatives(vec4(( _p + vec3(-9519.0, 9051.0, -123.0))  * pow(2.0, float(i)), _time)) * noiseScale * scale;
	}
	
	//compute curl
	vec3 curlVal = vec3( zNoisePotentialDerivatives[1] - yNoisePotentialDerivatives[2],
						 xNoisePotentialDerivatives[2] - zNoisePotentialDerivatives[0],
						 yNoisePotentialDerivatives[0] - xNoisePotentialDerivatives[1] );
	
	
	return curlVal;
}
