// Procedural world-space cloud layer for NOISE.
// Requires noise_common.glsl and noise_options.glsl to be included first.

#ifndef NOISE_CLOUDS_GLSL
#define NOISE_CLOUDS_GLSL

// The grayscale structural height field controls the local cloud ceiling:
// black is Y=300 and white is Y=500. These wider constants are only conservative
// ray bounds; actual density uses the local height map and a softly varying base.
const float noiseCloudHeightDark = 300.0;
const float noiseCloudHeightLight = 500.0;
const float noiseCloudBottom = 268.0;
const float noiseCloudTop = 504.0;

float noiseCloudHash(vec2 cell) {
	return fract(sin(dot(cell, vec2(127.1, 311.7))) * 43758.5453123);
}

float noiseCloudValue(vec2 position) {
	vec2 cell = floor(position);
	vec2 local = fract(position);
	vec2 curve = local * local * (3.0 - 2.0 * local);
	float a = noiseCloudHash(cell);
	float b = noiseCloudHash(cell + vec2(1.0, 0.0));
	float c = noiseCloudHash(cell + vec2(0.0, 1.0));
	float d = noiseCloudHash(cell + vec2(1.0, 1.0));
	return mix(mix(a, b, curve.x), mix(c, d, curve.x), curve.y);
}

float noiseCloudFbm(vec2 position) {
	float value = noiseCloudValue(position) * 0.57143;
	position = mat2(0.80, 0.60, -0.60, 0.80) * position * 2.03 + vec2(7.1, 3.7);
	value += noiseCloudValue(position) * 0.28571;
	position = mat2(0.80, 0.60, -0.60, 0.80) * position * 2.07 + vec2(5.4, 9.2);
	value += noiseCloudValue(position) * 0.14286;
	return value;
}

// Read the authored seamless grayscale PNG. The weather correction is zero at
// exact black/white, preserving the promised Y=300/Y=500 endpoints, and only
// nudges midtones so covered regions tend to rise as one coherent cloud system.
vec2 noiseCloudHeightBounds(vec2 weatherPosition, float weather) {
	vec2 heightUv = fract(weatherPosition * 0.31 + vec2(0.173, 0.417));
	float heightMap = textureLod(noiseCloudHeightMap, heightUv, 0.0).r;
	float endpointLock = 4.0 * heightMap * (1.0 - heightMap);
	heightMap = noiseSaturate(heightMap +
		(weather - 0.5) * 0.16 * endpointLock);
	float ceilingHeight = mix(noiseCloudHeightDark,
		noiseCloudHeightLight, heightMap);
	float undersideMap = noiseCloudValue(weatherPosition * 1.43 +
		vec2(-5.3, 12.8));
	float baseHeight = mix(270.0, 284.0, undersideMap);
	return vec2(baseHeight, ceilingHeight);
}

// A generated 64^3 RGBA Perlin-Worley volume is packed into an 8x8 atlas.
// One-pixel periodic gutters around every slice make these two bilinear reads
// equivalent to seamless trilinear 3D texture sampling.
vec4 noiseCloudVolumeSlice(vec2 position, float sliceIndex) {
	const float volumeResolution = 64.0;
	const float atlasTiles = 8.0;
	const float atlasTileSize = 66.0;
	const float atlasResolution = atlasTiles * atlasTileSize;
	float wrappedSlice = mod(sliceIndex, volumeResolution);
	vec2 tile = vec2(mod(wrappedSlice, atlasTiles),
		floor(wrappedSlice / atlasTiles));
	vec2 atlasPixel = tile * atlasTileSize + vec2(1.5) +
		fract(position) * volumeResolution;
	return textureLod(noisetex, atlasPixel / atlasResolution, 0.0);
}

vec4 noiseCloudVolume(vec3 position) {
	float slice = fract(position.z) * 64.0;
	float lowerSlice = floor(slice);
	float sliceFraction = fract(slice);
	vec4 lowerValue = noiseCloudVolumeSlice(position.xy, lowerSlice);
	vec4 upperValue = noiseCloudVolumeSlice(position.xy, lowerSlice + 1.0);
	return mix(lowerValue, upperValue,
		sliceFraction * sliceFraction * (3.0 - 2.0 * sliceFraction));
}

float noiseCloudHenyeyGreenstein(float cosineAngle, float anisotropy) {
	float anisotropySquared = anisotropy * anisotropy;
	return (1.0 - anisotropySquared) /
		max(pow(1.0 + anisotropySquared -
			2.0 * anisotropy * cosineAngle, 1.5), 1e-3);
}

float noiseCloudDensityLodHeight(vec3 worldPosition, vec2 wind,
	float rainStrength, float detailAmount, out float height01) {
	height01 = 0.0;
	if (worldPosition.y <= noiseCloudBottom || worldPosition.y >= noiseCloudTop)
		return 0.0;

	// Layout / weather is intentionally two-dimensional. It only says where a
	// cloud may exist and which height profile it should use; it never becomes
	// volume density by itself.
	vec2 weatherPosition = worldPosition.xz * 0.00185 + wind * 0.42;
	float weather = noiseCloudFbm(weatherPosition);
	vec2 heightBounds = noiseCloudHeightBounds(weatherPosition, weather);
	height01 = (worldPosition.y - heightBounds.x) /
		max(heightBounds.y - heightBounds.x, 1.0);
	if (height01 <= 0.0 || height01 >= 1.0) return 0.0;

	float coverageThreshold = mix(0.50, 0.37, rainStrength);
	float localCoverage = smoothstep(coverageThreshold, 0.82, weather);
	if (localCoverage <= 0.0005) return 0.0;

	float cloudType = noiseCloudValue(weatherPosition * 0.63 + vec2(17.2, -9.4));
	cloudType = noiseSaturate(cloudType + localCoverage * 0.07 + rainStrength * 0.16);
	float cumulusAmount = smoothstep(0.38, 0.78, cloudType);
	// Do not terminate every cloud on one mathematical plane. A slow secondary
	// weather signal varies the underside while separate stratus and cumulus
	// profiles stop every formation from becoming the same cotton ball.
	float bottomNoise = noiseCloudValue(weatherPosition * 1.91 +
		vec2(-5.3, 12.8));
	float localBottom = mix(0.0, 0.075, bottomNoise) *
		mix(1.0, 0.58, rainStrength);
	float bottomProfile = smoothstep(localBottom,
		localBottom + 0.085, height01);
	float stratusTop = mix(0.34, 0.54,
		noiseCloudValue(weatherPosition * 1.27 + vec2(8.1, 4.6)));
	float cumulusTop = mix(0.56, 0.92,
		smoothstep(0.24, 0.82, cloudType));
	float stratusProfile = bottomProfile *
		(1.0 - smoothstep(max(stratusTop - 0.13, 0.16), stratusTop, height01));
	float cumulusProfile = bottomProfile *
		(1.0 - smoothstep(max(cumulusTop - 0.25, 0.22), cumulusTop, height01));
	float heightProfile = mix(stratusProfile, cumulusProfile, cumulusAmount);
	if (heightProfile <= 0.0005) return 0.0;

	// R is a connected Perlin-Worley base. GBA contain inverted Worley octaves.
	// The base is remapped by local coverage, matching the separation used by
	// production cloud systems: layout, base volume, then edge detail.
	vec3 basePosition = vec3(worldPosition.xz * 0.0058 + wind * 1.30,
		worldPosition.y * 0.0115);
	vec4 baseNoise = noiseCloudVolume(basePosition);
	float baseWorley = dot(baseNoise.gba, vec3(0.625, 0.250, 0.125));
	float baseShape = noiseSaturate((baseNoise.r - (1.0 - baseWorley)) /
		max(baseWorley, 0.08));
	// A softened coverage remap fills the middle of a weather system before
	// adding edge erosion. Directly using 1-coverage produced separate pearls.
	float cohesiveCoverage = pow(localCoverage, 0.72);
	float densityEdge = mix(0.76, 0.055, cohesiveCoverage);
	float baseDensity = smoothstep(densityEdge,
		min(densityEdge + 0.19, 1.0), baseShape) * heightProfile *
		smoothstep(0.015, 0.22, localCoverage);
	if (baseDensity <= 0.001) return 0.0;

	// A higher-frequency Worley lookup is evaluated only inside the conservative
	// base shape. It is band-limited by the caller when distance or ray step can
	// no longer resolve it. Its average (0.5), rather than zero, is used during
	// the fade so switching detail off does not dilate the silhouette.
	float erosionSignal = 0.5;
	if (detailAmount > 0.001) {
		vec4 detailNoise = noiseCloudVolume(basePosition * 2.35 +
			vec3(5.7, -3.1, 8.6));
		// Prefer the broad cell octave. The old equal-ish octave mixture created
		// dozens of detached bright beads around every otherwise coherent mass.
		float detailWorley = dot(detailNoise.gba, vec3(0.78, 0.17, 0.05));
		float detailedErosion = mix(1.0 - detailWorley, detailWorley,
			smoothstep(0.14, 0.52, height01));
		erosionSignal = mix(0.5, detailedErosion,
			noiseSaturate(detailAmount));
	}
	float erosionStrength = mix(0.165, 0.125, rainStrength);
	// Nubis-style erosion belongs primarily to the isosurface. Keeping the dense
	// core coherent removes the cellular foam look and also stabilizes lighting.
	float edgeErosion = mix(0.30, 1.0,
		1.0 - smoothstep(0.24, 0.72, baseDensity));
	float density = baseDensity - (1.0 - erosionSignal) *
		erosionStrength * edgeErosion;
	return smoothstep(0.0, 0.18, density) *
		mix(0.62, 1.0, smoothstep(0.04, 0.30, height01));
}

float noiseCloudDensityLod(vec3 worldPosition, vec2 wind, float rainStrength,
	float detailAmount) {
	float ignoredHeight;
	return noiseCloudDensityLodHeight(worldPosition, wind, rainStrength,
		detailAmount, ignoredHeight);
}

float noiseCloudDensity(vec3 worldPosition, vec2 wind, float rainStrength) {
	return noiseCloudDensityLod(worldPosition, wind, rainStrength, 1.0);
}

vec4 noiseCloudLayer(
	vec3 worldRay,
	vec3 cameraPosition,
	vec3 sunDirectionWorld,
	float daylight,
	float rainStrength,
	float frameTimeCounter,
	vec2 screenPosition,
	int frameIndex,
	bool reflectionPass
) {
	#if CLOUD_QUALITY == 0
	return vec4(0.0);
	#else
	if (abs(worldRay.y) < 0.004) return vec4(0.0);

	float firstDistance = (noiseCloudBottom - cameraPosition.y) / worldRay.y;
	float secondDistance = (noiseCloudTop - cameraPosition.y) / worldRay.y;
	float entryDistance = max(min(firstDistance, secondDistance), 0.0);
	float exitDistance = max(firstDistance, secondDistance);
	float traceDistance;
	#if CLOUD_QUALITY == 1
	traceDistance = 1200.0;
	#elif CLOUD_QUALITY == 2
	traceDistance = 1600.0;
	#else
	traceDistance = 2100.0;
	#endif
	if (reflectionPass) traceDistance = min(traceDistance, 1200.0);
	if (exitDistance <= entryDistance || entryDistance > traceDistance)
		return vec4(0.0);
	exitDistance = min(exitDistance, traceDistance);

	float raySpan = exitDistance - entryDistance;
	int sampleCount;
	if (reflectionPass) {
		// Reflections cover fewer pixels and receive their own spatial filtering.
		sampleCount = int(clamp(ceil(raySpan / 42.0), 6.0, 14.0));
	} else {
		#if CLOUD_QUALITY == 1
		// Half-resolution Fast is still only 1/4 as many cloud pixels. A slightly
		// denser, centered march is much cheaper than trying to reconstruct holes
		// left by ten heavily dithered samples after the fact.
		sampleCount = int(clamp(ceil(raySpan / 24.0), 12.0, 32.0));
		#elif CLOUD_QUALITY == 2
		sampleCount = int(clamp(ceil(raySpan / 18.0), 16.0, 44.0));
		#else
		// At grazing angles the old fixed 24 samples produced 30-40 block
		// strides, aliasing directly against the erosion field. Ultra follows the
		// 64-128 potential-sample production range and keeps the fine step below
		// ten blocks; empty-space stepping prevents paying the maximum everywhere.
		sampleCount = int(clamp(ceil(raySpan / 9.0), 32.0, 112.0));
		#endif
	}

	// Roughly 1.3 blocks/second: visibly flowing in normal play without racing
	// across the sky during timelapses.
	vec2 wind = frameTimeCounter * vec2(0.00550, 0.00180);
	float segmentLength = raySpan / float(sampleCount);
	// Interleaved gradient noise turns the remaining coherent ray-step bands
	// into sub-pixel error. The phase is deliberately stable for every quality:
	// rotating the Ultra ray origin changed the density isosurface every frame
	// and looked like two cloud textures fighting at the silhouette.
	float spatialPhase = fract(52.9829189 * fract(dot(floor(screenPosition),
		vec2(0.06711056, 0.00583715))));
	#if CLOUD_QUALITY == 1
	// A full-pixel random ray offset becomes a 2x2 stipple after half-resolution
	// reconstruction and is amplified again by RCAS. Fast uses the segment
	// midpoint; its lower-frequency density field does not need stochastic AA.
	float sliceJitter = 0.5;
	#elif CLOUD_QUALITY == 2
	float sliceJitter = mix(0.44, 0.56, spatialPhase);
	#else
	// Ultra already resolves the volume with a <=9 block step, so only a small
	// static offset is needed to decorrelate neighbouring rays.
	float sliceJitter = mix(0.40, 0.60, spatialPhase);
	#endif
	vec3 accumulatedColor = vec3(0.0);
	float accumulatedAlpha = 0.0;
	vec3 nightShadow = vec3(0.018, 0.024, 0.038);
	vec3 nightLight = vec3(0.120, 0.132, 0.170);
	vec3 dayShadow = mix(vec3(0.31, 0.38, 0.47), vec3(0.19, 0.22, 0.27), rainStrength);
	vec3 dayLight = mix(vec3(1.08, 1.035, 0.94), vec3(0.57, 0.60, 0.64), rainStrength);
	vec3 shadowColor = mix(nightShadow, dayShadow, daylight);
	vec3 lightColor = mix(nightLight, dayLight, daylight);
	// Approximate the already rendered background sky for per-sample aerial
	// perspective. Both reference packs suppress distant cloud contrast instead
	// of carrying a white material unchanged all the way to the clip distance.
	float skyElevation = smoothstep(-0.12, 0.48, worldRay.y);
	vec3 cloudAirZenith = mix(vec3(0.025, 0.030, 0.040),
		vec3(0.050, 0.23, 0.54), daylight);
	vec3 cloudAirHorizon = mix(vec3(0.090, 0.098, 0.112),
		vec3(0.62, 0.74, 0.84), daylight);
	vec3 cloudAirColor = mix(cloudAirHorizon, cloudAirZenith, skyElevation);
	float viewLightCosine = dot(worldRay, sunDirectionWorld);
	float forwardPhase = noiseCloudHenyeyGreenstein(viewLightCosine, 0.55);
	float backwardPhase = noiseCloudHenyeyGreenstein(viewLightCosine, -0.20);
	float phase = clamp(mix(forwardPhase, backwardPhase, 0.16) * 0.15,
		0.32, 1.65);
	float sunForward = pow(max(viewLightCosine, 0.0), 10.0) * daylight;
	float distanceAlongRay = entryDistance + sliceJitter * segmentLength;
	int consecutiveEmptySamples = 0;

	for (int sampleIndex = 0; sampleIndex < 112; ++sampleIndex) {
		if (distanceAlongRay >= exitDistance || accumulatedAlpha > 0.985) break;
		vec3 worldPosition = cameraPosition + worldRay * distanceAlongRay;
		// Fine erosion must disappear before its voxels or the march interval become
		// sub-pixel. Low therefore becomes simpler but genuinely stable instead of
		// trying to reconstruct Ultra detail from a quarter of the samples.
		float detailDistanceFade;
		#if CLOUD_QUALITY == 1
		detailDistanceFade = 1.0 - smoothstep(320.0, 620.0,
			distanceAlongRay);
		#elif CLOUD_QUALITY == 2
		detailDistanceFade = 1.0 - smoothstep(480.0, 980.0,
			distanceAlongRay);
		#else
		detailDistanceFade = 1.0 - smoothstep(760.0, 1650.0,
			distanceAlongRay);
		#endif
		float stepBandLimit = 1.0 - smoothstep(11.0, 23.0, segmentLength);
		float detailAmount = detailDistanceFade * stepBandLimit;
		float height01;
		float density = noiseCloudDensityLodHeight(worldPosition, wind,
			rainStrength, detailAmount, height01);
		if (density <= 0.001) {
			consecutiveEmptySamples++;
			// The base volume has much larger features than a fine integration
			// step. After three confirmed empty samples, cross empty regions twice
			// as fast; the next potential hit automatically restores fine steps.
			#if CLOUD_QUALITY == 1
			// With only 12-32 samples, a doubled step can jump over an entire
			// resolved cloud lobe and creates isolated sky-coloured pixels.
			distanceAlongRay += segmentLength;
			#else
			distanceAlongRay += segmentLength *
				(consecutiveEmptySamples >= 3 ? 2.0 : 1.0);
			#endif
			continue;
		}
		consecutiveEmptySamples = 0;

		// Probe once towards the sun. Dense cloud in front of this sample casts
		// an internal shadow; together with the vertical gradient this makes the
		// layer read as separate illuminated masses rather than white fog.
		vec3 lightProbeWorld = worldPosition + sunDirectionWorld * 18.0 +
			vec3(0.0, 2.0, 0.0);
		float nearBlockingDensity = noiseCloudDensityLod(lightProbeWorld, wind,
			rainStrength, detailAmount * 0.35);
		float farBlockingDensity = 0.0;
		#if CLOUD_QUALITY >= 3
		float coneAngle = spatialPhase * 6.28318530718 + float(sampleIndex) * 2.39996323;
		vec3 coneOffset = vec3(cos(coneAngle), 0.35, sin(coneAngle)) * 7.0;
		vec3 farLightProbeWorld = worldPosition + sunDirectionWorld * 58.0 +
			coneOffset;
		farBlockingDensity = noiseCloudDensityLod(farLightProbeWorld, wind,
			rainStrength, detailAmount * 0.20);
		#endif
		float lightOpticalDepth = nearBlockingDensity * 1.35 +
			farBlockingDensity * 0.90;
		float lightVisibility = exp(-lightOpticalDepth *
			mix(1.75, 2.45, rainStrength));
		float topLighting = smoothstep(0.10, 0.92, height01);
		float baseShadow = 1.0 - smoothstep(0.02, 0.34, height01);
		float powder = 1.0 - exp(-density * 2.40);
		float ambientLighting = 0.16 + topLighting * 0.25 +
			(1.0 - density) * 0.08;
		float directLighting = lightVisibility * (0.43 + powder * 0.27) * phase;
		float lighting = ambientLighting + directLighting;
		lighting -= baseShadow * 0.14 +
			density * (1.0 - lightVisibility) * 0.12;
		lighting = clamp(lighting, 0.0, 1.0);
		vec3 sampleColor = mix(shadowColor, lightColor, lighting);
		float silverLining = sunForward * pow(1.0 - density, 2.0) * lightVisibility;
		sampleColor += vec3(1.0, 0.77, 0.48) * silverLining * 0.34;
		float airExtinction = 0.00052 * mix(1.55, 1.0,
			noiseSaturate(abs(worldRay.y)));
		float aerialPerspective = 1.0 - exp(-distanceAlongRay * airExtinction);
		sampleColor = mix(sampleColor, cloudAirColor,
			aerialPerspective * 0.78);

		float sampleAlpha = 1.0 - exp(-density * segmentLength * mix(0.060, 0.082, rainStrength));
		sampleAlpha *= 1.0 - aerialPerspective * 0.18;
		float visibleAlpha = (1.0 - accumulatedAlpha) * sampleAlpha;
		accumulatedColor += sampleColor * visibleAlpha;
		accumulatedAlpha += visibleAlpha;
		distanceAlongRay += segmentLength;
	}

	float horizonFade = smoothstep(0.012, 0.085, abs(worldRay.y));
	float farFade = 1.0 - smoothstep(traceDistance * 0.80,
		traceDistance, entryDistance);
	accumulatedAlpha *= horizonFade * farFade;
	if (accumulatedAlpha <= 0.0001) return vec4(0.0);
	return vec4(accumulatedColor / max(accumulatedAlpha, 0.0001), accumulatedAlpha);
	#endif
}

#endif
