/********************************************************
    © 2020 Continuum Graphics LLC. All Rights Reserved
 ********************************************************/

#if !defined _PHYSICALSKY_
#define _PHYSICALSKY_

#include "/../UserLib/Sky/SkyConstants.glsl"

/*******************************************************************************
- Phases
******************************************************************************/

#define phaseg0 0.25

float sky_rayleighPhase(float cosTheta) {
	const float atten = 3.0 / (16.0 * PI);
	return atten * (1.0 + cosTheta * cosTheta);
}

float sky_miePhase(float cosTheta, const float g) {
	const float gg = g * g;
	return (1.0 - gg) / (4.0 * PI * pow(1.0 + gg - 2.0 * g * cosTheta, 1.5));
}

vec2 sky_phase(float cosTheta, const float g) {
	return vec2(sky_rayleighPhase(cosTheta), sky_miePhase(cosTheta, g));
}

/*******************************************************************************
- Sky Functions
******************************************************************************/

float calculateStars(vec3 p, vec3 wSunVector) {
	p = p * RotateMatrix(vec3(0.0, 1.0, 0.0), wSunVector);

	const float res = 2048.0;
	const float starCoverage = 0.01;

	vec3 rp = p * (0.15 * res);

	vec3 id  = floor(rp);
	vec3 q   = id - rp + 0.5;
	vec3 rn  = hash33(id);
	float c  = 1.0 - smoothstep(0.0, 0.5, length(q));
		c *= step(rn.x, starCoverage * 0.1);

	return c * 16.0;
}

float CalculateSunSpot(float VdotL) {
	const float sunAngularSize = 0.545;
	const float sunRadius = radians(sunAngularSize);
	const float cosSunRadius = cos(sunRadius);
	const float sunLuminance = 1.0 / ((1.0 - cosSunRadius) * PI);

	return fstep(cosSunRadius, VdotL) * sunLuminance;
}

//Temporary moon disk untill we added a moon texture
float CalculateMoonSpot(float VdotL) {
	const float moonAngularSize = 0.545;
	const float moonRadius = radians(moonAngularSize);
	const float cosMoonRadius = cos(moonRadius);
	const float moonLuminance = 1.0 / ((1.0 - cosMoonRadius) * PI);

	return fstep(cosMoonRadius, VdotL) * moonLuminance;
}

/*******************************************************************************
- Sky Functions
******************************************************************************/

// No intersection if returned y component is < 0.0
vec2 rsi(vec3 position, vec3 direction, const float radius) {
	float PoD = dot(position, direction);
	const float radiusSquared = radius * radius;

	float delta = PoD * PoD + radiusSquared - dot(position, position);
	if (delta < 0.0) return vec2(-1.0);
		delta = sqrt(delta);

	return -PoD + vec2(-delta, delta);
}

vec3 sky_atmosphereDensity(float centerDistance) {
	vec2 rayleighMie = exp(centerDistance * -atmosphereInverseScaleHeights + atmosphereScaledPlanetRadius);

	float ozone = exp(-max(0.0, (35000.0 - centerDistance) - atmospherePlanetRadius) * 0.0002)
				* exp(-max(0.0, (centerDistance - 35000.0) - atmospherePlanetRadius) / 15000.0);

	return vec3(rayleighMie, ozone);
}

// I don't know if "thickness" is the right word, using it because Jodie uses it for that and I can't think of (or find) anything better.
vec3 sky_atmosphereThickness(vec3 position, vec3 direction, float rayLength, const float steps) {
	float stepSize  = rayLength / steps;
	vec3  increment = direction * stepSize;
	position += increment * 0.5;

	vec3 thickness = vec3(0.0);
	for (float i = 0.0; i < steps; ++i, position += increment) {
		thickness += sky_atmosphereDensity(length(position));
	}

	return thickness * stepSize;
}
vec3 sky_atmosphereThickness(vec3 position, vec3 direction, const float steps) {
	float rayLength = dot(position, direction);
		rayLength = sqrt(rayLength * rayLength + atmosphereRadiusSquared - dot(position, position)) - rayLength;

	return sky_atmosphereThickness(position, direction, rayLength, steps);
}

vec3 sky_atmosphereOpticalDepth(vec3 position, vec3 direction, float rayLength, const float steps) {
	return atmosphereAttenuationCoefficients * sky_atmosphereThickness(position, direction, rayLength, steps);
}
vec3 sky_atmosphereOpticalDepth(vec3 position, vec3 direction, const float steps) {
	return atmosphereAttenuationCoefficients * sky_atmosphereThickness(position, direction, steps);
}

vec3 sky_atmosphereTransmit(vec3 position, vec3 direction, const float steps) {
	return exp2(-sky_atmosphereOpticalDepth(position, direction, steps) * rLOG2);
}

vec3 GetSunColorZom() {
	return sky_phase(1.0, atmosphereMieG).x * sky_phase(1.0, atmosphereMieG).y * sky_atmosphereTransmit(vec3(0.0, atmospherePlanetRadius, 0.0), mat3(gbufferModelViewInverse) * (sunPositionVert * 0.01), 8) * skySunColor;
}

vec3 GetMoonColorZom() {
	return sky_atmosphereTransmit(vec3(0.0, atmospherePlanetRadius, 0.0), mat3(gbufferModelViewInverse) * (-sunPositionVert * 0.01), 8) * skyMoonColor;
}

vec3 sky_atmosphere(vec3 background, vec3 viewVector, vec3 upVector, vec3 sunVector, vec3 moonVector, vec3 sunIlluminance, vec3 moonIlluminance, const int iSteps, inout vec3 transmit) {
	//const int iSteps = 25; // For very high quality: 50 is enough, could get away with less if mie scale height was lower
	const int jSteps = 3;  // For very high quality: 10 is good, can probably get away with less

	vec3 viewPosition = (atmospherePlanetRadius + eyeAltitude) * upVector;

	vec2 aid = rsi(viewPosition, viewVector, atmosphereRadius);
	if (aid.y < 0.0) return background;
	vec2 pid = rsi(viewPosition, viewVector, atmospherePlanetRadius * 0.998);

	bool pi = pid.y >= 0.0;

	vec2 sd = vec2((pi && pid.x < 0.0) ? pid.y : max(aid.x, 0.0), (pi && pid.x > 0.0) ? pid.x : aid.y);

	float stepSize  = (sd.y - sd.x) / iSteps;
	vec3  increment = viewVector * stepSize;
	vec3  position  = viewVector * sd.x + (increment * 0.3 + viewPosition);

	vec2 phaseSun  = sky_phase(dot(viewVector, sunVector ), atmosphereMieG);
	vec2 phaseMoon = sky_phase(dot(viewVector, moonVector), atmosphereMieG);

	vec3 scatteringSun  = vec3(0.0);
	vec3 scatteringMoon = vec3(0.0);
	vec3 scatteringAmbient = vec3(0.0);
		transmit = vec3(1.0);

	for (int i = 0; i < iSteps; ++i, position += increment) {
		vec3 density          = sky_atmosphereDensity(length(position));
		if (density.y > 1e35) break;
		vec3 stepAirmass      = density * stepSize;
		vec3 stepOpticalDepth = atmosphereAttenuationCoefficients * stepAirmass;

		vec3 stepTransmit       = exp2(-stepOpticalDepth * rLOG2);
		vec3 stepTransmittedFraction = clamp((stepTransmit - 1.0) / -stepOpticalDepth, 0.0, 1.0);
		vec3 stepScatteringVisible   = transmit * stepTransmittedFraction;

		scatteringSun  += (atmosphereScatteringCoefficients * (stepAirmass.xy * phaseSun )) * stepScatteringVisible * sky_atmosphereTransmit(position, sunVector,  jSteps);
		scatteringMoon += (atmosphereScatteringCoefficients * (stepAirmass.xy * phaseMoon)) * stepScatteringVisible * sky_atmosphereTransmit(position, moonVector, jSteps);

		transmit  *= stepTransmit;
	}

	vec3 scattering = scatteringSun * sunIlluminance + scatteringMoon * moonIlluminance;

	return (!pi ? background * transmit : vec3(0.0)) + mix(scattering, vec3(dot(scattering, vec3(0.2125, 0.7154, 0.0721))), wetness);
}

#endif
