/*
    Physical atmosphere. Nishita single scattering: Rayleigh (sky blue, red
    sunsets) + Mie (white glow around the sun), raymarched at planet scale
    (R = 6371 km). Altitude-aware: high above the world the air thins, the
    sky darkens and the atmosphere becomes a blue rim at the horizon.
    No buffers; pure math helpers used by lib/sky.glsl (PHYSICAL_SKY).
*/

const float ATM_RPLANET = 6371.0e3;
const float ATM_RTOP    = 6471.0e3;
const float ATM_HR      = 8500.0;   // Rayleigh scale height
const float ATM_HM      = 1300.0;   // Mie scale height

vec3 atmBetaR() { return vec3(5.8e-6, 13.5e-6, 33.1e-6); }
float atmBetaM() { return 21.0e-6; }

// Exit distance of a sphere of radius R (origin = planet center)
float atmRaySphere(vec3 ro, vec3 rd, float R) {
    float b = dot(ro, rd);
    float c = dot(ro, ro) - R * R;
    float h = b * b - c;
    if (h < 0.0) return -1.0;
    h = sqrt(h);
    float t = -b - h;
    if (t > 0.0) return t;
    return -b + h;
}

// Optical depth (Rayleigh, Mie) from p towards the sun - 3 samples
vec2 atmLightOD(vec3 p, vec3 sunDir) {
    float t = atmRaySphere(p, sunDir, ATM_RTOP);
    if (t < 0.0) return vec2(1.0e9);
    float b = dot(p, sunDir);
    float c = dot(p, p) - ATM_RPLANET * ATM_RPLANET;
    if (b < 0.0 && b * b - c > 0.0) return vec2(1.0e9);

    float ds = t / 3.0;
    vec2 od = vec2(0.0);
    vec3 q = p + sunDir * ds * 0.5;
    for (int j = 0; j < 3; j++) {
        float h = max(length(q) - ATM_RPLANET, 0.0);
        od += vec2(exp(-h / ATM_HR), exp(-h / ATM_HM)) * ds;
        q += sunDir * ds;
    }
    return od;
}

/*
    Single scattering along the view ray, 10 samples.
    camH: camera altitude in meters. trans: view transmittance (out), used to
    attenuate sun/moon/stars seen through the air.
*/
vec3 atmScatter(vec3 dir, vec3 sunDir, float camH, out vec3 trans) {
    vec3 ro = vec3(0.0, ATM_RPLANET + 200.0 + camH, 0.0);
    trans = vec3(1.0);

    // [entry, exit] interval inside the atmosphere shell: from space the march
    // starts AT THE ENTRY point so the limb stays visible
    float bb = dot(ro, dir);
    float cc = dot(ro, ro) - ATM_RTOP * ATM_RTOP;
    float hh = bb * bb - cc;
    if (hh < 0.0) return vec3(0.0); // ray misses the atmosphere
    float hs = sqrt(hh);
    float tIn = max(-bb - hs, 0.0);
    float tOut = -bb + hs;
    if (tOut <= tIn) return vec3(0.0);

    // Ground stops the ray
    float cg = dot(ro, ro) - ATM_RPLANET * ATM_RPLANET;
    float hg = bb * bb - cg;
    if (hg > 0.0) {
        float tg = -bb - sqrt(hg);
        if (tg > tIn) tOut = min(tOut, tg);
    }

    float ds = (tOut - tIn) / 10.0;
    vec3 p = ro + dir * (tIn + ds * 0.5);
    vec2 odView = vec2(0.0);
    vec3 sumR = vec3(0.0);
    vec3 sumM = vec3(0.0);

    for (int i = 0; i < 10; i++) {
        float h = max(length(p) - ATM_RPLANET, 0.0);
        vec2 dens = vec2(exp(-h / ATM_HR), exp(-h / ATM_HM)) * ds;
        odView += dens;
        vec2 odSun = atmLightOD(p, sunDir);
        vec3 att = exp(-atmBetaR() * (odView.x + odSun.x)
                       - atmBetaM() * 1.1 * (odView.y + odSun.y));
        sumR += att * dens.x;
        sumM += att * dens.y;
        p += dir * ds;
    }

    trans = exp(-atmBetaR() * odView.x - atmBetaM() * 1.1 * odView.y);

    float mu = dot(dir, sunDir);
    float phR = 0.0596831 * (1.0 + mu * mu);
    float g = 0.78;
    float phM = 0.1193662 * (1.0 - g * g) * (1.0 + mu * mu)
              / ((2.0 + g * g) * pow(1.0 + g * g - 2.0 * g * mu, 1.5));

    return (sumR * atmBetaR() * phR * RAYLEIGH_STRENGTH + sumM * vec3(atmBetaM()) * phM * MIE_STRENGTH) * 28.0;
}

// Generous altitude curve: nearly invisible at normal play heights, the sky
// only darkens very high up - stratosphere arrives smoothly near build limit.
float atmCameraHeight(float camY) {
    float dh = max(camY - 63.0, 0.0);
    float high = max(dh - 220.0, 0.0);
    return (dh * 8.0 + high * high * 8.0) * ATMO_ALTITUDE_BOOST;
}
