#include "/Lib/Settings.glsl"
#include "/Lib/Utilities.glsl"
// MinecraftPT — Precomputed Atmosphere (analytic)
// A compact analytic Rayleigh/Mie sky model with precomputed transmittance.
// Provides sun/sky irradiance for the scene and sky radiance for the panorama.

#ifndef PRECOMPUTED_ATMOSPHERE_GLSL
#define PRECOMPUTED_ATMOSPHERE_GLSL

// Earth-like atmosphere constants (scaled)
const float atmosphereModel_bottom_radius = 6360.0;
const float atmosphereModel_top_radius = 6420.0;

const vec3 wavelengths = vec3(680.0, 550.0, 440.0); // nm
const vec3 betaRayleigh = vec3(5.802, 13.558, 33.1) * 1e-6;
const float betaMie = 3.996e-6;
const float mieG = 0.8;

vec3 RayleighPhase(float cosTheta){
	return 3.0 / (16.0 * 3.14159265) * (1.0 + cosTheta * cosTheta);
}

vec3 MiePhase(float cosTheta){
	float g = mieG;
	float g2 = g * g;
	return 3.0 / (8.0 * 3.14159265) * ((1.0 - g2) * (1.0 + cosTheta * cosTheta)) /
		((2.0 + g2) * pow(1.0 + g2 - 2.0 * g * cosTheta, 1.5));
}

// Ground intersection for ray from camera
float GroundIntersection(vec3 pos, vec3 dir){
	float b = dot(pos, dir);
	float c = dot(pos, pos) - atmosphereModel_bottom_radius * atmosphereModel_bottom_radius;
	float h = b * b - c;
	if (h > 0.0){
		float t = -b - sqrt(h);
		if (t > 0.0) return t;
	}
	return -1.0;
}

float AtmosphereIntersection(vec3 pos, vec3 dir){
	float b = dot(pos, dir);
	float c = dot(pos, pos) - atmosphereModel_top_radius * atmosphereModel_top_radius;
	float h = b * b - c;
	if (h > 0.0){
		float t = -b + sqrt(h);
		if (t > 0.0) return t;
	}
	return -1.0;
}

// Optical depth to space (transmittance approximation)
vec3 OpticalDepthToSpace(vec3 pos, vec3 dir){
	float tMax = AtmosphereIntersection(pos, dir);
	if (tMax < 0.0) return vec3(1e10);

	const int steps = 16;
	float dt = tMax / float(steps);
	vec3 opticalDepth = vec3(0.0);

	vec3 p = pos;
	for (int i = 0; i < steps; i++){
		float h = length(p) - atmosphereModel_bottom_radius;
		vec3 density = exp(-h / 8.0) * betaRayleigh + exp(-h / 1.2) * betaMie;
		opticalDepth += density * dt;
		p += dir * dt;
	}

	return opticalDepth;
}

// Sun/sky irradiance at a camera position (simplified single scattering)
void GetSunAndSkyIrradiance(vec3 camera, vec3 sunDir, vec3 moonDir,
	out vec3 colorSunlight, out vec3 colorMoonlight,
	out vec3 colorSunSkylight, out vec3 colorMoonSkylight){

	vec3 sunColor = vec3(1.0, 0.98, 0.92);
	vec3 moonColor = vec3(0.5, 0.6, 0.8) * 0.1;

	float sunElevation = sunDir.y;
	float moonElevation = moonDir.y;

	float sunStrength = curve(saturate(sunElevation * 15.0 + 0.5));
	float moonStrength = curve(saturate(-sunElevation * 5.0 + 0.5));

	// Sun disk color by elevation (sunset tint)
	vec3 sunsetTint = mix(vec3(1.0, 0.3, 0.1), vec3(1.0), curve(saturate(sunElevation * 3.0)));
	sunColor *= sunsetTint;

	colorSunlight = sunColor * sunStrength * SUNLIGHT_INTENSITY;
	colorMoonlight = moonColor * moonStrength;

	// Sky irradiance: ambient hemisphere light
	vec3 skyColorDay = vec3(0.55, 0.75, 1.0) * sunStrength * 0.8;
	vec3 skyColorNight = vec3(0.05, 0.08, 0.15) * 0.3;
	vec3 skyColorSunset = vec3(1.0, 0.5, 0.3) * curve(saturate(1.0 - abs(sunElevation) * 4.0)) * 0.3;

	colorSunSkylight = skyColorDay + skyColorSunset;
	colorMoonSkylight = skyColorNight;
}

// Sky radiance for a direction (used for the sky panorama and reflections)
vec3 GetSkyRadiance(vec3 dir, vec3 sunDir){
	dir = normalize(dir);
	sunDir = normalize(sunDir);

	float cosTheta = dot(dir, sunDir);

	// Horizon / elevation response
	float elevation = dir.y;

	// Day sky
	vec3 skyDay = vec3(0.4, 0.62, 0.9) * 0.6;
	skyDay = mix(skyDay, vec3(0.9, 0.95, 1.0), pow(saturate(elevation), 0.6));
	skyDay *= 1.0 - exp(-max(elevation, 0.0) * 3.0);

	// Rayleigh scattering glow around sun
	vec3 rayleigh = RayleighPhase(cosTheta) * betaRayleigh * 1.2;

	// Mie scattering (sun disk)
	float mie = MiePhase(cosTheta) * betaMie * 20.0;

	// Sun disk itself
	float sunDisk = exp(-(1.0 - cosTheta) / (2.0 * SUN_ANGULAR_RADIUS * SUN_ANGULAR_RADIUS));

	vec3 sunColor = vec3(1.0, 0.97, 0.9) * SUNLIGHT_INTENSITY;

	// Sunset tint
	float sunElevation = sunDir.y;
	sunColor *= mix(vec3(1.0, 0.3, 0.1), vec3(1.0), curve(saturate(sunElevation * 3.0)));

	float sunVisible = curve(saturate(sunElevation * 15.0 + 0.5));

	vec3 radiance = skyDay;
	radiance += sunColor * (rayleigh + mie) * sunVisible;
	radiance += sunColor * sunDisk * sunVisible * 80.0;

	// Night sky
	float night = 1.0 - curve(saturate(sunElevation * 10.0 + 0.5));
	radiance += vec3(0.02, 0.03, 0.06) * night;

	// Horizon haze
	float horizon = exp(-abs(elevation) * 4.0);
	radiance += vec3(0.9, 0.7, 0.5) * horizon * sunVisible * 0.1;

	return radiance;
}

#endif