#if !defined CLOUD_FUNCTIONS
	#define CLOUD_FUNCTIONS
	float get3DNoise(in vec3 position) {
		float flr = floor(position.z);
		vec2 coord = (position.xy * 0.015625) + (flr * 0.265625);
		vec2 noise = texture(noisetex, coord).xy;
		return mix(noise.x, noise.y, position.z - flr);
	}

	float fbm3D(vec3 position, in int octaves) {
		float amplitude = 1.0;
		float V = amplitude * get3DNoise(position); 

		float persistance = 0.63;
		float lacunarity = 3.15;
		
		float invWeight = 0.0;
		for (int i = 0; i < octaves; ++i) {
			V += amplitude * get3DNoise(position); 
			amplitude *= persistance; 
			position *= lacunarity;
			invWeight += amplitude;
		}
		
		return V / invWeight; 
	}

	float fbm2D(vec2 position, in int octaves) {
		float amplitude = 1.0;
		float V = amplitude * texture(noisetex, position).r; 

		float persistance = 0.63;
		float lacunarity = 3.25;
		
		float invWeight = 0.0;
		for (int i = 0; i < octaves; ++i) {
			V += amplitude * texture(noisetex, position).r; 
			amplitude *= persistance; 
			position *= lacunarity;
			invWeight += amplitude;
		}
		
		return V / invWeight; 
	}

	vec3 distortFBM3D(vec3 position) {
		float amplitude = 1.0;
		vec3 V = amplitude * vec3(get3DNoise(position), get3DNoise(position + vec3(3.33, 5.71, 1.96)), get3DNoise(position + vec3(7.77, 2.65, 4.37))); 

		float persistance = 0.35;
		float lacunarity = 3.25;
		
		float invWeight = 0.0;
		for (int i = 0; i < 2; ++i) {
			V += amplitude * vec3(get3DNoise(position), get3DNoise(position + vec3(3.33, 5.71, 1.96)), get3DNoise(position + vec3(7.77, 2.65, 4.37))); 
			amplitude *= persistance; 
			position *= lacunarity;
			invWeight += amplitude;
		}
		
		return V / invWeight; 
	}

	float remap(float value, float oldLow, float oldHight, float newLow, float newHight) {
		return newLow + (value - oldLow) * (newHight - newLow) / (oldHight - oldLow);
	}

	float heightAlter(float percentHeight, float localizedCoverage) {
		float returnValue =  saturate(remap(percentHeight, 0.0, 0.07, 0.0, 1.0));
		float stopHeight = saturate(localizedCoverage + 0.6); 
		returnValue *= saturate(remap(percentHeight, stopHeight * 0.2, stopHeight, 1.0, 0.0));

		returnValue = pow(returnValue, saturate(remap(percentHeight, 0.65, 0.95, 1.0, (1 - cloudAnvilAmount * mix(globalCoverage, 0.4, rainStrength)))));

		return returnValue;
	}

	float densityAlter(float percentHeight, float localizedCoverage) {
		float returnValue = percentHeight;
		returnValue *= saturate(remap(percentHeight, 0.0, 0.2, 0.0, 1.0));
		returnValue *= localizedCoverage * 2.0;

		return returnValue;
	}

	float calculateCloudShape(in vec3 position) {
		if(position.y > cloudsMaxAltitude || position.y < cloudsAltitude) return 0.0;
		vec3 cloudPosition = position / (cloudsThickness * 64.0);

		float localizedCoverage = saturate(fbm2D((cloudPosition.xz * rcp(10.0)) + (frameTimeCounter * 0.000004), 3).r * (3.0-2.5));
		float coverage = mix(globalCoverage, 0.4, rainStrength) + mix(localizedCoverage, 0.0, rainStrength);

		float height = position.y - cloudsAltitude;
		float normalizeHeight = height * (1.0 / cloudsThickness);
		float heightAlteration = heightAlter(normalizeHeight, localizedCoverage);
		float densityAlteration = densityAlter(normalizeHeight, localizedCoverage) * 2.5;
			  densityAlteration = pow(densityAlteration, 1.8);

		float mainNoise  = fbm2D(cloudPosition.xz + (frameTimeCounter * 0.0004), 3) * 0.4;
			mainNoise -= coverage;
			mainNoise += heightAlteration;
		if(mainNoise <= 0.0) return 0.0;

		cloudPosition = cloudPosition * 100.0;
		float detailNoise = fbm3D(cloudPosition, 4) * 0.4;

		float cloudNoise = mainNoise - detailNoise;

		float density = saturate(16.0 * cloudNoise * densityAlteration);

		return density;
	}

	float henyeyGreensteinPhase(float cosTheta, float g) {
		cFloat norm = 0.25/pi;

		float gg = g * g;
		return norm * ((1.0-gg) / pow(1.0+gg-2.0*g*cosTheta, 3.0/2.0));
	}

	float cloudsPhase(in float cosTheta, in float cn) {
		cFloat backScatterStrength = 0.2;
		cFloat peakStrength = 0.15;

		float forward = mix(henyeyGreensteinPhase(cosTheta, 0.80 * cn), henyeyGreensteinPhase(cosTheta, 0.90 * cn), peakStrength);
		float backward = henyeyGreensteinPhase(cosTheta, -0.50 * cn);

		return mix(forward, backward, backScatterStrength);
	}
#endif