#ifndef ATMOSPHERE_FUNCTIONS
#define ATMOSPHERE_FUNCTIONS
#include "Packages/com.unity.render-pipelines.universal/Shaders/LitInput.hlsl"
#include "AtmosphereConstants.hlsl"



float GetHazeFactorAtAltitude(float alt_km)
{
	if (alt_km < mieBaseHeightKm)
		alt_km = mieBaseHeightKm;
	float val = exp((mieBaseHeightKm - alt_km) / mieScaleHeightKm);
	return val;
}

float GetMieFactorAtAltitude(float alt_km)
{
	if (alt_km<0)
		alt_km = 0;
	float val = exp(-alt_km / RAYLEIGH_SCALE_HEIGHT_KM);
	return val;
}

float GetShortestDistanceToAltitude(float sine_elevation,float start_h_km,float finish_h_km)
{
	float RH		=planetRadiusKm+finish_h_km;
	float Rh		=planetRadiusKm+start_h_km;
	float cosine	=-sine_elevation;
	float b			=-2.0*Rh*cosine;
	float c			=Rh*Rh-RH*RH;
	float b24c		=b*b-4*c;
	if(b24c<0)
		return -1.0;
	float dist;
	float s=sqrt(b24c);
	//if(b+s>=0.0)
		dist=0.5*(-b+s);
	//else
	//	dist=0.5*(-b-s);
	return dist;
}

float GetDistanceToSpace(float sine_elevation,float h_km)
{
	return GetShortestDistanceToAltitude(sine_elevation,h_km, atmosphereThicknessKm);
}
//
//float3 Convert2DCoordTo3DCoord(float2 uv)
//{
//	float3 texc;
//	texc.x = uv
//}

float TexcToAltKm(float texc, float minOutputAltKm, float maxOutputAltKm)
{
	return minOutputAltKm + texc*texc*(maxOutputAltKm - minOutputAltKm);
}

float GetHazeOpticalLength(float sine_elevation,float h_km)
{
	float R=planetRadiusKm;
	float Rh=R+h_km;
	float RH=R+mieBaseHeightKm;
	float c=sqrt(1.0-sine_elevation*sine_elevation);
	float u=RH*RH-Rh*Rh*c*c;
	float U=R*R-Rh*Rh*c*c;
	// If the ray passes through the earth, infinite optical length.
	//if(sine_elevation<0&&U>0.0)
	//	return 1000000.0;
	float haze_opt_len=0.0;
	// If we have a solution, there exists a path through the constant haze area.
	if(sine_elevation<0&&u>0)
		haze_opt_len=2.0*sqrt(u);
	// Next, the inward path, if that exists.
	float Rmin=Rh*c;
	if(sine_elevation>0.0)
		Rmin=Rh;
	// But if the inward path goes into the constant area, include only the part outside that area.
	if(sine_elevation<0.0&&u>0.0)
		Rmin=RH;
	float h1=Rh-RH;
	float h2=Rmin-RH;
	float s=sine_elevation;
	float n=mieScaleHeightKm;
	// s<0, and h2<h1
	if(s<0.0)
		haze_opt_len+=n/s*(saturate(exp(-h1/n))-saturate(exp(-h2/n)));
	// Now the outward path, in this case h2 -> infinity
	// and elevation is reversed.
	if(s<0.0)
		s*=-1.0;
	if(s<0.01)
		s=0.01;
	haze_opt_len+=n/s*(exp(-abs(h2)/n));
	return haze_opt_len;
}

float4 GetSunlightFactor(float alt_km, float3 DirectionToLight)
{
	float sine = clamp(DirectionToLight.y, -1.0, 1.0);
	/*float2 table_texc;
	table_texc.y = (0.5 - 0.5*sine)*NUM_ELEVATIONS;
	table_texc.x = sqrt(max(0, alt_km - minOutputAltKm)/(maxOutputAltKm - minOutputAltKm))*NUM_ALTITUDES;
	table_texc += float2(0.5, 0.5);
	table_texc = float2(table_texc.x / (NUM_ALTITUDES*NUM_DISTANCES), table_texc.y / NUM_ELEVATIONS);
	table_texc.x = table_texc.x + (NUM_DISTANCES-1.0) / NUM_DISTANCES;*/

	//float3 loss = SAMPLE_TEXTURE2D(loss_texture, sampler_loss_texture, table_texc);
	//float vis = GetVisibleFactor(alt_km, sine);
	float2 table_texc = float2(DENSITY_TABLE_SIZE * (0.5 + 0.5 * sine), DENSITY_TABLE_SIZE * (alt_km / maxDensityAltKm));

	table_texc += float2(0.5, 0.5);
	table_texc = float2(table_texc.x / DENSITY_TABLE_SIZE, table_texc.y / DENSITY_TABLE_SIZE);

	float4 lookup = SAMPLE_TEXTURE2D(optical_depth_texture, sampler_optical_depth_texture, table_texc);
	float illuminated_length = lookup.x;
	float vis = lookup.y;
	float ozone_length = lookup.w;
	float haze_opt_len = GetHazeOpticalLength(sine, alt_km);

	float4 factor = float4(vis, vis, vis, vis);
	//factor.rgb *= loss;
	// test only
	//factor.rgb *= exp(-rayleigh * illuminated_length - hazeMie * haze_opt_len - ozone * ozone_length);
	
	return factor;
}

float CalcRayleighBeta(float cos0)
{
	return (1.0 + cos0*cos0); // Factor of (0.0596831) is applied in texture generation now.
}

float HenyeyGreenstein(float g, float cos0)
{
	float g2 = g*g;
	float u = abs(1.0 + g2 - 2.0*g*cos0);
	return (1.0 - g2)*pow(u, -1.5) / (4.0*PI);
}

half3 PrecalculatedInscatterFunction(half4 inscatter_factor, float BetaRayleigh, float BetaMie, float3 mieRayleighRatio)
{
	float3 mr = mieRayleighRatio.xyz;//*inscatter_factor.a
	float3 BetaTotal = (BetaRayleigh + BetaMie*mr) / (float3(1.0, 1.0, 1.0) + mr);
	half3 colour = BetaTotal*inscatter_factor.rgb;
	return colour;
}

half3 InscatterFunction(half4 inscatter_factor, float hazeEccentricity, float cos0, float3 mieRayleighRatio)
{
	float BetaRayleigh = CalcRayleighBeta(cos0);
	float BetaMie = HenyeyGreenstein(hazeEccentricity, cos0);		// Mie's phase function
	return PrecalculatedInscatterFunction(inscatter_factor, BetaRayleigh, BetaMie, mieRayleighRatio);
}

float3 GetSkylight(float alt_km, TEXTURE2D_PARAM(insc_texture, sampler_insc_texture))
{
	// The inscatter factor, at this altitude looking straight up, is given by:
	float2 inscTexc;
	inscTexc.x = sqrt(alt_km / maxOutputAltKm) * NUM_ALTITUDES / (NUM_ALTITUDES*NUM_DISTANCES) + 1.0 / NUM_DISTANCES;
	inscTexc.y = 0;
	float4 insc = SAMPLE_TEXTURE2D(insc_texture, sampler_insc_texture, inscTexc);
	float3 skylight = InscatterFunction(insc, hazeEccentricity, 0.0, mieRayleighRatio);
	return skylight;
}

float4 Skyl(TEXTURE2D_PARAM(insc_texture, sampler_insc_texture)
		,TEXTURE2D_PARAM(density_texture, sampler_density_texture)
		,TEXTURE2D_PARAM(blackbody_texture, sampler_blackbody_texture)
		,float3 previous_loss
		,float4 previous_skyl
		,float spaceDistKm
		,float viewAltKm
		,float dist_km
		,float prevDist_km
		,float sin_e
		,float cos_e)
{
	float maxd			=min(spaceDistKm,dist_km);
	float mind			=min(spaceDistKm,prevDist_km);
	float dist			=0.5*(mind+maxd);
	float stepLengthKm	=max(0.0,maxd-mind);
	float y				=planetRadiusKm+viewAltKm+dist*sin_e;
	float x				=dist*cos_e;
	float r				=sqrt(x*x+y*y);
	float alt_km		=r-planetRadiusKm;
	// lookups is: dens_factor,ozone_factor,haze_factor;
	float dens_texc		=(alt_km/maxDensityAltKm*(DENSITY_TABLE_SIZE -1.0)+0.5)/ DENSITY_TABLE_SIZE;
	float4 lookups		=SAMPLE_TEXTURE2D(density_texture, sampler_density_texture, float2(dens_texc,0.5));
	float dens_factor	=lookups.x;
	float ozone_factor	=lookups.y;
	float haze_factor	=GetHazeFactorAtAltitude(alt_km);
	float4 light		=float4(starlight+GetSkylight(alt_km, TEXTURE2D_ARGS(insc_texture, sampler_insc_texture)),0.0);
	float4 skyl			=light;
	float3 extinction	=dens_factor*rayleigh+haze_factor*hazeMie;
	float3 total_ext	=extinction+ozone*ozone_factor;
	float3 loss			=exp(-extinction*stepLengthKm);
	skyl.rgb			*=float3(1.0,1.0,1.0)-loss;
	float mie_factor	=exp(-skyl.w*stepLengthKm*haze_factor*hazeMie.x);
	skyl.w				=saturate((1.0-mie_factor)/(1.0-total_ext.x+0.0001));
#if 1//def BLACKBODY
	float dens_dist		=dens_factor*stepLengthKm;
	float emis_ext		=exp(-emissivity*dens_dist);
	float3 bb;
	//float T				=seaLevelTemperatureK*lookups.w;
	bb.xyz = float3(54.48783, 117.76968, 293.04865);//Blackbody(blackbody_texture,T);

	skyl				*=emis_ext;
	bb					*=1.0-emis_ext;
	skyl.rgb			+=bb;
#endif
	skyl.rgb			*=previous_loss.rgb;
	skyl.rgb			+=previous_skyl.rgb;
	float lossw			=1.0;
	skyl.w				=lossw*(1.0-previous_skyl.w)*skyl.w+previous_skyl.w;
	return skyl;
}





#endif // ATMOSPHERE_FUNCTIONS
