// UTILITIES ///////////////////////////////
float ClampCosine(float mu) 
{
	return clamp(mu, float(-1.0), float(1.0));
}

float ClampDistance(float d) 
{
	return max(d, 0.0);
}

float ClampRadius(in AtmosphereParameters atmosphere, float r) 
{
	return clamp(r, atmosphere._bottom_radius, atmosphere._top_radius);
}

float SafeSqrt(float a) 
{
	return sqrt(max(a, 0.0));
}

// COMMON //////////////////////////////////////
float DistanceToTopAtmosphereBoundary(in AtmosphereParameters atmosphere, float r, float mu) 
{
	//assert(r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	float discriminant = r * r * (mu * mu - 1.0) + atmosphere._top_radius * atmosphere._top_radius;
	return ClampDistance(-r * mu + SafeSqrt(discriminant));
}

float DistanceToBottomAtmosphereBoundary(in AtmosphereParameters atmosphere, float r, float mu) 
{
	//assert(r >= atmosphere.bottom_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	float discriminant = r * r * (mu * mu - 1.0) + atmosphere._bottom_radius * atmosphere._bottom_radius;
	return ClampDistance(-r * mu - SafeSqrt(discriminant));
}

bool RayIntersectsGround( in AtmosphereParameters atmosphere, float r, float mu ) 
{
	//assert(r >= atmosphere.bottom_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	return mu < 0.0 && r * r * (mu * mu - 1.0) + atmosphere._bottom_radius * atmosphere._bottom_radius >= 0.0; //* m2; //??
}

float GetLayerDensity(in DensityProfileLayer layer, float altitude) 
{
	float density = layer._exp_term * exp(layer._exp_scale * altitude) + layer._linear_term * altitude + layer._constant_term;
	return clamp(density, 0.0, 1.0);
}

float GetProfileDensity(in DensityProfile profile, float altitude) 
{
	return altitude < profile._layers[0]._width ? GetLayerDensity(profile._layers[0], altitude) : GetLayerDensity(profile._layers[1], altitude);
}

float GetUnitRangeFromTextureCoord(float u, int texture_size) 
{
	return (u - 0.5 / float(texture_size)) / (1.0 - 1.0 / float(texture_size));
}

float GetTextureCoordFromUnitRange(float x, int texture_size) 
{
	return 0.5 / float(texture_size) + x * (1.0 - 1.0 / float(texture_size));
}

/* 
-------------------------------------------------- 

TRANSMITTANCE

-----------------------------------------------------
*/

// COMPUTE /////////////////////////////////////////

float ComputeOpticalLengthToTopAtmosphereBoundary(in AtmosphereParameters atmosphere, in DensityProfile profile, float r, float mu) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
  
	// Number of intervals for the numerical integration.
	const int SAMPLE_COUNT = 500;
 
	// The integration step, i.e. the length of each integration interval.
	float dx = DistanceToTopAtmosphereBoundary(atmosphere, r, mu) / float(SAMPLE_COUNT);
	
	// Integration loop.
	float result = 0.0;// * m;
  
	for (int i = 0; i <= SAMPLE_COUNT; ++i) 
	{
		float d_i = float(i) * dx;
    
		// Distance between the current sample point and the planet center.
		float r_i = sqrt(d_i * d_i + 2.0 * r * mu * d_i + r * r);
    
		// Number density at the current sample point (divided by the number density
		// at the bottom of the atmosphere, yielding a dimensionless number).
		float y_i = GetProfileDensity(profile, r_i - atmosphere._bottom_radius);
    
		// Sample weight (from the trapezoidal rule).
		float weight_i = (i == 0 || i == SAMPLE_COUNT ? 0.5 : 1.0);
    
		result += y_i * weight_i * dx;
	}
	
	return result;
}

vec3 ComputeTransmittanceToTopAtmosphereBoundary( in AtmosphereParameters atmosphere, float r, float mu) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	float opticalLengthRay=ComputeOpticalLengthToTopAtmosphereBoundary(atmosphere, atmosphere._rayleigh_density, r, mu);
	float opticalLengthMie=ComputeOpticalLengthToTopAtmosphereBoundary(atmosphere, atmosphere._mie_density, r, mu);
	float opticalLengthAbs=ComputeOpticalLengthToTopAtmosphereBoundary(atmosphere, atmosphere._absorption_density, r, mu);
	
	return exp( -(atmosphere._rayleigh_scattering * opticalLengthRay + atmosphere._mie_extinction * opticalLengthMie + atmosphere._absorption_extinction * opticalLengthAbs) );
}

/////// PRECOMPUTE ///////////////////////////////////////////

void GetRMuFromTransmittanceTextureUv(in AtmosphereParameters atmosphere, in vec2 uv, out float r, out float mu) 
{
	//assert(uv.x >= 0.0 && uv.x <= 1.0);
	//assert(uv.y >= 0.0 && uv.y <= 1.0);
  
	float x_mu = GetUnitRangeFromTextureCoord(uv.x, TRANSMITTANCE_TEXTURE_WIDTH);
	float x_r = GetUnitRangeFromTextureCoord(uv.y, TRANSMITTANCE_TEXTURE_HEIGHT);
  
	// Distance to top atmosphere boundary for a horizontal ray at ground level.
	float H = sqrt(atmosphere._top_radius * atmosphere._top_radius - atmosphere._bottom_radius * atmosphere._bottom_radius);
	
	// Distance to the horizon, from which we can compute r:
	float rho = H * x_r;
	r = sqrt(rho * rho + atmosphere._bottom_radius * atmosphere._bottom_radius);
  
	// Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
	// and maximum values over all mu - obtained for (r,1) and (r,mu_horizon) -
	// from which we can recover mu:
	float d_min = atmosphere._top_radius - r;
	float d_max = rho + H;
	float d = d_min + x_mu * (d_max - d_min);
	
	mu = d == 0.0 /*m*/? 1.0 : (H * H - rho * rho - d * d) / (2.0 * r * d);
	mu = ClampCosine(mu);
}

vec3 ComputeTransmittanceToTopAtmosphereBoundaryTexture(in AtmosphereParameters atmosphere, in vec2 frag_coord) 
{
  const vec2 TRANSMITTANCE_TEXTURE_SIZE = vec2(TRANSMITTANCE_TEXTURE_WIDTH, TRANSMITTANCE_TEXTURE_HEIGHT);
  
  float r;
  float mu;
  
  GetRMuFromTransmittanceTextureUv(atmosphere, frag_coord / TRANSMITTANCE_TEXTURE_SIZE, r, mu);
  return ComputeTransmittanceToTopAtmosphereBoundary(atmosphere, r, mu);
}

// LOOKUP ///////////////////////////////////////////////////////////////

vec2 GetTransmittanceTextureUvFromRMu(in AtmosphereParameters atmosphere, float r, float mu) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
  
	// Distance to top atmosphere boundary for a horizontal ray at ground level.
	float H = sqrt(atmosphere._top_radius * atmosphere._top_radius - atmosphere._bottom_radius * atmosphere._bottom_radius);
  
	// Distance to the horizon.
	float rho = SafeSqrt(r * r - atmosphere._bottom_radius * atmosphere._bottom_radius);
  
	// Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
	// and maximum values over all mu - obtained for (r,1) and (r,mu_horizon).
	float d = DistanceToTopAtmosphereBoundary(atmosphere, r, mu);
	float d_min = atmosphere._top_radius - r;
	float d_max = rho + H;
	float x_mu = (d - d_min) / (d_max - d_min);
	float x_r = rho / H;
  
	return vec2(GetTextureCoordFromUnitRange(x_mu, TRANSMITTANCE_TEXTURE_WIDTH), GetTextureCoordFromUnitRange(x_r, TRANSMITTANCE_TEXTURE_HEIGHT));
}

vec3 GetTransmittanceToTopAtmosphereBoundary(in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, float r, float mu) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
  
	vec2 uv = GetTransmittanceTextureUvFromRMu(atmosphere, r, mu);
  
	return vec3(texture(transmittance_texture, uv));
}

vec3 GetTransmittance( in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, float r, float mu, float d, bool ray_r_mu_intersects_ground) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	//assert(d >= 0.0 * m);

	float r_d = ClampRadius(atmosphere, sqrt(d * d + 2.0 * r * mu * d + r * r));
	float mu_d = ClampCosine((r * mu + d) / r_d);

	if (ray_r_mu_intersects_ground) 
		return min(GetTransmittanceToTopAtmosphereBoundary(atmosphere, transmittance_texture, r_d, -mu_d)/GetTransmittanceToTopAtmosphereBoundary(atmosphere, transmittance_texture, r, -mu),vec3(1.0));
	else 
		return min(GetTransmittanceToTopAtmosphereBoundary(atmosphere, transmittance_texture, r, mu)/GetTransmittanceToTopAtmosphereBoundary(atmosphere, transmittance_texture, r_d, mu_d),vec3(1.0));
}

vec3 GetTransmittanceToSun( in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, float r, float mu_s) 
{
	float sin_theta_h = atmosphere._bottom_radius / r;
	float cos_theta_h = -sqrt(max(1.0 - sin_theta_h * sin_theta_h, 0.0));
  
	return GetTransmittanceToTopAtmosphereBoundary(atmosphere, transmittance_texture, r, mu_s) * 
	smoothstep(-sin_theta_h * atmosphere._sun_angular_radius /*/ rad*/, sin_theta_h * atmosphere._sun_angular_radius /*/ rad*/, mu_s - cos_theta_h);
}


/* 
-------------------------------------------------- 

SINGLE SCATTERING

-----------------------------------------------------
*/

// COMPUTE ///////////////////////////////////////////////////////////////////

void ComputeSingleScatteringIntegrand( in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, float r, float mu, float mu_s, float nu, float d, bool ray_r_mu_intersects_ground, out vec3 rayleigh, out vec3 mie) 
{
	float r_d = ClampRadius(atmosphere, sqrt(d * d + 2.0 * r * mu * d + r * r));
	float mu_s_d = ClampCosine((r * mu_s + d * nu) / r_d);
  
	vec3 transmittance =  	GetTransmittance(atmosphere, transmittance_texture, r, mu, d, ray_r_mu_intersects_ground) *  
							GetTransmittanceToSun(atmosphere, transmittance_texture, r_d, mu_s_d);
							
	rayleigh = transmittance * GetProfileDensity(atmosphere._rayleigh_density, r_d - atmosphere._bottom_radius);
	mie = transmittance * GetProfileDensity(atmosphere._mie_density, r_d - atmosphere._bottom_radius);
}

float DistanceToNearestAtmosphereBoundary(in AtmosphereParameters atmosphere, float r, float mu, bool ray_r_mu_intersects_ground) 
{
	if (ray_r_mu_intersects_ground)
		return DistanceToBottomAtmosphereBoundary(atmosphere, r, mu);
	else
		return DistanceToTopAtmosphereBoundary(atmosphere, r, mu);
}

void ComputeSingleScattering(in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, float r, float mu, float mu_s, float nu, bool ray_r_mu_intersects_ground, out vec3 rayleigh, out vec3 mie) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	//assert(mu_s >= -1.0 && mu_s <= 1.0);
	//assert(nu >= -1.0 && nu <= 1.0);

	// Number of intervals for the numerical integration.
	const int SAMPLE_COUNT = 50;
  
	// The integration step, i.e. the length of each integration interval.
	float dx = DistanceToNearestAtmosphereBoundary(atmosphere, r, mu, ray_r_mu_intersects_ground) / float(SAMPLE_COUNT);
	
	// Integration loop.
	vec3 rayleigh_sum = vec3(0.0);
	vec3 mie_sum = vec3(0.0);
	
	for (int i = 0; i <= SAMPLE_COUNT; ++i) 
	{
		float d_i = float(i) * dx;
    
		// The Rayleigh and Mie single scattering at the current sample point.
		vec3 rayleigh_i;
		vec3 mie_i;
    
		ComputeSingleScatteringIntegrand(atmosphere, transmittance_texture, r, mu, mu_s, nu, d_i, ray_r_mu_intersects_ground, rayleigh_i, mie_i);
    
		// Sample weight (from the trapezoidal rule).
		float weight_i = (i == 0 || i == SAMPLE_COUNT) ? 0.5 : 1.0;
		
		rayleigh_sum += rayleigh_i * weight_i;
		mie_sum += mie_i * weight_i;
	}
  
	rayleigh = rayleigh_sum * dx * atmosphere._solar_irradiance * atmosphere._rayleigh_scattering;
	mie = mie_sum * dx * atmosphere._solar_irradiance * atmosphere._mie_scattering;
}

float RayleighPhaseFunction(float nu) 
{
	float k = 3.0 / (16.0 * PI /*sr*/);
	return k * (1.0 + nu * nu);
}

float MiePhaseFunction(float g, float nu) 
{
	float k = 3.0 / (8.0 * PI /*sr*/) * (1.0 - g * g) / (2.0 + g * g);
	return k * (1.0 + nu * nu) / pow(1.0 + g * g - 2.0 * g * nu, 1.5);
}

// PRECOMPUTE //////////////////////////////////////////////////////

void GetRMuMuSNuFromScatteringTextureUvwz(in AtmosphereParameters atmosphere, in vec4 uvwz, out float r, out float mu, out float mu_s, out float nu, out bool ray_r_mu_intersects_ground) 
{
	//assert(uvwz.x >= 0.0 && uvwz.x <= 1.0);
	//assert(uvwz.y >= 0.0 && uvwz.y <= 1.0);
	//assert(uvwz.z >= 0.0 && uvwz.z <= 1.0);
	//assert(uvwz.w >= 0.0 && uvwz.w <= 1.0);

	// Distance to top atmosphere boundary for a horizontal ray at ground level.
	float H = sqrt(atmosphere._top_radius * atmosphere._top_radius - atmosphere._bottom_radius * atmosphere._bottom_radius);
  
	// Distance to the horizon.
	float rho =  H * GetUnitRangeFromTextureCoord(uvwz.w, SCATTERING_TEXTURE_R_SIZE);
	r = sqrt(rho * rho + atmosphere._bottom_radius * atmosphere._bottom_radius);

	if (uvwz.z < 0.5) 
	{
		// Distance to the ground for the ray (r,mu), and its minimum and maximum
		// values over all mu - obtained for (r,-1) and (r,mu_horizon) - from which
		// we can recover mu:
		float d_min = r - atmosphere._bottom_radius;
		float d_max = rho;
		float d = d_min + (d_max - d_min) * GetUnitRangeFromTextureCoord(1.0 - 2.0 * uvwz.z, SCATTERING_TEXTURE_MU_SIZE / 2);
		mu = d == 0.0 /*m*/ ? float(-1.0) :
        ClampCosine(-(rho * rho + d * d) / (2.0 * r * d));
		ray_r_mu_intersects_ground = true;
	} 
	else 
	{
		// Distance to the top atmosphere boundary for the ray (r,mu), and its
		// minimum and maximum values over all mu - obtained for (r,1) and
		// (r,mu_horizon) - from which we can recover mu:
		float d_min = atmosphere._top_radius - r;
		float d_max = rho + H;
		float d = d_min + (d_max - d_min) * GetUnitRangeFromTextureCoord(2.0 * uvwz.z - 1.0, SCATTERING_TEXTURE_MU_SIZE / 2);
		mu = d == 0.0 /*m*/ ? float(1.0) :
        ClampCosine((H * H - rho * rho - d * d) / (2.0 * r * d));
		ray_r_mu_intersects_ground = false;
	}

	float x_mu_s = GetUnitRangeFromTextureCoord(uvwz.y, SCATTERING_TEXTURE_MU_S_SIZE);
	float d_min = atmosphere._top_radius - atmosphere._bottom_radius;
	float d_max = H;
	
	float A =-2.0 * atmosphere._mu_s_min * atmosphere._bottom_radius / (d_max - d_min);
	float a = (A - x_mu_s * A) / (1.0 + x_mu_s * A);
	float d = d_min + min(a, A) * (d_max - d_min);
	mu_s = d == 0.0 /*m*/ ? float(1.0) :
	ClampCosine((H * H - d * d) / (2.0 * atmosphere._bottom_radius * d));
	nu = ClampCosine(uvwz.x * 2.0 - 1.0);
}

void GetRMuMuSNuFromScatteringTextureFragCoord(in AtmosphereParameters atmosphere, in vec3 frag_coord, out float r, out float mu, out float mu_s, out float nu, out bool ray_r_mu_intersects_ground) 
{
	const vec4 SCATTERING_TEXTURE_SIZE = vec4(SCATTERING_TEXTURE_NU_SIZE - 1, SCATTERING_TEXTURE_MU_S_SIZE, SCATTERING_TEXTURE_MU_SIZE, SCATTERING_TEXTURE_R_SIZE);
  
	float frag_coord_nu = floor(frag_coord.x / float(SCATTERING_TEXTURE_MU_S_SIZE));
	float frag_coord_mu_s = mod(frag_coord.x, float(SCATTERING_TEXTURE_MU_S_SIZE));
	vec4 uvwz = vec4(frag_coord_nu, frag_coord_mu_s, frag_coord.y, frag_coord.z) / SCATTERING_TEXTURE_SIZE;
	
	GetRMuMuSNuFromScatteringTextureUvwz(atmosphere, uvwz, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
  
	// Clamp nu to its valid range of values, given mu and mu_s.
	nu = clamp(nu, mu * mu_s - sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s)), mu * mu_s + sqrt((1.0 - mu * mu) * (1.0 - mu_s * mu_s)));
}

void ComputeSingleScatteringTexture(in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, in vec3 frag_coord, out vec3 rayleigh, out vec3 mie) 
{
	float r;
	float mu;
	float mu_s;
	float nu;
	bool ray_r_mu_intersects_ground;
  
	GetRMuMuSNuFromScatteringTextureFragCoord(atmosphere, frag_coord, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
  
	ComputeSingleScattering(atmosphere, transmittance_texture, r, mu, mu_s, nu, ray_r_mu_intersects_ground, rayleigh, mie);
}

// LOOKUP //////////////////////////////////////////////////////

vec4 GetScatteringTextureUvwzFromRMuMuSNu(in AtmosphereParameters atmosphere, float r, float mu, float mu_s, float nu, bool ray_r_mu_intersects_ground) 
{
	// assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	// assert(mu >= -1.0 && mu <= 1.0);
	// assert(mu_s >= -1.0 && mu_s <= 1.0);
	// assert(nu >= -1.0 && nu <= 1.0);

	// Distance to top atmosphere boundary for a horizontal ray at ground level.
	float H = sqrt(atmosphere._top_radius * atmosphere._top_radius - atmosphere._bottom_radius * atmosphere._bottom_radius);
  
	// Distance to the horizon.
	float rho = SafeSqrt(r * r - atmosphere._bottom_radius * atmosphere._bottom_radius);
	float u_r = GetTextureCoordFromUnitRange(rho / H, SCATTERING_TEXTURE_R_SIZE);

	// Discriminant of the quadratic equation for the intersections of the ray
	// (r,mu) with the ground (see RayIntersectsGround).
	float r_mu = r * mu;
	float discriminant = r_mu * r_mu - r * r + atmosphere._bottom_radius * atmosphere._bottom_radius;
  
	float u_mu;
	if (ray_r_mu_intersects_ground) 
	{
		// Distance to the ground for the ray (r,mu), and its minimum and maximum
		// values over all mu - obtained for (r,-1) and (r,mu_horizon).
		float d = -r_mu - SafeSqrt(discriminant);
		float d_min = r - atmosphere._bottom_radius;
		float d_max = rho;
		u_mu = 0.5 - 0.5 * GetTextureCoordFromUnitRange(d_max == d_min ? 0.0 : (d - d_min) / (d_max - d_min), SCATTERING_TEXTURE_MU_SIZE / 2);
	} 
	else 
	{
		// Distance to the top atmosphere boundary for the ray (r,mu), and its
		// minimum and maximum values over all mu - obtained for (r,1) and
		// (r,mu_horizon).
		float d = -r_mu + SafeSqrt(discriminant + H * H);
		float d_min = atmosphere._top_radius - r;
		float d_max = rho + H;
		u_mu = 0.5 + 0.5 * GetTextureCoordFromUnitRange((d - d_min) / (d_max - d_min), SCATTERING_TEXTURE_MU_SIZE / 2);
	}

	float d = DistanceToTopAtmosphereBoundary(atmosphere, atmosphere._bottom_radius, mu_s);
	float d_min = atmosphere._top_radius - atmosphere._bottom_radius;
	float d_max = H;
	float a = (d - d_min) / (d_max - d_min);
	float A = -2.0 * atmosphere._mu_s_min * atmosphere._bottom_radius / (d_max - d_min);
	float u_mu_s = GetTextureCoordFromUnitRange(max(1.0 - a / A, 0.0) / (1.0 + a), SCATTERING_TEXTURE_MU_S_SIZE);
	float u_nu = (nu + 1.0) / 2.0;
  
	return vec4(u_nu, u_mu_s, u_mu, u_r);
}

vec3 GetScattering(in AtmosphereParameters atmosphere,in sampler3D scattering_texture, float r, float mu, float mu_s, float nu, bool ray_r_mu_intersects_ground) 
{
	vec4 uvwz = GetScatteringTextureUvwzFromRMuMuSNu(atmosphere, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
	float tex_coord_x = uvwz.x * float(SCATTERING_TEXTURE_NU_SIZE - 1);
	float tex_x = floor(tex_coord_x);
	float lerp = tex_coord_x - tex_x;
	vec3 uvw0 = vec3((tex_x + uvwz.y) / float(SCATTERING_TEXTURE_NU_SIZE), uvwz.z, uvwz.w);
	vec3 uvw1 = vec3((tex_x + 1.0 + uvwz.y) / float(SCATTERING_TEXTURE_NU_SIZE), uvwz.z, uvwz.w);
  
	return mix(texture(scattering_texture, uvw0).rgb, texture(scattering_texture, uvw1).rgb, lerp);//vec3(texture(scattering_texture, uvw0) * (1.0 - lerp) + texture(scattering_texture, uvw1) * lerp);
}

vec3 GetScattering(in AtmosphereParameters atmosphere, in sampler3D single_rayleigh_scattering_texture, in sampler3D single_mie_scattering_texture, in sampler3D multiple_scattering_texture,
					float r, float mu, float mu_s, float nu, bool ray_r_mu_intersects_ground, int scattering_order) 
{
	if (scattering_order == 1) 
	{
		vec3 rayleigh = GetScattering(atmosphere, single_rayleigh_scattering_texture, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
		vec3 mie = GetScattering(atmosphere, single_mie_scattering_texture, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
    
		return rayleigh * RayleighPhaseFunction(nu) + mie * MiePhaseFunction(atmosphere._mie_phase_function_g, nu);
	} 
	else 
	{
		return GetScattering(atmosphere, multiple_scattering_texture, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
	}
}

/* 
-------------------------------------------------- 

MULTIPLE SCATTERING

-----------------------------------------------------
*/

///////// COMPUTATION //////////////////////////////////////////////////

vec3 GetIrradiance(in AtmosphereParameters atmosphere, in sampler2D irradiance_texture, float r, float mu_s);

vec3 ComputeScatteringDensity(	in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, in sampler3D single_rayleigh_scattering_texture,
								in sampler3D single_mie_scattering_texture, in sampler3D multiple_scattering_texture, in sampler2D irradiance_texture,
								float r, float mu, float mu_s, float nu, int scattering_order) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	//assert(mu_s >= -1.0 && mu_s <= 1.0);
	//assert(nu >= -1.0 && nu <= 1.0);
	//assert(scattering_order >= 2);

	// Compute unit direction vectors for the zenith, the view direction omega and
	// and the sun direction omega_s, such that the cosine of the view-zenith
	// angle is mu, the cosine of the sun-zenith angle is mu_s, and the cosine of
	// the view-sun angle is nu. The goal is to simplify computations below.
	vec3 zenith_direction = vec3(0.0, 0.0, 1.0);
	vec3 omega = vec3(sqrt(1.0 - mu * mu), 0.0, mu);
	float sun_dir_x = omega.x == 0.0 ? 0.0 : (nu - mu * mu_s) / omega.x;
	float sun_dir_y = sqrt(max(1.0 - sun_dir_x * sun_dir_x - mu_s * mu_s, 0.0));
	vec3 omega_s = vec3(sun_dir_x, sun_dir_y, mu_s);

	const int SAMPLE_COUNT = 16;
	const float dphi = PI / float(SAMPLE_COUNT);
	const float dtheta = PI / float(SAMPLE_COUNT);
	vec3 rayleigh_mie = vec3(0.0);// * watt_per_cubic_meter_per_sr_per_nm);

	// Nested loops for the integral over all the incident directions omega_i.
	for (int l = 0; l < SAMPLE_COUNT; ++l) 
	{
		float theta = (float(l) + 0.5) * dtheta;
		float cos_theta = cos(theta);
		float sin_theta = sin(theta);
		bool ray_r_theta_intersects_ground = RayIntersectsGround(atmosphere, r, cos_theta);

		// The distance and transmittance to the ground only depend on theta, so we
		// can compute them in the outer loop for efficiency.
		float distance_to_ground = 0.0;// * m;
		vec3 transmittance_to_ground = vec3(0.0);
		vec3 ground_albedo = vec3(0.0);
		
		if (ray_r_theta_intersects_ground) 
		{
			distance_to_ground = DistanceToBottomAtmosphereBoundary(atmosphere, r, cos_theta);
			transmittance_to_ground = GetTransmittance(atmosphere, transmittance_texture, r, cos_theta, distance_to_ground, true /* ray_intersects_ground */);
			ground_albedo = atmosphere._ground_albedo;
		}

		for (int m = 0; m < 2 * SAMPLE_COUNT; ++m) 
		{
			float phi = (float(m) + 0.5) * dphi;
			vec3 omega_i = vec3(cos(phi) * sin_theta, sin(phi) * sin_theta, cos_theta);
			float domega_i = (dtheta /*/ rad*/) * (dphi /*/ rad*/) * sin(theta) /* sr*/;

			// The radiance L_i arriving from direction omega_i after n-1 bounces is
			// the sum of a term given by the precomputed scattering texture for the
			// (n-1)-th order:
			float nu1 = dot(omega_s, omega_i);
			vec3 incident_radiance = GetScattering(atmosphere, single_rayleigh_scattering_texture, single_mie_scattering_texture,  multiple_scattering_texture, r, omega_i.z, mu_s, nu1, ray_r_theta_intersects_ground, scattering_order - 1);

			// and of the contribution from the light paths with n-1 bounces and whose
			// last bounce is on the ground. This contribution is the product of the
			// transmittance to the ground, the ground albedo, the ground BRDF, and
			// the irradiance received on the ground after n-2 bounces.
			vec3 ground_normal = normalize(zenith_direction * r + omega_i * distance_to_ground);
			vec3 ground_irradiance = GetIrradiance(atmosphere, irradiance_texture, atmosphere._bottom_radius, dot(ground_normal, omega_s));
			incident_radiance += transmittance_to_ground * ground_albedo * vec3(1.0 / (PI /*sr*/)) * ground_irradiance;

			// The radiance finally scattered from direction omega_i towards direction
			// -omega is the product of the incident radiance, the scattering
			// coefficient, and the phase function for directions omega and omega_i
			// (all this summed over all particle types, i.e. Rayleigh and Mie).
			float nu2 = dot(omega, omega_i);
			float rayleigh_density = GetProfileDensity(atmosphere._rayleigh_density, r - atmosphere._bottom_radius);
			float mie_density = GetProfileDensity(atmosphere._mie_density, r - atmosphere._bottom_radius);
			
			rayleigh_mie += incident_radiance * (atmosphere._rayleigh_scattering * rayleigh_density * RayleighPhaseFunction(nu2) + atmosphere._mie_scattering * mie_density * MiePhaseFunction(atmosphere._mie_phase_function_g, nu2)) * domega_i;
		}
	}
	return rayleigh_mie;
}

vec3 ComputeMultipleScattering(in AtmosphereParameters atmosphere, in sampler2D transmittance_texture, in sampler3D scattering_density_texture,
								float r, float mu, float mu_s, float nu,
								bool ray_r_mu_intersects_ground) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu >= -1.0 && mu <= 1.0);
	//assert(mu_s >= -1.0 && mu_s <= 1.0);
	//assert(nu >= -1.0 && nu <= 1.0);

	// Number of intervals for the numerical integration.
	const int SAMPLE_COUNT = 50;
  
	// The integration step, i.e. the length of each integration interval.
	float dx = DistanceToNearestAtmosphereBoundary(atmosphere, r, mu, ray_r_mu_intersects_ground) / float(SAMPLE_COUNT);
	
	// Integration loop.
	vec3 rayleigh_mie_sum = vec3 (0.0);// * watt_per_square_meter_per_sr_per_nm);
  
	for (int i = 0; i <= SAMPLE_COUNT; ++i) 
	{
		float d_i = float(i) * dx;

		// The r, mu and mu_s parameters at the current integration point (see the
		// single scattering section for a detailed explanation).
		float r_i = ClampRadius(atmosphere, sqrt(d_i * d_i + 2.0 * r * mu * d_i + r * r));
		float mu_i = ClampCosine((r * mu + d_i) / r_i);
		float mu_s_i = ClampCosine((r * mu_s + d_i * nu) / r_i);

		// The Rayleigh and Mie multiple scattering at the current sample point.
		vec3 rayleigh_mie_i = 	GetScattering(atmosphere, scattering_density_texture, r_i, mu_i, mu_s_i, nu, ray_r_mu_intersects_ground) *
								GetTransmittance(atmosphere, transmittance_texture, r, mu, d_i, ray_r_mu_intersects_ground) * dx;
		
		// Sample weight (from the trapezoidal rule).
		float weight_i = (i == 0 || i == SAMPLE_COUNT) ? 0.5 : 1.0;
		rayleigh_mie_sum += rayleigh_mie_i * weight_i;
	}	
	return rayleigh_mie_sum;
}

///////// PRE-COMPUTATION //////////////////////////////////////////////////

vec3 ComputeScatteringDensityTexture(
										in AtmosphereParameters atmosphere,
										in sampler2D transmittance_texture,
										in sampler3D single_rayleigh_scattering_texture,
										in sampler3D single_mie_scattering_texture,
										in sampler3D multiple_scattering_texture,
										in sampler2D irradiance_texture,
										in vec3 frag_coord, int scattering_order) 
{
	float r;
	float mu;
	float mu_s;
	float nu;
	bool ray_r_mu_intersects_ground;
	
	GetRMuMuSNuFromScatteringTextureFragCoord(atmosphere, frag_coord, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
  
	return ComputeScatteringDensity(atmosphere, transmittance_texture, single_rayleigh_scattering_texture, single_mie_scattering_texture,
									multiple_scattering_texture, irradiance_texture, r, mu, mu_s, nu, scattering_order);
}

vec3 ComputeMultipleScatteringTexture(
										in AtmosphereParameters atmosphere,
										in sampler2D transmittance_texture,
										in sampler3D scattering_density_texture,
										in vec3 frag_coord, out float nu) 
{
	float r;
	float mu;
	float mu_s;
	bool ray_r_mu_intersects_ground;
	
	GetRMuMuSNuFromScatteringTextureFragCoord(atmosphere, frag_coord, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
  
	return ComputeMultipleScattering(atmosphere, transmittance_texture, scattering_density_texture, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
}

/* 
-------------------------------------------------- 

GROUND IRRANDIANCE

-----------------------------------------------------
*/

///////////// COMPUTATION ////////////////////////////////////////////////////////////////
vec3 ComputeDirectIrradiance(
								in AtmosphereParameters atmosphere,
								in sampler2D transmittance_texture,
								float r, float mu_s) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu_s >= -1.0 && mu_s <= 1.0);

	float alpha_s = atmosphere._sun_angular_radius /* / rad*/;
	
	// Approximate average of the cosine factor mu_s over the visible fraction of
	// the Sun disc.
	float average_cosine_factor = mu_s < -alpha_s ? 0.0 : (mu_s > alpha_s ? mu_s : (mu_s + alpha_s) * (mu_s + alpha_s) / (4.0 * alpha_s));

	return atmosphere._solar_irradiance * GetTransmittanceToTopAtmosphereBoundary(atmosphere, transmittance_texture, r, mu_s) * average_cosine_factor;
}

vec3 ComputeIndirectIrradiance(
								in AtmosphereParameters atmosphere,
								in sampler3D single_rayleigh_scattering_texture,
								in sampler3D single_mie_scattering_texture,
								in sampler3D multiple_scattering_texture,
								float r, float mu_s, int scattering_order) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu_s >= -1.0 && mu_s <= 1.0);
	//assert(scattering_order >= 1);

	const int SAMPLE_COUNT = 32;
	const float dphi = PI / float(SAMPLE_COUNT);
	const float dtheta = PI / float(SAMPLE_COUNT);

	vec3 result = vec3(0.0);// * watt_per_square_meter_per_nm);
	vec3 omega_s = vec3(sqrt(1.0 - mu_s * mu_s), 0.0, mu_s);
  
	for (int j = 0; j < SAMPLE_COUNT / 2; ++j) 
	{
		float theta = (float(j) + 0.5) * dtheta;
    
		for (int i = 0; i < 2 * SAMPLE_COUNT; ++i) 
		{
			float phi = (float(i) + 0.5) * dphi;
			vec3 omega = vec3(cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta));
      
			float domega = (dtheta /* / rad*/) * (dphi /* / rad*/) * sin(theta) /*sr*/;
			float nu = dot(omega, omega_s);
      
			result += GetScattering(atmosphere, single_rayleigh_scattering_texture, single_mie_scattering_texture, multiple_scattering_texture, 
									r, omega.z, mu_s, nu, false /* ray_r_theta_intersects_ground */,
									scattering_order) *
									omega.z * domega;
		}
	}
	return result;
}

///////////// PRE-COMPUTATION ////////////////////////////////////////////////////////////////

void GetRMuSFromIrradianceTextureUv(in AtmosphereParameters atmosphere, in vec2 uv, out float r, out float mu_s) 
{
	//assert(uv.x >= 0.0 && uv.x <= 1.0);
	//assert(uv.y >= 0.0 && uv.y <= 1.0);
	float x_mu_s = GetUnitRangeFromTextureCoord(uv.x, IRRADIANCE_TEXTURE_WIDTH);
	float x_r = GetUnitRangeFromTextureCoord(uv.y, IRRADIANCE_TEXTURE_HEIGHT);
  
	r = atmosphere._bottom_radius + x_r * (atmosphere._top_radius - atmosphere._bottom_radius);
	mu_s = ClampCosine(2.0 * x_mu_s - 1.0);
}

const vec2 IRRADIANCE_TEXTURE_SIZE = vec2(IRRADIANCE_TEXTURE_WIDTH, IRRADIANCE_TEXTURE_HEIGHT);

vec3 ComputeDirectIrradianceTexture(
										in AtmosphereParameters atmosphere,
										in sampler2D transmittance_texture,
										in vec2 frag_coord) 
{
	float r;
	float mu_s;
	
	GetRMuSFromIrradianceTextureUv(atmosphere, frag_coord / IRRADIANCE_TEXTURE_SIZE, r, mu_s);
	
	return ComputeDirectIrradiance(atmosphere, transmittance_texture, r, mu_s);
}

vec3 ComputeIndirectIrradianceTexture(
														in AtmosphereParameters atmosphere,
														in sampler3D single_rayleigh_scattering_texture,
														in sampler3D single_mie_scattering_texture,
														in sampler3D multiple_scattering_texture,
														in vec2 frag_coord, int scattering_order) 
{
	float r;
	float mu_s;
  
	GetRMuSFromIrradianceTextureUv(atmosphere, frag_coord / IRRADIANCE_TEXTURE_SIZE, r, mu_s);
  
	return ComputeIndirectIrradiance(atmosphere, single_rayleigh_scattering_texture, single_mie_scattering_texture, multiple_scattering_texture, r, mu_s, scattering_order);
}

///////////// LOOK-UP ////////////////////////////////////////////////////////////////

vec2 GetIrradianceTextureUvFromRMuS(in AtmosphereParameters atmosphere, float r, float mu_s) 
{
	//assert(r >= atmosphere.bottom_radius && r <= atmosphere.top_radius);
	//assert(mu_s >= -1.0 && mu_s <= 1.0);
	float x_r = (r - atmosphere._bottom_radius) / (atmosphere._top_radius - atmosphere._bottom_radius);
	float x_mu_s = mu_s * 0.5 + 0.5;
	
	return vec2(GetTextureCoordFromUnitRange(x_mu_s, IRRADIANCE_TEXTURE_WIDTH), GetTextureCoordFromUnitRange(x_r, IRRADIANCE_TEXTURE_HEIGHT));
}

vec3 GetIrradiance(
					in AtmosphereParameters atmosphere,
					in sampler2D irradiance_texture,
					float r, float mu_s) 
{
	vec2 uv = GetIrradianceTextureUvFromRMuS(atmosphere, r, mu_s);

	return vec3(texture2D(irradiance_texture, uv).rgb);
}


/*
-----------------------------------------------------

RENDERING

-----------------------------------------------------
*/

#ifdef COMBINED_SCATTERING_TEXTURES
	vec3 GetExtrapolatedSingleMieScattering(in AtmosphereParameters atmosphere, in vec4 scattering) 
	{
		if (scattering.r == 0.0) 
			return vec3(0.0);
			
		return scattering.rgb * scattering.a / scattering.r * (atmosphere._rayleigh_scattering.r / atmosphere._mie_scattering.r) * (atmosphere._mie_scattering / atmosphere._rayleigh_scattering);
	}
#endif

vec3 GetCombinedScattering(
							in AtmosphereParameters atmosphere,
							in sampler3D scattering_texture,
							in sampler3D single_mie_scattering_texture,
							float r, float mu, float mu_s, float nu,
							bool ray_r_mu_intersects_ground,
							out vec3 single_mie_scattering) 
{
	vec4 uvwz = GetScatteringTextureUvwzFromRMuMuSNu(atmosphere, r, mu, mu_s, nu, ray_r_mu_intersects_ground);
  
	float tex_coord_x = uvwz.x * float(SCATTERING_TEXTURE_NU_SIZE - 1);
	float tex_x = floor(tex_coord_x);
	float lerp = tex_coord_x - tex_x;
	vec3 uvw0 = vec3((tex_x + uvwz.y) / float(SCATTERING_TEXTURE_NU_SIZE), uvwz.z, uvwz.w);
	vec3 uvw1 = vec3((tex_x + 1.0 + uvwz.y) / float(SCATTERING_TEXTURE_NU_SIZE), uvwz.z, uvwz.w);

#ifdef COMBINED_SCATTERING_TEXTURES
	vec4 combined_scattering = texture3D(scattering_texture, uvw0) * (1.0 - lerp) + texture3D(scattering_texture, uvw1) * lerp;
	vec3 scattering = vec3(combined_scattering);
	single_mie_scattering = GetExtrapolatedSingleMieScattering(atmosphere, combined_scattering);
#else
	vec3 scattering = vec3(texture3D(scattering_texture, uvw0) * (1.0 - lerp) + texture3D(scattering_texture, uvw1) * lerp);
	single_mie_scattering = vec3(texture3D(single_mie_scattering_texture, uvw0) * (1.0 - lerp) + texture3D(single_mie_scattering_texture, uvw1) * lerp);
#endif
  
	return scattering;
}

/////////////////////////////////////////////////////////////
///// SKY ///////////////////////////////////////////////////
/////////////////////////////////////////////////////////////
vec3 GetSkyRadiance(
					in AtmosphereParameters atmosphere,
					in sampler2D transmittance_texture,
					in sampler3D scattering_texture,
					in sampler3D single_mie_scattering_texture,
					vec3 camera, in vec3 view_ray, float shadow_length,
					in vec3 sun_direction, out vec3 transmittance) 
{
	// Compute the distance to the top atmosphere boundary along the view ray,
	// assuming the viewer is in space (or NaN if the view ray does not intersect
	// the atmosphere).
	float r = length(camera);
	float rmu = dot(camera, view_ray);
	float distance_to_top_atmosphere_boundary = -rmu - sqrt(rmu * rmu - r * r + atmosphere._top_radius * atmosphere._top_radius);
  
	// If the viewer is in space and the view ray intersects the atmosphere, move
	// the viewer to the top atmosphere boundary (along the view ray):
	if (distance_to_top_atmosphere_boundary > 0.0) // * m ) 
	{
		camera = camera + view_ray * distance_to_top_atmosphere_boundary;
		r = atmosphere._top_radius;
		rmu += distance_to_top_atmosphere_boundary;
	} 
	else if (r > atmosphere._top_radius) 
	{
		// If the view ray does not intersect the atmosphere, simply return 0.
		transmittance = vec3(1.0);
		return vec3(0.0 /* watt_per_square_meter_per_sr_per_nm*/);
	}
  
	// Compute the r, mu, mu_s and nu parameters needed for the texture lookups.
	float mu = rmu / r;
	float mu_s = dot(camera, sun_direction) / r;
	float nu = dot(view_ray, sun_direction);
	bool ray_r_mu_intersects_ground = RayIntersectsGround(atmosphere, r, mu);

	transmittance = ray_r_mu_intersects_ground ? vec3(0.0) : GetTransmittanceToTopAtmosphereBoundary(atmosphere, transmittance_texture, r, mu);
	
	vec3 single_mie_scattering;
	vec3 scattering;
	
	if (shadow_length == 0.0)// * m) 
	{
		scattering = GetCombinedScattering(atmosphere, scattering_texture, single_mie_scattering_texture, r, mu, mu_s, nu, ray_r_mu_intersects_ground, single_mie_scattering);
	} 
	else 
	{
		// Case of light shafts (shadow_length is the total length noted l in our
		// paper): we omit the scattering between the camera and the point at
		// distance l, by implementing Eq. (18) of the paper (shadow_transmittance
		// is the T(x,x_s) term, scattering is the S|x_s=x+lv term).
		float d = shadow_length;
		float r_p = ClampRadius(atmosphere, sqrt(d * d + 2.0 * r * mu * d + r * r));
		float mu_p = (r * mu + d) / r_p;
		float mu_s_p = (r * mu_s + d * nu) / r_p;

		scattering = GetCombinedScattering(atmosphere, scattering_texture, single_mie_scattering_texture, r_p, mu_p, mu_s_p, nu, ray_r_mu_intersects_ground, single_mie_scattering);
		vec3 shadow_transmittance = GetTransmittance(atmosphere, transmittance_texture, r, mu, shadow_length, ray_r_mu_intersects_ground);
		
		scattering = scattering * shadow_transmittance;
		single_mie_scattering = single_mie_scattering * shadow_transmittance;
  }
  
  return scattering * RayleighPhaseFunction(nu) + single_mie_scattering * MiePhaseFunction(atmosphere._mie_phase_function_g, nu);
}

/////////////////////////////////////////////////////////////
// AERIAL PERSPECTIVE ///////////////////////////////////////
/////////////////////////////////////////////////////////////
vec3 GetSkyRadianceToPoint(
							in AtmosphereParameters atmosphere,
							in sampler2D transmittance_texture,
							in sampler3D scattering_texture,
							in sampler3D single_mie_scattering_texture,
							vec3 camera, in vec3 point, float shadow_length,
							in vec3 sun_direction, out vec3 transmittance) 
{
	// Compute the distance to the top atmosphere boundary along the view ray,
	// assuming the viewer is in space (or NaN if the view ray does not intersect
	// the atmosphere).
	vec3 view_ray = normalize(point - camera);
	float r = length(camera);
	float rmu = dot(camera, view_ray);
	float distance_to_top_atmosphere_boundary = -rmu - sqrt(rmu * rmu - r * r + atmosphere._top_radius * atmosphere._top_radius);
	
	// If the viewer is in space and the view ray intersects the atmosphere, move
	// the viewer to the top atmosphere boundary (along the view ray):
	if (distance_to_top_atmosphere_boundary > 0.0)// * m) 
	{
		camera = camera + view_ray * distance_to_top_atmosphere_boundary;
		r = atmosphere._top_radius;
		rmu += distance_to_top_atmosphere_boundary;
	}

	// Compute the r, mu, mu_s and nu parameters for the first texture lookup.
	float mu = rmu / r;
	float mu_s = dot(camera, sun_direction) / r;
	float nu = dot(view_ray, sun_direction);
	float d = length(point - camera);
	bool ray_r_mu_intersects_ground = RayIntersectsGround(atmosphere, r, mu);

	transmittance = GetTransmittance(atmosphere, transmittance_texture, r, mu, d, ray_r_mu_intersects_ground);

	vec3 single_mie_scattering;
	vec3 scattering = GetCombinedScattering(atmosphere, scattering_texture, single_mie_scattering_texture, r, mu, mu_s, nu, ray_r_mu_intersects_ground, single_mie_scattering);

	// Compute the r, mu, mu_s and nu parameters for the second texture lookup.
	// If shadow_length is not 0 (case of light shafts), we want to ignore the
	// scattering along the last shadow_length meters of the view ray, which we
	// do by subtracting shadow_length from d (this way scattering_p is equal to
	// the S|x_s=x_0-lv term in Eq. (17) of our paper).
	d = max(d - shadow_length, 0.0);// * m);
	float r_p = ClampRadius(atmosphere, sqrt(d * d + 2.0 * r * mu * d + r * r));
	float mu_p = (r * mu + d) / r_p;
	float mu_s_p = (r * mu_s + d * nu) / r_p;

	vec3 single_mie_scattering_p;
	vec3 scattering_p = GetCombinedScattering(atmosphere, scattering_texture, single_mie_scattering_texture, r_p, mu_p, mu_s_p, nu, ray_r_mu_intersects_ground, single_mie_scattering_p);

	// Combine the lookup results to get the scattering between camera and point.
	vec3 shadow_transmittance = transmittance;
	if (shadow_length > 0.0)// * m) 
	{
		// This is the T(x,x_s) term in Eq. (17) of our paper, for light shafts.
		shadow_transmittance = GetTransmittance(atmosphere, transmittance_texture, r, mu, d, ray_r_mu_intersects_ground);
	}
  
	scattering = scattering - shadow_transmittance * scattering_p;
	single_mie_scattering = single_mie_scattering - shadow_transmittance * single_mie_scattering_p;

	#ifdef COMBINED_SCATTERING_TEXTURES
		single_mie_scattering = GetExtrapolatedSingleMieScattering(atmosphere, vec4(scattering, single_mie_scattering.r));
	#endif

	// Hack to avoid rendering artifacts when the sun is below the horizon.
	single_mie_scattering = single_mie_scattering * smoothstep(float(0.0), float(0.01), mu_s);

	return scattering * RayleighPhaseFunction(nu) + single_mie_scattering * MiePhaseFunction(atmosphere._mie_phase_function_g, nu);
}

/////////////////////////////////////////////////////////////
// GROUND ///////////////////////////////////////////////////
/////////////////////////////////////////////////////////////
vec3 GetSunAndSkyIrradiance( 
							in AtmosphereParameters atmosphere, 
							in sampler2D transmittance_texture, 
							in sampler2D irradiance_texture, 
							in vec3 point, in vec3 normal, in vec3 sun_direction,
							out vec3 sky_irradiance) 
{
	float r = length(point);
	float mu_s = dot(point, sun_direction) / r;

	// Indirect irradiance (approximated if the surface is not horizontal).
	sky_irradiance = GetIrradiance(atmosphere, irradiance_texture, r, mu_s) * vec3((1.0 + dot(normal, point) / r) * 0.5);

	// Direct irradiance.
	return atmosphere._solar_irradiance * GetTransmittanceToSun(atmosphere, transmittance_texture, r, mu_s) * max(dot(normal, sun_direction), 0.0);
}