vec2 rsi(vec3 pos, vec3 dir, float r) {
    float b     = dot(pos, dir);
    float det   = b*b - dot(pos, pos) + r*r;

    if (det < 0.0) return vec2(-1.0);

        det     = sqrt(det);

    return vec2(-b) + vec2(-det, det);
}

float rayleighPhase(float cosTheta) {
    return 3.0/(16.0*PI)*(1.0 + cosTheta*cosTheta);
}

float miePhase(float cosTheta, float g) {
    return (1-g*g)/(4.0*PI*pow(1.0 + g*g - 2*g*cosTheta, 1.5));
}

// float rayleighPhase(float cosTheta) {
// 	return (cosTheta * cosTheta + 1.0) * 0.75 / PI;
// }

// float miePhase(float cosTheta, float g) {
// 	float gg = g * g;
// 	float p1 = (0.375 * (1.0 - gg)) / (PI * (2.0 + gg));
// 	float p2 = (cosTheta * cosTheta + 1.0) * pow(-2.0 * g * cosTheta + 1.0 + gg, -1.5);
// 	return p1 * p2;
// }

float miePhaseCS(float cosTheta) {
	float g = mieg;
  	float gg = g*g;
  	float p1 = 3.0*(1.0 - gg) / (PI*(2.0 + gg));
  	float p2 = (1.0 + cosTheta*cosTheta) / pow((1.0 + gg - 2.0 * g * cosTheta), 3.0/2.0);
  	float phase = p1*p2;
        phase  /= 8.0;
  	return max0(phase);
}

float ozoneDensity(float elevation) {
	const float ozonePeakAlt = 3e4;
const float ozoneFalloff = 6e3;
const float rOzoneFalloff = 1.0 / ozoneFalloff;
const float ozoneFalloffTop = 18e3;
const float rOzoneFalloffTop = 1.0 / ozoneFalloffTop;
	// float elevationKm = elevation*0.001;
    // float o1 = 25.0 *     exp(( 0.0 - elevationKm) /   8.0) * 0.5;
    // float o2 = 30.0 * pow(exp((18.0 - elevationKm) /  80.0), elevationKm - 18.0);
    // float o3 = 75.0 * pow(exp((25.3 - elevationKm) /  35.0), elevationKm - 25.3);
    // float o4 = 50.0 * pow(exp((30.0 - elevationKm) / 150.0), elevationKm - 30.0);
    // return 2.0*(o1 + o2 + o3 + o4) / 134.628;
	float d   = exp(-max(0.0, (ozonePeakAlt - elevation) - planetRadius) * rOzoneFalloff);
    d  *= exp(-max(0.0, (elevation - ozonePeakAlt) - planetRadius) * rOzoneFalloffTop);
	return d*0.7;
}

vec3 getAtmosDensity(float h) {
    return vec3(exp(-h*rScaleHeights), ozoneDensity(h))*vec3(1.0, 1., 1.0);//*1.5
}

float getElevation(vec3 rayPos) {
	return length(rayPos) - planetRadius;
}

vec3 getAtmosOpticalDepth(vec3 rayPos, vec3 rayDir, float dist, const int steps) {
	const float rSteps = rcp(float(steps));

	float dd = dist * rSteps;
	vec3 rayStep = rayDir*dd;
	rayPos = rayPos + rayStep*(0.5 * (max0(rayDir.y) * 0.5 + 0.5));

	vec3 od = vec3(0.0);
	for(int i = 0; i < steps; i++) {
		vec3 density = getAtmosDensity(getElevation(rayPos));
		od += density * dd;

		rayPos += rayStep;
	}
	return atmosExtinctionMat * od;
}

vec3 getAtmosOpticalDepth(vec3 rayPos, vec3 rayDir, const int steps) {
	return getAtmosOpticalDepth(rayPos, rayDir, rsi(rayPos, rayDir, planetRadius + atmosHeight).y, steps);
}

vec3 getAtmosOpticalDepth(float cosTheta, float h) {
	return atmosExtinctionMat * texture2D(colortex1, vec2(cosTheta*0.5 + 0.5, h/atmosHeight)).rgb*10e3;
}