/*
    Sky: gradients, stars, sun/moon discs, aurora, 2D clouds for reflections.
    getSkyBase      : cheap analytic gradient (fog, reflections).
    getSkyCelestial : base + stars + sun/moon (no clouds).
    getSkyPhysical  : full Rayleigh+Mie sky with real altitude (PHYSICAL_SKY).
    getSkyAll       : celestial + 2D clouds, used by reflections and glass.
    Requires common.glsl, lighting.glsl (and atmosphere.glsl if PHYSICAL_SKY).
*/

// Analytic atmospheric gradient (fog, reflections - cheap)
vec3 getSkyBase(vec3 dir, vec3 sunDir) {
    float sh = sunDir.y;
    float dayF  = smoothstep(-0.05, 0.25, sh);
    // Asymmetric: below the horizon the sunset red dies fast (fixes /time set night)
    float duskF = clamp(1.0 - abs(sh) * 4.0, 0.0, 1.0) * smoothstep(-0.14, -0.04, sh);
    float horiz = pow(1.0 - clamp(dir.y, 0.0, 1.0), 3.0);

    // CONFIG_SKY dials (neutral at 1.0): horizon-band brightness vs zenith and
    // night luminance scale the base colors before the zenith/horizon blend.
    vec3 dayZen = vec3(0.13, 0.40, 0.92);
    vec3 dayHor = vec3(0.62, 0.80, 0.98) * SKY_HORIZON_LIFT;
    vec3 ngtZen = vec3(0.004, 0.007, 0.018) * SKY_NIGHT_BRIGHT;
    vec3 ngtHor = vec3(0.012, 0.018, 0.040) * (SKY_NIGHT_BRIGHT * SKY_HORIZON_LIFT);

    vec3 sky = mix(mix(ngtZen, dayZen, dayF), mix(ngtHor, dayHor, dayF), horiz);

    float sunAmt = pow(clamp(dot(dir, sunDir) * 0.5 + 0.5, 0.0, 1.0), 3.0);
    // Punchier horizon band: redder/more saturated, glow x1.4
    sky = mix(sky, vec3(1.0, 0.34, 0.09), clamp(duskF * sunAmt * horiz * 1.45, 0.0, 1.0));

    float glow = pow(clamp(dot(dir, sunDir), 0.0, 1.0), 8.0);
    sky += vec3(1.0, 0.80, 0.50) * glow * 0.22 * dayF;
    sky += vec3(1.0, 0.36, 0.10) * glow * 0.85 * duskF;

    sky = mix(sky, vec3(luma(sky)) * vec3(0.8, 0.85, 0.95), rainStrength * 0.85);

    // Final per-channel tint (CONFIG_SKY): sky dome, reflections and fogs all
    // follow, this gradient being the shared source. 1.0 = exact legacy sky.
    return sky * 1.22 * vec3(SKY_TINT_R, SKY_TINT_G, SKY_TINT_B);
}

// Procedural 2D clouds (reflections, glass, light mode)
#if CLOUD_HIGH_V3 == 1
// Smaller cells (2 scales), sharpened edges, crisper sun relief -> matches the
// volumetric cumulus inside reflections. Cost is nearly identical.
vec4 getClouds(vec3 dir, vec3 sunDir) {
    if (dir.y < 0.02) return vec4(0.0);

    vec2 p = dir.xz / dir.y * 2.4;
    p += cloudWorldTime() * vec2(0.012, 0.004); // world time: these 2D clouds follow /time set like the volumetric ones

    float n = fbm(p * 1.1) * 0.72 + fbm3(p * 3.1 + 11.0) * 0.31;
    float cov = smoothstep(0.52 - 0.25 * rainStrength, 0.80, n);
    if (cov <= 0.0) return vec4(0.0);

    vec2 sunOff = normalize(sunDir.xz + vec2(1e-4)) * 0.20;
    float n2 = fbm((p + sunOff) * 1.1) * 0.72 + fbm3((p + sunOff) * 3.1 + 11.0) * 0.31;
    float shade = clamp((n - n2) * 5.0 + 0.5, 0.0, 1.0);

    vec3 bright = getLightColor() * 0.50 + getAmbientColor() * 0.80;
    vec3 dark   = getAmbientColor() * 0.45;
    vec3 col = mix(bright, dark, clamp(shade * 0.75 + cov * 0.25, 0.0, 1.0));

    float a = cov * smoothstep(0.02, 0.12, dir.y) * 0.85;
    return vec4(col, a);
}
#else
vec4 getClouds(vec3 dir, vec3 sunDir) {
    if (dir.y < 0.02) return vec4(0.0);

    vec2 p = dir.xz / dir.y * 2.4;
    p += cloudWorldTime() * vec2(0.012, 0.004); // world time: these 2D clouds follow /time set like the volumetric ones

    float n = fbm(p * 0.7);
    n = n * 0.7 + fbm(p * 2.3 + 11.0) * 0.3;
    float cov = smoothstep(0.50 - 0.25 * rainStrength, 0.82, n);
    if (cov <= 0.0) return vec4(0.0);

    vec2 sunOff = normalize(sunDir.xz + vec2(1e-4)) * 0.25;
    float n2 = fbm((p + sunOff) * 0.7);
    float shade = clamp((n - n2) * 4.0 + 0.5, 0.0, 1.0);

    vec3 bright = getLightColor() * 0.42 + getAmbientColor() * 0.75;
    vec3 dark   = getAmbientColor() * 0.55;
    vec3 col = mix(bright, dark, clamp(shade * 0.7 + cov * 0.3, 0.0, 1.0));

    float a = cov * smoothstep(0.02, 0.12, dir.y) * 0.9;
    return vec4(col, a);
}
#endif

// Celestial disc: square (vanilla style) or round (realistic), per SQUARE_SUN
float celestialDisc(vec3 dir, vec3 bodyDir, float size) {
    float dd = dot(dir, bodyDir);
    if (dd < 0.5) return 0.0;
    #ifdef SQUARE_SUN
        vec3 axis = (abs(bodyDir.y) > 0.95) ? vec3(1.0, 0.0, 0.0) : vec3(0.0, 1.0, 0.0);
        vec3 right = normalize(cross(bodyDir, axis));
        vec3 up2 = cross(right, bodyDir);
        vec2 ang = vec2(dot(dir, right), dot(dir, up2)) / dd;
        float sq = max(abs(ang.x), abs(ang.y));
        return 1.0 - smoothstep(size * 0.92, size, sq);
    #else
        // Atmospheric refraction: near the horizon the disc grows and
        // flattens vertically, like real sunsets
        vec3 axisR = (abs(bodyDir.y) > 0.95) ? vec3(1.0, 0.0, 0.0) : vec3(0.0, 1.0, 0.0);
        vec3 rightR = normalize(cross(bodyDir, axisR));
        vec3 upR = cross(rightR, bodyDir);
        float horizS = 1.0 - clamp(abs(bodyDir.y) * 4.0, 0.0, 1.0);
        vec2 angR = vec2(dot(dir, rightR), dot(dir, upR)) / max(dd, 0.5);
        angR.y *= 1.0 + horizS * 0.50;
        float sizeR = size * (1.0 + horizS * 0.30);
        return 1.0 - smoothstep(sizeR * 0.90, sizeR, length(angR));
    #endif
}

#ifdef AURORA
/*
    Aurora: thin vertical curtains that ripple and snake.
    The sky is projected onto a plane; vertical "folds" (the curtain edge) are
    defined by an animated sinuous line. Several layers at different phases
    give the characteristic thickness and shimmer. Gradient green (bottom) ->
    cyan -> magenta (top).
*/
vec3 aurora(vec3 dir, float t, float nightF) {
    if (dir.y < 0.04 || nightF <= 0.0) return vec3(0.0);
    // Fades out in "space" view (camera > 350 blocks): from orbit the aurora
    // would sit below the atmosphere, not wrap all around the viewer
    float highFade = 1.0 - clamp((cameraPosition.y - 350.0) / 300.0, 0.0, 1.0);
    if (highFade <= 0.0) return vec3(0.0);
    // Variable activity: some nights blaze, others show nothing
    float act = smoothstep(0.30, 0.70, vnoise(vec2(t * 0.0035, 3.7)));
    if (act <= 0.001) return vec3(0.0);

    // Projection: x = azimuth, height modulates intensity
    vec2 p = dir.xz / (dir.y + 0.25);
    float hgrad = clamp((dir.y - 0.05) * 1.6, 0.0, 1.0);

    vec3 acc = vec3(0.0);
    // 5 stacked curtains, each with its own sinuosity and color
    for (int i = 0; i < 5; i++) {
        float fi = float(i);
        // Snaking curtain line (sum of animated sines)
        float wob = sin(p.x * 1.3 + t * 0.25 + fi * 1.7) * 0.55
                  + sin(p.x * 3.1 - t * 0.40 + fi * 4.2) * 0.22
                  + sin(p.x * 6.7 + t * 0.15 + fi) * 0.10;
        float curtainX = p.x * 0.5 + wob + fi * 0.6;
        // Distance to the fold -> thin bright edge
        float edge = exp(-pow(fract(curtainX) - 0.5, 2.0) * 60.0);
        // Vertical rays shimmering along the curtain
        float ray = 0.55 + 0.45 * sin(p.x * 24.0 + t * 0.8 + fi * 10.0
                                      + sin(p.x * 5.0 - t) * 3.0);
        // Color by height: deep green -> cyan -> magenta
        vec3 cc = mix(vec3(0.10, 1.00, 0.42),
                      mix(vec3(0.15, 0.85, 0.95), vec3(0.85, 0.25, 0.95), hgrad),
                      smoothstep(0.0, 0.7, hgrad));
        acc += cc * edge * ray * (0.5 + 0.5 * hash12(vec2(fi, floor(curtainX))));
    }
    // Fade: weak at the horizon, peak mid-height, dissolves toward zenith
    float vfade = smoothstep(0.04, 0.22, dir.y) * smoothstep(1.0, 0.45, dir.y);
    return acc * 0.085 * act * nightF * vfade * AURORA_BRIGHTNESS * highFade;
}
#endif

// Point star inside a 3D cell: sharp core + soft halo (legacy STARS_V2=0 path)
vec3 starInCell(vec3 cube, float scale, float thresh) {
    vec3 g = cube * scale;
    vec3 cell = floor(g);
    float s = hash13(cell);
    if (s < thresh) return vec3(0.0);
    // Random sub-cell position of the star
    vec3 off = vec3(hash13(cell + 1.3), hash13(cell + 2.7), hash13(cell + 4.1));
    float d = length(fract(g) - off);
    float h2 = hash13(cell + 17.3);
    float mag = (s - thresh) / (1.0 - thresh);          // brightness scales with rarity
    float core = exp(-d * d * 900.0) + exp(-d * d * 120.0) * 0.35; // core + halo
    float tw = 0.55 + 0.45 * sin(frameTimeCounter * (1.5 + h2 * 4.0) + h2 * 50.0);
    vec3 tint = mix(vec3(0.60, 0.74, 1.35), vec3(1.35, 0.92, 0.58), h2 * h2);
    return tint * core * tw * (0.25 + mag * mag * 2.2);
}

#if STARS_V2 == 1
// ==================== STAR FIELD V2 (STARS_V2) ====================
// Power-law magnitudes (rare bright stars, myriad faint ones), subtle stellar
// colors (blue-white -> white -> orange), clusters & voids (low-frequency
// density sampled at the CELL CENTER so no star is cut in half), per-star
// twinkle, and a gaussian core of FIXED ANGULAR width (~1.4 px): faint stars
// are DARKER, never sub-pixel -> no crawling under TAA.

// Rotate a vector about a unit axis (Rodrigues formula)
vec3 starRotate(vec3 v, vec3 ax, float ang) {
    float ca = cos(ang);
    return v * ca + cross(ax, v) * sin(ang) + ax * dot(ax, v) * (1.0 - ca);
}

// Stellar tint from a "temperature" index t:
// 0 = blue-white (O/B), ~0.4 = white (A/F), 1 = orange (K/M).
// Saturation kept deliberately low, as seen by the naked eye.
vec3 starTintV2(float t) {
    vec3 c = (t < 0.4)
        ? mix(vec3(0.60, 0.71, 1.00), vec3(1.00, 0.98, 0.95), t * 2.5)
        : mix(vec3(1.00, 0.98, 0.95), vec3(1.00, 0.76, 0.52), (t - 0.4) * 1.6667);
    return mix(vec3(1.0), c, 0.55);
}

// One star layer over 3D cells. keep = fraction of populated cells;
// clusterAmt = weight of the cluster/void field on that density;
// floorB/tailB = floor and tail of the magnitude power law.
vec3 starLayerV2(vec3 sd, float scale, float keep, float clusterAmt,
                 float floorB, float tailB, float twAmp) {
    vec3 g = sd * scale;
    vec3 cell = floor(g);

    // Clusters & voids: density sampled at the cell center (constant per cell)
    float keepEff = keep;
    if (clusterAmt > 0.0) {
        vec3 cc = (cell + 0.5) / scale;
        float dens = vnoise3(cc * 2.5) * 0.65 + vnoise3(cc * 6.0 + 11.0) * 0.35;
        keepEff *= mix(1.0, 0.45 + smoothstep(0.28, 0.75, dens) * 1.30, clusterAmt);
    }
    float s = hash13(cell);
    if (s > keepEff) return vec3(0.0);

    // Sub-cell position, confined to limit clipping at cell borders
    vec3 off = vec3(hash13(cell + 1.3), hash13(cell + 2.7), hash13(cell + 4.1)) * 0.7 + 0.15;
    float dAng = length(fract(g) - off) / scale; // ANGULAR distance (rad)

    // Power-law magnitude: most sit at the floor, bright ones are rare;
    // the step() kicker adds a "first magnitude" star (~1 in 300 populated)
    float hm = hash13(cell + 8.9);
    float flux = floorB + pow(hm, 9.0) * tailB + step(0.9965, hm) * tailB * 0.85;

    // Temperature: bright stars skew blue-white (realistic)
    float ht = hash13(cell + 13.7);
    float t = mix(ht, ht * 0.45, clamp(flux * 0.8, 0.0, 1.0));

    // Gaussian core of fixed angular width (~1.4 px at 1440p / FOV 70, scaled
    // by STAR_SIZE) + a soft halo reserved for bright stars (credible PSF that
    // feeds the bloom). No renormalization on purpose: a wider star IS a bit
    // brighter when integrated, which is what makes it read as bigger.
    float sigA = 7.0e-4 * STAR_SIZE;
    float sig2 = 2.0 * sigA * sigA;
    float x2 = dAng * dAng;
    float core = exp(-x2 / sig2);
    float halo = exp(-x2 / (sig2 * 9.0)) * 0.10 * clamp(flux - 0.4, 0.0, 2.0);

    // Discreet twinkle, per-star phase and speed
    float h2 = hash13(cell + 17.3);
    float tw = 1.0 - twAmp * (0.5 + 0.5 * sin(frameTimeCounter * (0.9 + h2 * 2.6) + h2 * 61.0));

    return starTintV2(t) * (core + halo) * flux * tw;
}
#endif

// Deep sky: multi-scale stars, Milky Way, shooting stars
vec3 skyStars(vec3 dir, float visibility) {
    if (visibility <= 0.0) return vec3(0.0);

    // Sky vault rotation about the north axis. The mod is taken on the exact
    // rotation PERIOD (2pi/0.000262 ~ 23984.7 ticks) BEFORE multiplying: on an
    // old world (worldDay>5000) float ULP reached 8-16 ticks and the vault
    // stuttered. The wrap is invisible (the jump is exactly 2pi).
    float rotA = mod(float(worldDay) * 24000.0 + float(worldTime), 23984.68) * 0.000262; // continuous time (no jump at dawn)
    #if STARS_V2 == 1
    // V2: the vault turns about a tilted CELESTIAL POLE (~38 deg above the
    // north horizon), like a real mid-latitude sky - stars trace circles
    // centered on the pole, not on the zenith.
    vec3 sd = starRotate(dir, normalize(vec3(0.0, 0.616, -0.788)), rotA);
    #else
    mat3 spin = mat3(cos(rotA), 0.0, sin(rotA), 0.0, 1.0, 0.0, -sin(rotA), 0.0, cos(rotA));
    vec3 sd = spin * dir;
    #endif

    #if STARS_V2 == 1
    // Stronger twinkle near the horizon (more air crossed)
    float horizF = 1.0 - smoothstep(0.0, 0.35, dir.y);
    float twAmp = 0.10 + 0.24 * horizF;
    // 3 layers: rare bright stars (uniformly spread, like real nearby stars),
    // medium and fine ones grouped into clusters/voids
    vec3 col = (starLayerV2(sd, 24.0, 0.032, 0.0, 0.12, 3.0, twAmp)
             +  starLayerV2(sd, 60.0, 0.090, 1.0, 0.060, 0.90, twAmp)
             +  starLayerV2(sd, 105.0, 0.140, 1.0, 0.038, 0.30, twAmp)) * STAR_BRIGHTNESS;
    #else
    // Legacy star field: 3 densities (rare large, medium, fine)
    vec3 col = (starInCell(sd, 26.0, 0.94)
             + starInCell(sd, 55.0, 0.96) * 0.7
             + starInCell(sd, 110.0, 0.975) * 0.4) * STAR_BRIGHTNESS;
    #endif

    // Milky Way: galactic axis, projected plane
    vec3 mwN = normalize(vec3(0.36, 0.28, 0.89));
    float dN = dot(sd, mwN);
    float band = exp(-dN * dN * 5.5);                // wide gaussian band
    vec3 mwU = normalize(cross(mwN, vec3(0.0, 1.0, 0.0)));
    vec3 mwV = cross(mwN, mwU);
    vec2 mwp = vec2(dot(sd, mwU), dot(sd, mwV));

    // Multi-octave cloud structure along the band
    vec2 q = vec2(mwp.x * 1.2, mwp.y * 5.0);
    float neb = fbmTex(q * 1.5 + 7.0) * 0.6 + fbmTex(q * 4.0 - 3.0) * 0.3 + fbmTex(q * 9.0) * 0.1;
    // Dark dust lanes striping the band.
    // PERF: the dust fbm is skipped when the Milky Way is off
    // (MILKYWAY_BRIGHTNESS=0.0 by default -> folded at compile time)
    float dust = (MILKYWAY_BRIGHTNESS > 0.001) ? smoothstep(0.40, 0.72, fbmTex(q * 2.2 + vec2(0.0, 4.0))) : 0.0;
    // Galactic bulge: warm swelling near the center
    float core = exp(-dot(mwp, mwp) * 1.4);
    float mw = band * smoothstep(0.25, 0.85, neb) * (1.0 - dust * 0.7);
    // Tint: blue-violet arms -> golden cream core
    vec3 mwCol = mix(vec3(0.42, 0.50, 0.95), vec3(1.10, 0.90, 0.68), core);
    col += (mwCol * mw * 0.5 + mwCol * core * band * 0.35) * MILKYWAY_BRIGHTNESS;
    #if STARS_V2 == 1
    // Denser myriad of fine stars INSIDE the band + a subtle "unresolved" glow
    // (the Milky Way IS a myriad of stars). Independent of MILKYWAY_BRIGHTNESS,
    // which only drives the nebulosity above.
    col += starLayerV2(sd, 150.0, 0.170, 0.6, 0.030, 0.12, twAmp) * band * 2.6 * STAR_BRIGHTNESS;
    col += mwCol * band * (0.25 + neb * 0.75) * 0.028 * STAR_BRIGHTNESS;
    #else
    // Denser star myriad inside the galactic band
    col += starInCell(sd, 150.0, 0.965) * band * 1.6 * STAR_BRIGHTNESS;
    #endif

    #ifdef SHOOTING_STARS
    {
        // Shooting star on a GREAT CIRCLE path (no planar-projection warping).
        // Spawn window ~2.86 s, lifetime 0.86 s, variable magnitude.
        float sB = floor(frameTimeCounter * 0.35);
        float sT = fract(frameTimeCounter * 0.35);
        if (hash12(vec2(sB, 11.3)) > 0.55) {
            float ph = sT * 3.326; // 0..1 over the lifetime (0.86 s), >1 after (off)
            if (ph < 1.0) {
                // Great-circle axis + orthonormal basis (random per window)
                vec3 axN = normalize(vec3(hash12(vec2(sB, 5.1)) - 0.5,
                                          0.25 + hash12(vec2(sB, 9.7)) * 0.6,
                                          hash12(vec2(sB, 3.3)) - 0.5));
                vec3 axU = normalize(cross(axN, vec3(0.0, 1.0, 0.0)) + vec3(1e-4, 0.0, 0.0));
                vec3 axV = cross(axN, axU);
                float a0g = hash12(vec2(sB, 7.7)) * 6.2831853;
                float sweep = 0.55 + hash12(vec2(sB, 2.9)) * 0.5; // arc ~31-60 deg
                float aH = a0g + ph * sweep;
                vec3 headP = axU * cos(aH) + axV * sin(aH);
                // ANGULAR distances: off the circle plane / along the arc
                float perpA = dot(sd, axN);
                float aP = atan(dot(sd, axV), dot(sd, axU));
                float dA = aH - aP;
                dA = dA - 6.2831853 * floor(dA / 6.2831853 + 0.5); // wrap [-pi,pi]
                float mag = 0.6 + hash12(vec2(sB, 4.2)) * 0.9;     // variable magnitude
                float fade = sin(clamp(ph, 0.0, 1.0) * 3.14159);
                // Min angular width ~1.5 px at 1440p (avoids TAA shimmer);
                // tapered trail with no hard step (smooth onset at the head).
                float trail = exp(-perpA * perpA * 250000.0)
                            * exp(-max(dA, 0.0) * 22.0)
                            * smoothstep(-0.004, 0.004, dA)
                            * (1.0 - smoothstep(0.20, 0.35, dA));
                vec3 relH = sd - headP;
                float headG = exp(-dot(relH, relH) * 90000.0)
                            + exp(-dot(relH, relH) * 9000.0) * 0.35;
                col += vec3(0.95, 0.97, 1.15) * (trail * 1.1 + headG * 2.6) * fade * mag;
            }
        }
    }
    #endif

    #if STARS_V2 == 1
    // Atmospheric extinction: near the horizon stars dim and warm slightly
    // (air mass), then fade out.
    // NOTE: fade on dir.y (the REAL horizon) - sd.y is tilted by the V2 pole.
    col *= mix(vec3(1.0), vec3(0.80, 0.62, 0.45), horizF * 0.55);
    col *= 1.0 - horizF * 0.45;
    return col * visibility * smoothstep(-0.06, 0.12, dir.y);
    #else
    return col * visibility * smoothstep(-0.06, 0.12, sd.y);
    #endif
}

/*
    Realistic moon: spherical disc with procedural maria, limb darkening and
    REAL PHASES (the terminator follows moonPhase: full, quarters, crescents,
    new moon). Stays square vanilla-style under SQUARE_SUN.
*/
vec3 moonSurface(vec3 dir, vec3 moonDir, float size) {
    float dd = dot(dir, moonDir);
    // Soft halo around the moon
    vec3 glow = vec3(0.55, 0.65, 0.90) * pow(clamp(dd, 0.0, 1.0), 900.0) * 0.75;
    #ifdef MOON_HALO
        // 22-degree halo: ice ring on humid nights
        float ring = exp(-pow((dd - 0.927) * 55.0, 2.0));
        glow += vec3(0.55, 0.62, 0.85) * ring * 0.05
              * (0.25 + 0.75 * clamp(wetness + rainStrength * 0.5, 0.0, 1.0));
    #endif
    if (dd < 0.5) return vec3(0.0);

    #ifdef SQUARE_SUN
        return glow + vec3(0.75, 0.82, 1.0) * celestialDisc(dir, moonDir, size) * 1.4;
    #else
        vec3 axis = (abs(moonDir.y) > 0.95) ? vec3(1.0, 0.0, 0.0) : vec3(0.0, 1.0, 0.0);
        vec3 right = normalize(cross(moonDir, axis));
        vec3 up2 = cross(right, moonDir);
        vec2 ang = vec2(dot(dir, right), dot(dir, up2)) / dd / size;
        float r2 = dot(ang, ang);
        float disc = 1.0 - smoothstep(0.80, 1.0, r2);
        if (disc <= 0.0) return glow;

        // The disc is a sphere: local normal + sun direction rotating with the
        // phase -> correct terminator
        vec3 n = vec3(ang, sqrt(max(1.0 - r2, 0.0)));
        float pa = float(moonPhase) * 0.7853982;
        vec3 Lm = vec3(sin(pa), 0.0, cos(pa));
        float lit = smoothstep(-0.05, 0.20, dot(n, Lm));

        // Lunar maria (dark patches) + limb darkening
        float maria = vnoise(ang * 3.5 + 7.3) * 0.45 + vnoise(ang * 8.0 + 2.1) * 0.25;
        vec3 surf = vec3(0.85, 0.88, 0.97) * (0.85 - maria * 0.45);
        surf *= 0.75 + 0.25 * n.z;

        return glow + surf * lit * disc * 2.4;
    #endif
}

// Analytic atmosphere + stars + sun/moon (NO clouds) - light fallback
vec3 getSkyCelestial(vec3 dir, vec3 sunDir) {
    vec3 sky = getSkyBase(dir, sunDir);
    float dayF = smoothstep(-0.05, 0.25, sunDir.y);

    // Day gate (audited): once sunDir.y >= 0.25 the visibility below is
    // EXACTLY 0 and skyStars early-outs before any noise work; the aurora
    // call further down early-outs on nightF <= 0 the same way. Daytime sky
    // and reflection pixels pay neither - no extra sunHeight gate needed.
    sky += skyStars(dir, (1.0 - dayF) * (1.0 - rainStrength) * 0.85);

    float sunDisc = celestialDisc(dir, sunDir, SUN_SIZE);
    sky += vec3(1.0, 0.93, 0.82) * sunDisc * 25.0 * dayF * (1.0 - rainStrength * 0.9);

    sky += moonSurface(dir, -sunDir, SUN_SIZE * 0.8)
         * (1.0 - dayF * 0.8) * (1.0 - rainStrength * 0.9);

    #ifdef AURORA
        sky += aurora(dir, frameTimeCounter, (1.0 - dayF) * (1.0 - rainStrength));
    #endif

    return sky;
}

#ifdef PHYSICAL_SKY
// Full physical sky: Rayleigh+Mie, real altitude, stars, celestial discs
vec3 getSkyPhysical(vec3 dir, vec3 sunDir) {
    float camH = atmCameraHeight(cameraPosition.y);

    vec3 trans;
    vec3 sky = atmScatter(dir, sunDir, camH, trans) * ATMO_INTENSITY;

    // Moonlight: second pass, faint and cold, at dusk/night
    float sh = sunDir.y;
    if (sh < 0.15) {
        vec3 mtrans;
        sky += atmScatter(dir, -sunDir, camH, mtrans)
             * vec3(0.55, 0.75, 1.0) * 0.035 * ATMO_INTENSITY;
    }

    // Stars: visible when the air lets them through (night or altitude).
    // Day gate (audited): daytime skyGlow drives starVis to exactly 0 and
    // skyStars early-outs; the aurora nightFA gate below does the same.
    float skyGlow = luma(sky);
    float starVis = luma(trans) * clamp(1.0 - skyGlow * 6.0, 0.0, 1.0) * (1.0 - rainStrength);
    sky += skyStars(dir, starVis * 1.1);

    // Sun: attenuated by the air crossed -> red at sunset, crisp at altitude
    float sunDisc = celestialDisc(dir, sunDir, SUN_SIZE);
    sky += vec3(1.0, 0.98, 0.94) * sunDisc * 40.0 * trans;

    // Moon boosted, slightly blue-tinted silver
    sky += moonSurface(dir, -sunDir, SUN_SIZE * 0.8) * trans * 1.85 * vec3(0.97, 1.0, 1.08);

    // More dramatic sunrise/sunset: soft glow amplification around the body +
    // pushed saturation at dusk, applied ON TOP of the PHYSICAL result
    // (accentuation, not replacement).
    float duskB = clamp(1.0 - abs(sh) * 4.0, 0.0, 1.0) * smoothstep(-0.12, -0.02, sh);
    if (duskB > 0.001) {
        float nearSun = pow(clamp(dot(dir, sunDir), 0.0, 1.0), 3.0);
        // Warm TINT rather than intensity GAIN (a gain boost clips to white
        // through ACES) + strong saturation push.
        sky *= mix(vec3(1.0), vec3(1.12, 0.78, 0.48), duskB * nearSun * 0.85);
        vec3 satS = sky - vec3(luma(sky));
        sky = max(sky + satS * duskB * 0.55, vec3(0.0));
    }

    #ifdef AURORA
        float nightFA = (1.0 - smoothstep(-0.10, 0.10, sh)) * (1.0 - rainStrength);
        sky += aurora(dir, frameTimeCounter, nightFA);
    #endif

    // Rain: grey sky
    sky = mix(sky, vec3(luma(sky)) * vec3(0.8, 0.85, 0.95), rainStrength * 0.85);

    // CONFIG_SKY creator tint: applied on the physical dome too, so it never
    // diverges from the getSkyBase consumers (reflections, fogs).
    sky *= vec3(SKY_TINT_R, SKY_TINT_G, SKY_TINT_B);
    return sky;
}
#endif

// Full sky with 2D clouds (used by reflections and glass)
vec3 getSkyAll(vec3 dir, vec3 sunDir) {
    vec3 sky = getSkyCelestial(dir, sunDir);

    #ifdef CLOUDS_2D
        vec4 cl = getClouds(dir, sunDir);
        sky = mix(sky, cl.rgb, cl.a);
    #endif

    return sky;
}
