

vec2 rsi(vec3 r0, vec3 rd, float sr) {
    // ray-sphere intersection that assumes
    // the sphere is centered at the origin.
    // No intersection when result.x > result.y
    float a = dot(rd, rd);
    float b = 2.0 * dot(rd, r0);
    float c = dot(r0, r0) - (sr * sr);
    float d = (b*b) - 4.0*a*c;
    if (d < 0.0) return vec2(1e5,-1e5);
    return vec2(
        (-b - sqrt(d))/(2.0*a),
        (-b + sqrt(d))/(2.0*a)
    );
}

#define PLANET_RADIUS 6365e3
#define ATMOSPHERE_RADIUS 6481e3

const vec3 RAYLEIGH_COEFF  = vec3(6.05e-6, 1.25e-5, 3.02e-5);//vec3(6.05e-6, 1.25e-5, 3.02e-5);
const vec3 MIE_COEFF       = vec3(8.8e-6);
const vec3 OZONE_COEFF     = vec3(4.113e-7, 3.120e-7, 1.340e-8)*2.0;//vec3(3.58768136e-7, 3.14953105e-7, 5.43182681e-8)*2.0;
float RAYLEIGH_SCALE = mix(8.5e3, 13.5e3, sunset);
const float MIE_SCALE      = 1.2e3;
const float OZONE_HEIGHT   = 35e3;

const float densitySeaLevel = 1.225;
const mat3x3 EXCTINCT_MAT = mat3x3(RAYLEIGH_COEFF, MIE_COEFF * 1.11, OZONE_COEFF);
vec2 scaleHeights = vec2(RAYLEIGH_SCALE, MIE_SCALE);
vec2 inverseScaleHeights = 1.0 / vec2(RAYLEIGH_SCALE, MIE_SCALE);

const float g = 0.74;


vec3 approx_r_ho(vec3 rd) {
    // Approximate densities 
    // x - rayleigh
    // y - mie
    // z - ozone
    float y = rd.y * PLANET_RADIUS;
    vec4 atmosRadius = vec4(RAYLEIGH_SCALE, MIE_SCALE, OZONE_HEIGHT, OZONE_HEIGHT + RAYLEIGH_SCALE) + PLANET_RADIUS;
    vec4 z = sqrt(atmosRadius*atmosRadius + y*y - PLANET_RADIUS*PLANET_RADIUS) - y;

    return vec3(z.xy, (z.w - z.z) * 18.0);
}

float chapmanApproximation(in float x, in float h, in float mu) {
    float c = sqrt(x - h);
    float ch = c * expf(-h);

    if(mu >= 0.0) {
        return ch / (c * mu + 1.0);
    } else {
        float x0 = sqrt(1.0 - mu * mu) * (x + h);
        float c0 = sqrt(x0);

        return 2.0 * c0 * expf(x - x0) - ch / (1.0 - c * mu);
    }
}

vec3 calculateAnalyticalDensity(in float h, float mu) {
    float r = scaleHeights.x * chapmanApproximation(PLANET_RADIUS * inverseScaleHeights.x, h * inverseScaleHeights.x, mu) * densitySeaLevel;
    float m = scaleHeights.y * chapmanApproximation(PLANET_RADIUS * inverseScaleHeights.y, h * inverseScaleHeights.y, mu);
    float o = scaleHeights.x * chapmanApproximation((PLANET_RADIUS+30.0) * inverseScaleHeights.x, (h+30.0) * inverseScaleHeights.x, mu) * (mu*40.0);
    // float o = scaleHeights.x * chapmanApproximation((PLANET_RADIUS) * inverseScaleHeights.x, h * inverseScaleHeights.x, mu) * (mu*40.0);
    return vec3(r, m, o);
}

vec3 calculateAtmosphereTransmittanceAnalytical(in vec3 rayVector, in vec3 position, in mat3x3 ec, out vec3 density) {
    float rl = length(position);
    float h = rl - PLANET_RADIUS;
    vec3 r = position / rl;

    float cosChiRay = dot(r, rayVector);

    density = calculateAnalyticalDensity(h, cosChiRay);
    return expf(-ec * density);
}

vec3 rayleighPhase(float cos_theta, vec3 coeff) {
    return (3.0*coeff*(1.0 + cos_theta*cos_theta))/(16*PI);
}


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


vec3 F_ex(float dist) {
    return expf(-(RAYLEIGH_COEFF + MIE_COEFF)*dist);
}

vec3 calcSunTransmitted(vec3 dirToSun, vec3 L_0) {
     vec3 sun_r_ho = approx_r_ho(dirToSun);
     vec3 airdensity_InScatter = RAYLEIGH_COEFF*sun_r_ho.x + MIE_COEFF*sun_r_ho.y + OZONE_COEFF*sun_r_ho.z; // The probability of a scatter event 
     vec3 sunExtinction = expf(-airdensity_InScatter);
     return L_0 * sunExtinction;
}

vec3 calcSky(vec3 rd, vec3 dirToSun, vec3 E_sun) {
    float cos_theta = dot(rd, dirToSun);

    vec3 rPhase = rayleighPhase(cos_theta, RAYLEIGH_COEFF);
    vec3 mPhase = miePhase(cos_theta);

    vec3 r_ho = approx_r_ho(rd); 

    vec3 airdensity = EXCTINCT_MAT * r_ho; // Airmass                                        
    vec3 transmittance = expf(-airdensity);

    vec3 sun_r_ho = approx_r_ho(dirToSun); 
    vec3 airdensity_InScatter = EXCTINCT_MAT * sun_r_ho; // The probability of a scatter event 
    vec3 transmittance_InScatter = expf(-airdensity_InScatter);

    vec3 phases = (r_ho.x*rPhase + r_ho.y*mPhase);

    return E_sun * ATMOSPHERE_INTENSITY * phases * abs(transmittance_InScatter - transmittance + 1e-35) / (abs(airdensity_InScatter - airdensity) + 1e-35);
}

vec3 calcSky2(vec3 rd, vec3 dirToSun, vec3 E_sun) {
    rd.y = max(rd.y, -0.0);
    dirToSun.y = max(dirToSun.y, -0.1);
    float cos_theta = dot(rd, dirToSun);

    vec3 r0 = vec3(0, PLANET_RADIUS + 1.0, 0); // Eye is 1 meters above sea level
    float rl = length(r0);
    float h = rl - PLANET_RADIUS;
    vec3 r = r0 / rl;
    float cosChiRay = dot(r, rd);
    float cosChiSunRay = dot(r, dirToSun);

    vec3 rPhase = rayleighPhase(cos_theta, RAYLEIGH_COEFF);
    vec3 mPhase = miePhase(cos_theta);

    // vec3 r_ho = approx_r_ho(rd); 
    vec3 r_ho = calculateAnalyticalDensity(h, cosChiRay);

    
    

    vec3 airdensity = EXCTINCT_MAT * r_ho; // Airmass                                        
    vec3 transmittance = expf(-airdensity);

    vec3 sun_r_ho = calculateAnalyticalDensity(h, cosChiSunRay);//approx_r_ho(dirToSun); 
    vec3 airdensity_InScatter = EXCTINCT_MAT * sun_r_ho; // The probability of a scatter event 
    vec3 transmittance_InScatter = expf(-airdensity_InScatter);

    vec3 phases = r_ho.x*rPhase + r_ho.y*mPhase;

    return E_sun * ATMOSPHERE_INTENSITY * phases * abs(transmittance_InScatter - transmittance) / (abs(airdensity_InScatter - airdensity) + 1e-16);
}

vec3 sunDisk(vec3 sunDir, vec3 viewDir) {
    return float(1.0 - dot(sunDir, viewDir) < 0.00017*SUN_SIZE) * shadowLightIrradiance * 1000.0;
}