in vec2 uv;
in vec2 uvn;

layout (location = 0) out vec4 C;

uniform int prev_frame_tex;
uniform int pt_frame;


vec3 car_to_pol(float th, float phi){
	return vec3(
	sin(th)*sin(phi),
	sin(th)*cos(phi),
	cos(th)
	);
}


vec3 sample_hemi(){
	// pdf(w)       = 1./pi
	// pdf(the,phi) = sin(the)/pi

	// marg dens theta     - sin(the)
	// cond dens phi/theta - (sin(the)/pi) / sin(the) = 1./pi

	// --- integrate & inverse
	// x_theta = cos(the) = acos(the)
	// x_phi   = phi/2pi  = phi*2pi

	vec2 x = hash_v2();
	x.x = acos(x.x);
	x.y *= tau;
	return car_to_pol(x.x,x.y);
}
vec3 sample_sph(){
	vec3 s = sample_hemi();
	if(hash_f() < 0.5)
	s = -s;
	return s;
}
vec4 qsqr(vec4 a){return vec4(a.x*a.x-a.y*a.y-a.z*a.z-a.w*a.w, 2.0*a.x*a.y, 2.0*a.x*a.z, 2.0*a.x*a.w);}
float sdfMandelbulb(vec3 _z, out vec3 diffuseColor) {
	// return sdMandelbulb(_z).w;
	vec3 c = vec3(-1., .37, .0);
	vec3 z = _z;
//	z *= rotationMatrix(vec2f(-0.2, .01).xxy);
	z.xy *= rot(0.1);
	z.yz *= rot(0.1);
	vec3 orbitTrap = vec3(1);
	float r = length(z);
	float dr = 1.;
	float xxyy;
	for (int i = 0; i < 35; i++) {
		xxyy = z.x*z.x+z.y*z.y;
		r = sqrt(xxyy+z.z*z.z);
		orbitTrap = min(orbitTrap, abs(z*3.));
		if (r>5.) {break;}
		dr = r*2.*dr + 1.;
		z = vec3(vec2(z.x*z.x-z.y*z.y, 2.*z.x*z.y)*(1.-z.z*z.z/xxyy), -2.*z.z*sqrt(xxyy)) + c;
	}
	diffuseColor = orbitTrap*11115.0;
	float dd = 0.1*log(r)*r/dr - 0.002;
	// dd = min(dd, length(_z) - 0.0001);
	return dd;
}

vec4 distanceestimator(vec3 pos){
	if (true){
		vec3 c = vec3(1);
		float d = sdfMandelbulb(pos, c);
		return vec4(c,d);
	}
	vec4 z = vec4(pos, 0.0);
	float md2 = 1.0;
	float mz2 = dot(z, z);
	vec4 orbitTrap = vec4(1.0);
	for(int i = 0; i < 8; i++){
		md2 *= 4.0*mz2;
		z = qsqr(z)+vec4(-0.5, 0.35, 0.5, 0.0);
		//z = qsqr(z)+vec4(-0.5, 0.5, 0.25, 0.0);
		orbitTrap = min(abs(z), orbitTrap);
		mz2 = dot(z,z);
		if(mz2 > 4.0) break;}
	float sdf = 0.25*sqrt(mz2/md2)*log(mz2);
	return vec4(orbitTrap.rgb, sdf);
}

#define max_volumedensity 70.0

vec3 rgb2hsv_2(vec3 c){
	vec4 K=vec4(0.,-1./3.,2./3.,-1.),
	p=mix(vec4(c.bg ,K.wz),vec4(c.gb,K.xy ),step(c.b,c.g)),
	q=mix(vec4(p.xyw,c.r ),vec4(c.r ,p.yzx),step(p.x,c.r));
	float d=q.x-min(q.w,q.y),
	e=1e-10;
	return vec3(abs(q.z+(q.w-q.y)/(6.*d+e)),d/(q.x+e),q.x);
}

// 3D Volumetric Density Function
vec4 densityfunction(vec3 pathposition){
	float density = 0.0;
	vec4 distanceestimation = distanceestimator(pathposition);
	if(distanceestimation.w < 0.0){
		density = max_volumedensity;
		distanceestimation.rgb /= max(1.- pow(rgb2hsv_2(distanceestimation.rgb).y,2),0.01);
		distanceestimation.xyz = 1.-exp(-distanceestimation.xyz*44.);
		distanceestimation.rgb *= 0.4;
//		distanceestimation.rgb *= distanceestimation.rgb;
//		distanceestimation.rgb *= 0.5;
	}
	return vec4(distanceestimation.rgb, density);
}



vec4 get_vol(vec3 p){
	if(false){
		float d = length(p) - 0.5;
		return vec4(float(d < 0.)*4.0,0.1,0.1,0.1);
	} else {
		vec4 d = densityfunction(p);
		return d.wxyz;
	}
//	return vec4(*4.,);
}

float blue_noise(vec2 _u, uint s) {
	uint f = F;
//	f = 125;
	float v = 0.;
	for (int k=0; k<9; k++){
		uvec2 offs = uvec2(k%3-1,k/3-1);
		uvec2 u = uvec2(_u) + offs;
		v += hash_f_s(hashi(u.x + f*52u) + hashi(u.y + f*2) + hashi(f) + s*12);
	}
	//return       1.125*hash(U)- v/8.  + .5; // some overbound, closer contrast
	return .9 *( 1.125*hash_f_s(hashi(uint(_u.x) + f*52u) + hashi(uint(_u.y) + f*2) + hashi(f)+ s*12)- v/8.) + .5; //
	//return .75*( 1.125*hash(U)- v/8.) + .5; // very slight overbound
	//return .65*( 1.125*hash(U)- v/8.) + .5; // dimmed, but histo just fit without clamp. flat up to .5 +- .23
}


float henyey_brdf(float cos_theta, float aniso){
	float denom = 1 + aniso * aniso + 2 * aniso * cos_theta;
	return (1./(4.*pi)) * (1.0 - aniso * aniso) / (denom * sqrt(denom));
}
//float henyey_brdf(float cos_theta, float aniso){
//	return 0.25 * (1./pi) * (1.0 - aniso*aniso) / pow(1.0 + aniso * (aniso - 2.0 * cos_theta), 1.5);
//}

vec3 sample_henyey(float g) {
	float cosTheta;
	if (abs(g) < 1e-3)
		cosTheta = 1. - 2. * hash_f();
	else {
		float sqrTerm = (1. - g * g) /
		(1. - g + 2. * g * hash_f());
		cosTheta = (1. + g * g - sqrTerm * sqrTerm) / (2. * g);
	}

	return car_to_pol(acos(clamp(cosTheta, -1., 1.0)) , hash_f()*tau);
//	return henyey_brdf(cosTheta, g);
}

float henyey_greenstein_pdf(float cos_theta, float g){
	float nom = 1.-g*g;
	float den = 1.+g*g + 2.*g*cos_theta;
	den *= sqrt(den);
	return nom/max(den,0.001) * (1./(4.*pi));
}

float henyeyGreensteinIntegral(float cosTheta, float a) {
	if(a == 0.0) return cosTheta;
	return (1.0 - a * a) / (a * sqrt(1.0 + a * (a - 2.0 * cosTheta)));
}

vec3 volumeBRDFSample(float g) {
	float cosTheta = 1.0 + henyeyGreensteinIntegral(-1.0, g) - henyeyGreensteinIntegral(hash_f() * 2.0 - 1.0, g);
	float r = sqrt(max(1.0 - cosTheta * cosTheta, 0.0));
	float a = hash_f() * tau;
	return vec3(r * cos(a), r * sin(a), cosTheta);
}

float volumeBRDFPDF(float cos_theta, float g){
	return 0.25 * (1/pi) * (1.0 - g * g) / pow(1.0 + g * (g - 2.0 * cos_theta), 1.5);
}

float rayleigh_brdf(float cos_theta){
	return 3.0 * (1.0 + cos_theta * cos_theta) / (pi * 16.0);
}

vec3 rayleigh_sample(){
	vec2 r = vec2(tau * hash_f(),
	acos(2.0 * hash_f() - 1.0));
	vec2 c = cos(r), s = sin(r);
	return vec3(s.y * s.x, s.y * c.x, c.y);
}


float anis = 0.8;

float standardObserver1931_w_min = 380.0f;
float standardObserver1931_w_max = 780.0f;

float trace_vol(inout vec3 p, inout vec4 vol_samp, vec3 rd, float wv){
//	float rayleigh_fac = wv*wv*wv*wv * 2.0e-11;
	float rayleigh_fac = wv*wv*wv*wv;
	rayleigh_fac = 1/(rayleigh_fac);
	rayleigh_fac /= 2.0e-13;
	rayleigh_fac *= 0.005;
//	rayleigh_fac *= pow(10,13);
//	rayleigh_fac *= 2.0;
//	rayleigh_fac = 20e-11/rayleigh_fac;
//	rayleigh_fac = 20e-11/rayleigh_fac;
//	rayleigh_fac = 1-rayleigh_fac;

//	float rayleigh_fac = standardObserver1931_w_min/(wv*wv*wv*wv);
//	rayleigh_fac = 1.-rayleigh_fac;

//	float mv_dens = max_volumedensity / rayleigh_fac;
	float mv_dens = max_volumedensity * rayleigh_fac;
	for(int i = 0; i < 450; i++){
		float t = -log(1.-hash_f()) / mv_dens;

		p = p + rd * t;
		vol_samp = get_vol(p);
		float dens = vol_samp.x;
//		dens = max_volumedensity;
//		dens *= 2511111111./(wv*wv*wv*wv);
//		dens *= max_volumedensity;

		// step size fac
//		dens *= mv_dens;

//		dens = max(dens,0.00001);
//		dens += 0.00001;
		// normalize
		dens /= max_volumedensity;
		dens *= 1.0;
//		dens *= rayleigh_fac;
		if(dens > 1.0){
//			return - 1.0;
		}
		bool hit = dens > hash_f();

		if(hit){
			return 1.0;
		}
		if(dot(p,p) > 6){
			return -1.0;
		}
	}
	return -1.0;
}



// Spectrum to xyz approx function from Sloan http://jcgt.org/published/0002/02/01/paper.pdf
// Inputs:  Wavelength in nanometers
float xFit_1931( float wave ) {
	float t1 = (wave-442.0)*((wave<442.0)?0.0624:0.0374),
	t2 = (wave-599.8)*((wave<599.8)?0.0264:0.0323),
	t3 = (wave-501.1)*((wave<501.1)?0.0490:0.0382);
	return 0.362*exp(-0.5*t1*t1) + 1.056*exp(-0.5*t2*t2)- 0.065*exp(-0.5*t3*t3);
}
float yFit_1931( float wave ) {
	float t1 = (wave-568.8)*((wave<568.8)?0.0213:0.0247),
	t2 = (wave-530.9)*((wave<530.9)?0.0613:0.0322);
	return 0.821*exp(-0.5*t1*t1) + 0.286*exp(-0.5*t2*t2);
}
float zFit_1931( float wave ) {
	float t1 = (wave-437.0)*((wave<437.0)?0.0845:0.0278),
	t2 = (wave-459.0)*((wave<459.0)?0.0385:0.0725);
	return 1.217*exp(-0.5*t1*t1) + 0.681*exp(-0.5*t2*t2);
}

#define xyzFit_1931(w) vec3( xFit_1931(w), yFit_1931(w), zFit_1931(w) )

vec3 XYZtosRGB( vec3 XYZ ) { mat3 m = mat3 ( 3.2404542, -1.5371385, -0.4985314, -0.9692660,  1.8760108,  0.0415560, 0.0556434, -0.2040259,  1.0572252 ); return XYZ * m; }
vec3 sRGBtoXYZ( vec3 RGB ) { mat3 m = mat3(  	0.4124564,  0.3575761, 0.1804375, 0.2126729,  0.7151522, 0.0721750, 0.0193339,  0.1191920, 0.9503041 ); return RGB * m; }


vec3 wv_to_xyz( float fWavelength ){ return xyzFit_1931(fWavelength); }

float t1(float x, float l, float m, float n, float o){ return n*(m*(x - l))/sqrt(1. + pow(m*(x-l),2.)) + o; }

float r1(float x){ return min( t1(x,589.,0.28,0.48,0.484) + 1. - t1(x,406.,0.078,0.17,0.83), 1. - t1(x,677.,0.048,0.325,0.349) ); }

float b1(float x){ return min( max( 1. - t1(x,491.5,0.145,0.475,0.508), t1(x,682.,0.0455,0.168,0.171) ), t1(x,408.5,0.083,0.331,0.645) ); }
float g1(float x){ return 1. - r1(x) - b1(x); }
float rgb_to_spectrum(vec3 rgb, float wv){ return r1(wv)*rgb.r + b1(wv)*rgb.b + g1(wv)*rgb.g; }

vec4 getLightStrengthD65(vec4 wv) {
	vec4 x = clamp(wv, 380.0, 780.0);
	return -10906.5 +
		x * (88.3842 +
		x * (-0.280601 +
		x * (0.000441847 +
		(-3.46017*1e-7+1.07868*1e-10 * x) * x)));
}


float A(float x){
	return 0.5*pow(x,0.2);
}
float B(float x){
	return 1. - 0.5*pow(1.-x,0.1-0.7*pow(1.-x,1.5));
}
float M(float x){
	return 0.31*pow(x,0.47-0.1*x - 4592.*pow(max((x-0.45)*0.26,0.),3.5)) + 0.205;
}

float inverse_CDF_Y(float x){
	return mix(
	mix(A(x),M(x),clamp((x-0.1)*13.,0.,1.)),
	B(x),
	clamp((x-0.45)*3.,0.,1.)
	);
}


// Spectrum with its integral
const vec3 cie[48] = vec3[48]( vec3( 0.000000f, 0.000000f, 0.000001f ),    vec3( 0.000006f, 0.000001f, 0.000026f ),    vec3( 0.000160f, 0.000017f, 0.000705f ),    vec3( 0.002362f, 0.000253f, 0.010482f ),    vec3( 0.019110f, 0.002004f, 0.086011f ),    vec3( 0.084736f, 0.008756f, 0.389366f ),    vec3( 0.204492f, 0.021391f, 0.972542f ),    vec3( 0.314679f, 0.038676f, 1.553480f ),    vec3( 0.383734f, 0.062077f, 1.967280f ),    vec3( 0.370702f, 0.089456f, 1.994800f ),    vec3( 0.302273f, 0.128201f, 1.745370f ),    vec3( 0.195618f, 0.185190f, 1.317560f ),    vec3( 0.080507f, 0.253589f, 0.772125f ),    vec3( 0.016172f, 0.339133f, 0.415254f ),    vec3( 0.003816f, 0.460777f, 0.218502f ),    vec3( 0.037465f, 0.606741f, 0.112044f ),    vec3( 0.117749f, 0.761757f, 0.060709f ),    vec3( 0.236491f, 0.875211f, 0.030451f ),    vec3( 0.376772f, 0.961988f, 0.013676f ),    vec3( 0.529826f, 0.991761f, 0.003988f ),    vec3( 0.705224f, 0.997340f, 0.000000f ),    vec3( 0.878655f, 0.955552f, 0.000000f ),    vec3( 1.014160f, 0.868934f, 0.000000f ),    vec3( 1.118520f, 0.777405f, 0.000000f ),    vec3( 1.123990f, 0.658341f, 0.000000f ),    vec3( 1.030480f, 0.527963f, 0.000000f ),    vec3( 0.856297f, 0.398057f, 0.000000f ),    vec3( 0.647467f, 0.283493f, 0.000000f ),    vec3( 0.431567f, 0.179828f, 0.000000f ),    vec3( 0.268329f, 0.107633f, 0.000000f ),    vec3( 0.152568f, 0.060281f, 0.000000f ),    vec3( 0.081261f, 0.031800f, 0.000000f ),    vec3( 0.040851f, 0.015905f, 0.000000f ),    vec3( 0.019941f, 0.007749f, 0.000000f ),    vec3( 0.009577f, 0.003718f, 0.000000f ),    vec3( 0.004553f, 0.001768f, 0.000000f ),    vec3( 0.002175f, 0.000846f, 0.000000f ),    vec3( 0.001045f, 0.000407f, 0.000000f ),    vec3( 0.000508f, 0.000199f, 0.000000f ),    vec3( 0.000251f, 0.000098f, 0.000000f ),    vec3( 0.000126f, 0.000050f, 0.000000f ),    vec3( 0.000065f, 0.000025f, 0.000000f ),    vec3( 0.000033f, 0.000013f, 0.000000f ),    vec3( 0.000018f, 0.000007f, 0.000000f ),    vec3( 0.000009f, 0.000004f, 0.000000f ),    vec3( 0.000005f, 0.000002f, 0.000000f ),    vec3( 0.000003f, 0.000001f, 0.000000f ),    vec3( 0.000002f, 0.000001f, 0.000000f )     );

const float invSpecIntegral = 1.0 / 116.643;

float lambdaSample() {
	float y = hash_f(), g;
	int i;
	for (i = 0; i < 47 && y >= 0.0; i++) {
		g = cie[i].y * invSpecIntegral;
		y -= g * 10.0;
	}
	return float(355 + i * 10) + y / g;
}

float lambdaPdf(float l) {
	int index = int(l * 0.1 - 35.5);
	if(index < 0 || index > 47) {
		return 0.0;
	}
	return cie[index].y * invSpecIntegral;
}



float blackbody(float l, float t) {
	const float h = 6.62607004e-16; const float k = 1.38064852e-5; const float c = 299792458e9;
	float a = 2.0 * h * c * c;
	float b = h * c / (l * k * t);
	return a / (l*l*l*l*l  * (exp(b) - 1.0));
}

float inters_sphere( vec3 ro, vec3 rd, vec3 ce, float ra ) {
	vec3 oc = ro - ce;
	float b = dot( oc, rd );
	float c = dot( oc, oc ) - ra*ra;
	float h = b*b - c;
	if( h<0.0 ) {return -1;} // no intersection
	h = sqrt( h );

	float da = -b-h;
	float db = -b+h;
	float d;
	if(da < 0. && db < 0.){
		return -1;
	} else if(da > 0. && db > 0.){
		d = min(da,db);
	} else if (da > 0.){
		d = da;
	} else {
		d = db;
	}

	return d;
}
float earth_diam_km = 12.742;
float max_atmo_height_km = 1.0;
float oxy_nitr_base_dens = 5.39*pow(10.0,18.0) + 2.01*pow(10.0,19.0);

// Temperature and pressure for standard air
const float Ts = 288.15f;   // Kelvin
const float Ps = 101325.0f; // Pa

// Temp K, 0-86km
float[] standard_atmosphere_temperature_lut = float[](
288.15f, 281.65f, 275.15f, 268.65f, 262.15f, 255.65f, 249.15f, 242.65f, 236.15f, 229.65f, 223.15f, 216.65f, 216.65f, 216.65f, 216.65f, 216.65f, 216.65f, 216.65f, 216.65f, 216.65f, 216.65f, 217.65f, 218.65f, 219.65f, 220.65f, 221.65f, 222.65f, 223.65f, 224.65f, 225.65f, 226.65f, 227.65f, 228.65f, 231.45f, 234.25f, 237.05f, 239.85f, 242.65f, 245.45f, 248.25f, 251.05f, 253.85f, 256.65f, 259.45f, 262.25f, 265.05f, 267.85f, 270.65f, 270.65f, 270.65f, 270.65f, 270.65f, 267.85f, 265.05f, 262.25f, 259.45f, 256.65f, 253.85f, 251.05f, 248.25f, 245.45f, 242.65f, 239.85f, 237.05f, 234.25f, 231.45f, 228.65f, 225.85f, 223.05f, 220.25f, 217.45f, 214.65f, 212.65f, 210.65f, 208.65f, 206.65f, 204.65f, 202.65f, 200.65f, 198.65f, 196.65f, 194.65f, 192.65f, 190.65f, 188.65f, 186.946f, 186.946f
);

/**
 * Rayleigh volume scattering coefficient beta_s (km^-1) for standard air
 * (Ps = 1013.25 mbars, Ts = 288.15K) and for wavelengths from 200nm to 4000nm.
 * A. Bucholtz 1995. Rayleigh-scattering calculations for the terrestrial atmosphere.
 * http://augerlal.lal.in2p3.fr/pmwiki/uploads/Bucholtz.pdf
 */
// 10 increment 200-800km
// 100 increment 800-2000km
// 200 increment 2000-3000km
// 500 increment 3000-4000km
const float[] rayleigh_volume_scattering_lut = float[]( 9.202e-1f, 7.225e-1f, 5.781e-1f, 4.691e-1f, 3.859e-1f, 3.208e-1f, 2.690e-1f, 2.277e-1f, 1.940e-1f, 1.665e-1f, 1.437e-1f, 1.249e-1f, 1.090e-1f, 9.557e-2f, 8.425e-2f, 7.450e-2f, 6.619e-2f, 5.901e-2f, 5.275e-2f, 4.734e-2f, 4.261e-2f, 3.845e-2f, 3.479e-2f, 3.156e-2f, 2.870e-2f, 2.616e-2f, 2.389e-2f, 2.187e-2f, 2.005e-2f, 1.842e-2f, 1.696e-2f, 1.564e-2f, 1.444e-2f, 1.336e-2f, 1.238e-2f, 1.149e-2f, 1.067e-2f, 9.927e-3f, 9.247e-3f, 8.624e-3f, 8.053e-3f, 7.530e-3f, 7.049e-3f, 6.605e-3f, 6.196e-3f, 5.819e-3f, 5.470e-3f, 5.146e-3f, 4.847e-3f, 4.568e-3f, 4.310e-3f, 4.069e-3f, 3.846e-3f, 3.637e-3f, 3.442e-3f, 3.261e-3f, 3.090e-3f, 2.931e-3f, 2.781e-3f, 2.641e-3f, 2.510e-3f, 1.561e-3f, 1.022e-3f, 6.964e-4f, 4.909e-4f, 3.560e-4f, 2.644e-4f, 2.005e-4f, 1.548e-4f, 1.214e-4f, 9.656e-5f, 7.775e-5f, 6.331e-5f, 4.322e-5f, 3.050e-5f, 2.214e-5f, 1.646e-5f, 1.249e-5f, 6.737e-6f, 3.948e-6f );


/**
 * Pressure (Pa) as a function of height (km) from the US Standard Atmosphere.
 * US COESA 1976. Standard Atmosphere, 1976. US Government Printing Office.
 * http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf
 */
// Pressure Pa, 0-86km
const float[] standard_atmosphere_pressure_lut = float[]( 101325.0f, 89874.6f, 79495.2f, 70108.5f, 61640.2f, 54019.9f, 47181.0f, 41060.7f, 35599.8f, 30742.5f, 26436.3f, 22632.1f, 19330.4f, 16510.4f, 14101.8f, 12044.6f, 10287.5f, 8786.68f, 7504.84f, 6410.01f, 5474.89f, 4677.89f, 3999.79f, 3422.43f, 2930.49f, 2511.02f, 2153.09f, 1847.46f, 1586.29f, 1362.96f, 1171.87f, 1008.23f, 868.019f, 748.228f, 646.122f, 558.924f, 484.317f, 420.367f, 365.455f, 318.220f, 277.522f, 242.395f, 212.030f, 185.738f, 162.937f, 143.135f, 125.910f, 110.906f, 97.7545f, 86.1623f, 75.9448f, 66.9389f, 58.9622f, 51.8668f, 45.5632f, 39.9700f, 35.0137f, 30.6274f, 26.7509f, 23.3296f, 20.3143f, 17.6606f, 15.3287f, 13.2826f, 11.4900f, 9.92203f, 8.55275f, 7.35895f, 6.31992f, 5.41717f, 4.63422f, 3.95642f, 3.37176f, 2.86917f, 2.43773f, 2.06792f, 1.75140f, 1.48092f, 1.25012f, 1.05351f, 0.88628f, 0.74428f, 0.623905f, 0.522037f, 0.435981f, 0.36342f, 0.302723f );

// 360-830
const float[] desert_dust_absorption_cross_section = float[]( 4.1273e-16f, 4.1287e-16f, 4.1298e-16f, 4.1312e-16f, 4.1345e-16f, 4.1373e-16f, 4.1380e-16f, 4.1492e-16f, 4.1562e-16f, 4.1699e-16f, 4.1868e-16f, 4.1810e-16f, 4.2080e-16f, 4.1986e-16f, 4.2291e-16f, 4.2402e-16f, 4.2776e-16f, 4.3358e-16f, 4.3448e-16f, 4.4695e-16f, 4.4645e-16f, 4.5019e-16f, 4.6188e-16f, 4.5954e-16f, 4.6552e-16f, 4.7575e-16f, 4.7034e-16f, 4.6794e-16f, 4.7937e-16f, 4.8121e-16f, 4.7506e-16f, 4.8124e-16f, 4.9612e-16f, 4.9827e-16f, 4.9147e-16f, 4.9394e-16f, 5.0729e-16f, 5.1219e-16f, 5.0358e-16f, 4.9739e-16f, 5.2279e-16f, 5.6035e-16f, 5.9223e-16f, 6.1249e-16f, 6.3642e-16f, 6.4141e-16f, 6.6077e-16f, 6.7588e-16f );
const float[] desert_dust_scattering_cross_section = float[]( 3.4215e-16f, 3.4240e-16f, 3.4266e-16f, 3.4288e-16f, 3.4308e-16f, 3.4320e-16f, 3.4319e-16f, 3.4296e-16f, 3.4234e-16f, 3.4110e-16f, 3.4069e-16f, 3.4022e-16f, 3.3966e-16f, 3.3900e-16f, 3.3824e-16f, 3.3598e-16f, 3.3303e-16f, 3.2921e-16f, 3.2435e-16f, 3.1835e-16f, 3.1454e-16f, 3.1042e-16f, 3.0589e-16f, 3.0092e-16f, 2.9559e-16f, 2.9413e-16f, 2.9275e-16f, 2.9137e-16f, 2.8999e-16f, 2.8859e-16f, 2.8490e-16f, 2.8105e-16f, 2.7707e-16f, 2.7296e-16f, 2.6875e-16f, 2.6727e-16f, 2.6584e-16f, 2.6441e-16f, 2.6301e-16f, 2.6143e-16f, 2.4192e-16f, 2.1835e-16f, 1.9026e-16f, 1.5663e-16f, 1.1720e-16f, 1.1560e-16f, 1.1457e-16f, 1.1362e-16f );
void get_fl_ceil_fr(float altitude, int max_range, out int fl_alt, out int ceil_alt, out float fr_alt ){
	fl_alt = int(altitude);
	ceil_alt = fl_alt + 1;
	fl_alt = min(fl_alt, max_range);
	ceil_alt = min(ceil_alt, max_range);
	fr_alt = fract(altitude);
}

void get_fl_ceil_fr_map_range(
	float altitude,
	float min_range,
	float max_range,
	out int fl_alt,
	out int ceil_alt,
	out float fr_alt
){
	altitude -= min_range;
	altitude /= max_range;
	fl_alt = int(altitude);
	ceil_alt = fl_alt + 1;
//	fl_alt = min(fl_alt, max_range);
//	ceil_alt = min(ceil_alt, max_range);
	fr_alt = fract(altitude);
}

float get_atmosphere_temp(float altitude){
	int fl_alt; int ceil_alt; float fr_alt;
	get_fl_ceil_fr(altitude, 86, fl_alt, ceil_alt, fr_alt );
	float a = standard_atmosphere_temperature_lut[fl_alt]; float b = standard_atmosphere_temperature_lut[ceil_alt];
	return mix(a,b,fr_alt);
}
float get_atmosphere_pressure(float altitude){
	int fl_alt; int ceil_alt; float fr_alt;
	get_fl_ceil_fr(altitude, 86, fl_alt, ceil_alt, fr_alt );
	float a = standard_atmosphere_pressure_lut[fl_alt]; float b = standard_atmosphere_pressure_lut[ceil_alt];
	return mix(a,b,fr_alt);
}
float get_rayleigh_scattering_lut(float altitude){
	// 10 increment 200-800km
	// 100 increment 800-2000km
	// 200 increment 2000-3000km
	// 500 increment 3000-4000km
	float min_range = 200;
	float max_range = 800;
	int offs = 0;
	if(altitude >= 3000){
		min_range = 3000; max_range = 4000;
		offs += 39 + 12 + 5;
	} else if (altitude >= 2000){
		min_range = 2000; max_range = 3000;
		offs += 39 + 12;
	} else if (altitude >= 800){
		min_range = 800; max_range = 2000;
		offs += 39;
	}
//	if(altitude >= 3000)

	int fl_alt; int ceil_alt; float fr_alt;
	get_fl_ceil_fr_map_range( altitude, min_range, max_range, fl_alt, ceil_alt, fr_alt );

	fl_alt += offs;
	ceil_alt += offs;
	ceil_alt = min(ceil_alt, rayleigh_volume_scattering_lut.length());
	fl_alt = min(ceil_alt, rayleigh_volume_scattering_lut.length());
	float a = rayleigh_volume_scattering_lut[fl_alt]; float b = rayleigh_volume_scattering_lut[ceil_alt];
	return mix(a,b,fr_alt);
}

float get_molecular_density_at_alt(float alt_km){
	return oxy_nitr_base_dens * exp(-alt_km/7.994);
}

float get_molecular_scattering(float height, float wl) {
//	height *= 1e-3f; // To km
//	float beta_s = lut_lerp(rayleigh_volume_scattering_lut, wl); // km^-1
	float beta_s = get_rayleigh_scattering_lut(height);
	float T = get_atmosphere_temp(height);
	float P = get_atmosphere_pressure(height);
	float beta = beta_s * (P / Ps) * (Ts / T); // km^-1
	return beta * 1e-3f; // m^-1
}

//float
//GuimeraAtmosphere::get_molecular_absorption(float height, float wl) const
//{
//	float sigma_a = lut_lerp(ozone_cross_section_lut, wl) * 1e-4f; // m^2 / molecules
//	// 1 Dobson = 2.6867e20 molecules / m^2
//	float total_ozone = ozone_mean_monthly_dobson[_month] * 2.6867e20f; // molecules / m^-2
//	float density = get_ozone_height_distribution(height) * total_ozone; // molecules / m^-2
//	density /= 9e3f; // m^-3
//	return sigma_a * density; // m^-1
//}

float trace_atmo(inout vec3 p, vec3 rd, float wv, bool to_sun){
	float d_to_earth_km = inters_sphere( p, rd, vec3(0), earth_diam_km/2.0);
	float d_to_atmo_end_km = abs(inters_sphere( p, rd, vec3(0), earth_diam_km/2.0 + max_atmo_height_km));

	float dist_km = d_to_earth_km > 0.0 && d_to_earth_km < d_to_atmo_end_km ? d_to_earth_km : d_to_atmo_end_km;

	if(to_sun && d_to_earth_km < 0){
		return -1.0;
	}

	float rayleigh_fac = wv*wv*wv*wv;
	rayleigh_fac = 1/(rayleigh_fac);
	rayleigh_fac /= 2.0e-13;
	rayleigh_fac *= 0.005;

	rayleigh_fac = 1/(wv*wv*wv*wv);
	float mv_dens = max_volumedensity * rayleigh_fac;
	float step_sz = 0.01;
	for(int i = 0; i < 1450; i++){
//		float t = -log(1.-hash_f()) / mv_dens;
		p = p + rd * step_sz;
		float alt_km = length(p) - earth_diam_km/2.;

		if(alt_km > max_atmo_height_km || alt_km < 0.0){
			if(to_sun){
				return 1.0;
			} else {
				return -1.0;
			}
		}

		float scatt_pct = get_molecular_scattering(alt_km, 0.0);
		float dens = scatt_pct;
		dens *= 1000;
		dens *= 100000000000.;
//		dens *= 100.;
//		dens = get_molecular_density_at_alt(alt_km)*0.00000000000000000001;
		dens = exp(-alt_km/7.994);
		dens *= 56141111114.;
//		dens *= 0.4;

//		dens *= 0.2;
		dens *= rayleigh_fac;
//		dens =
//		dens *= 0.01;
//		bool hit = dens > hash_f();
		bool hit = hash_f() > (exp(-dens*step_sz));

		if(hit){
			return 1.0;
		}
	}
	return -1.0;
}




void main(){
//	vec4 C = vec4(0);
	seed = 125125125;
	seed += F;
	if(U.x > R.x/2){
		seed += hashi(uint(U.x) + F*52u) + hashi(uint(U.y) + F*2);
	} else {
		seed += uint(blue_noise(U, 3)*1000000);
	}


	bool debug = U.x/R.x > 0.5;
	bool sample_Y = true;
//	sample_Y = false;

	if(pt_frame < 10){
//		sample_Y = false;
	}


	float wv = mix(
		standardObserver1931_w_min,
		standardObserver1931_w_max,
//		debug ? inverse_CDF_Y(1.-hash_f()) : hash_f()
		hash_f()
	);
	if(sample_Y){
		wv = lambdaSample();
	}

	float wv_idx = (wv - standardObserver1931_w_min)/(standardObserver1931_w_max - standardObserver1931_w_min);


	vec3 ro;
	vec3 rd;
	get_ro_and_rd(ro, rd, uv, false);

	vec3 col = vec3(0);


    vec4 prev_frame = tex(prev_frame_tex, uvn);
//	col = B(U).xxx;

	vec3 p = ro;
//	vec3 thr = vec3(1);
	float thr = 1.0;
	float delta_dist = 0.005;

	p += rd*delta_dist * blue_noise(U, 0)*4.;

	float sun_col = blackbody(wv, 6400.0)*0.005;


	float acc = 0.;

	vec3 sun_pos = normalize(vec3(1,1.6,0))*1000;

	if(false){
//		float t = 0.;
//		for(int i = 0; i < 1450; i++){
//			vec4 vol_samp = get_vol(p);
//			float absorbance = exp(-delta_dist*vol_samp.x);
//			float transmittance = 1.0-absorbance;
//			//		thr *= 1 - d;
//			//		thr = mix(vec3(0), thr,1.-d);
//			if(hash_f() > absorbance){
//				rd = sample_sph();
//				thr *= vol_samp.yzw;
//				thr /= tau;
//			}
//			if(t > 4){
//				break;
//			}
//			//		if(d > 0.){
//			//			thr = vec3(0);
//			//		}
//
//			t += delta_dist;
//			p += rd * delta_dist;
//		}
	} else {
		float t = 0.;
		if(false){
			for (int bnc = 0; bnc < 40; bnc++){
				//			vec3 q = p;
				vec4 vol_samp;
				float t = trace_vol(p, vol_samp, rd, wv);
				bool hit = t >= 0;

				if (hit){
					vec3 prev_rd = rd;
					//				thr *= vol_samp.yzw;
					thr *= vol_samp.y;

					{
						// NEE
						vec3 q = p;
						vec4 nee_vol_samp;
						//					vec3 light_dir = normalize(vec3(1));
						vec3 light_dir = normalize(sun_pos - p);

						//					float brdf = volumeBRDFPDF(dot(prev_rd,light_dir), anis);

						float brdf = rayleigh_brdf(dot(prev_rd, light_dir));

						float trace_chance = brdf;
						trace_chance = max(trace_chance, 0.0);
						//					trace_chance = pow(trace_chance, 4.);
						float t = trace_vol(q, nee_vol_samp, light_dir, wv);
						bool nee_hit = t < 0.;
						if (nee_hit){
							acc += thr * brdf* sun_col;
						}
					}

					if (U.x/R.x < 0.5){
						if (false){
							vec3 s = sample_henyey(anis);
							float brdf = henyey_brdf(s.z, anis);
							float pdf = henyey_greenstein_pdf(s.z, anis);
							mat3 m = get_orth_mat(vec3(0), rd);
							rd = m * s;
						} else {
							vec3 s = rayleigh_sample();
							float brdf = rayleigh_brdf(s.z);
							float pdf = rayleigh_brdf(s.z);
							mat3 m = get_orth_mat(vec3(0), rd);
							rd = m * s;
						}

						bnc++;
					} else {
						rd = sample_sph();
						//					float brdf = volumeBRDFPDF(dot(prev_rd,rd), anis);
						//					brdf = max(brdf,0.);
						float brdf = rayleigh_brdf(dot(prev_rd, rd));
						thr *= brdf/(1/(4.*pi));
					}
				}

				//			if(luma(thr) < 0.001){
				if (thr < 0.005){
					break;
				}
			}
		}
	}

//	col += thr * sample_env_map(rd)*5.;
//	vec3 env_samp_rgb = sample_env_map(rd);
//	env_samp_rgb *= 0.5;
//	float env = rgb_to_spectrum( env_samp_rgb, wv)*1.0;
//	env_samp_rgb *= 5.0;



	float d_to_earth_km = inters_sphere( ro, rd, vec3(0), earth_diam_km/2.0);
	float d_to_atmo_end_km = abs(inters_sphere( ro, rd, vec3(0), earth_diam_km/2.0 + max_atmo_height_km));

//	vec3 p = ro;

	if(d_to_earth_km >0.){
//		col += 0.5;
	}

	debug = all(ivec2(U) == ivec2(200,200));
	if(debug){
		dbg_print_buff[print_idx].cnt = 0;
	}

	int bnc;
	for(bnc = 0; bnc < 116; bnc++){
		float t = trace_atmo(p, rd, wv, false);
		bool hit = t > 0.0;
		if(hit){
			vec3 prev_rd = rd;
//			thr *= vol_samp.y;
			if(debug){
				print(bnc);
			}
			{
				// NEE
				vec3 q = p;
				vec4 nee_vol_samp;
				//					vec3 light_dir = normalize(vec3(1));
				vec3 light_dir = normalize(sun_pos - p);

				//					float brdf = volumeBRDFPDF(dot(prev_rd,light_dir), anis);

				float brdf = rayleigh_brdf(dot(prev_rd, light_dir));

//				float t = trace_vol(q, nee_vol_samp, light_dir, wv);
				float t = trace_atmo(q, light_dir, wv, true);
				bool nee_hit = t < 0.;
				if (nee_hit){
					acc += thr * brdf* sun_col*0.04;
				}
			}

			vec3 s = rayleigh_sample();
			float brdf = rayleigh_brdf(s.z);
			float pdf = rayleigh_brdf(s.z);
			mat3 m = get_orth_mat(vec3(0), rd);
			rd = m * s;

			thr *= 0.5;
//			col += 0.1;
		} else {
			break;
		}
	}
//	acc += thr * env;
	vec3 xyz = wv_to_xyz(wv);
	col = acc * xyz;




//	col *= getLightStrengthD65(wv.xxxx).x*0.0066;
	col *= 2.;
	if(sample_Y){
		col /= lambdaPdf(wv);
		//		col /= xyz.y;
		col *= 0.0024;
	}
	//	col += thr * sample_env_map(rd)*5.;



//	col += fract(uv.xyx + T);

	if(pt_frame == 0){
		prev_frame.rgb = vec3(0);
	}
    C.rgb = mix(prev_frame.rgb, col, 1./(1. + float(pt_frame)));
	C.w = 1.0;

//	C += 5.;
//	C.rgb = fract(U.xyx);
}
