#ifndef ASTRAL_COSMIC_RIFT_GLSL
#define ASTRAL_COSMIC_RIFT_GLSL

uniform sampler2D astralCosmosCycleAtlas;
uniform sampler2D astralRiftCataclysm;

/* Aurora Astral - twin-ribbon cosmic band
 *
 * The geometry deliberately follows Aurora Astral's horizon projection and
 * horizontal advection.  Two independently phased aurora waves drive the
 * blue and violet boundaries, so the opening expands and contracts locally
 * instead of moving through the sky as one rigid piece.
 *
 * The dimensional galaxy is an independent direction-space volume evaluated
 * on the sky sphere.  The signed mask between the moving boundaries is the
 * only connection between galaxy and rift geometry.  No coordinate stretch,
 * fit or edge distortion is used; every cosmic layer is clipped at the mask.
 *
 * RGB is premultiplied radiance and alpha is displaced sky energy:
 *     sky = sky * (1.0 - rift.a) + rift.rgb;
 */

float rift_hash12(vec2 p) {
    vec3 p3 = fract(vec3(p.xyx) * 0.1031);
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.x + p3.y) * p3.z);
}

vec2 rift_hash22(vec2 p) {
    vec3 p3 = fract(vec3(p.xyx)
        * vec3(0.1031, 0.1030, 0.0973));
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.xx + p3.yz) * p3.zy);
}

float rift_noise(vec2 p) {
    vec2 cell = floor(p);
    vec2 f = fract(p);
    f = f * f * (3.0 - 2.0 * f);
    return mix(
        mix(rift_hash12(cell),
            rift_hash12(cell + vec2(1.0, 0.0)), f.x),
        mix(rift_hash12(cell + vec2(0.0, 1.0)),
            rift_hash12(cell + vec2(1.0, 1.0)), f.x), f.y);
}

float rift_fbm(vec2 p) {
    float result = 0.0;
    float weight = 0.5;
    #if COSMIC_RIFT_QUALITY == 1
        const int riftOctaves = 3;
    #elif COSMIC_RIFT_QUALITY == 2
        const int riftOctaves = 4;
    #else
        const int riftOctaves = 5;
    #endif
    for (int octave = 0; octave < riftOctaves; octave++) {
        result += rift_noise(p) * weight;
        p = mat2(1.82, -0.63, 0.63, 1.82) * p
            + vec2(13.17, 7.91);
        weight *= 0.5;
    }
    return result;
}

float rift_ridged(vec2 p) {
    return 1.0 - abs(rift_fbm(p) * 2.0 - 1.0);
}

float rift_hash13(vec3 p3) {
    p3 = fract(p3 * 0.1031);
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.x + p3.y) * p3.z);
}

float rift_noise3(vec3 p) {
    vec3 cell = floor(p);
    vec3 f = fract(p);
    f = f * f * (3.0 - 2.0 * f);

    float n000 = rift_hash13(cell);
    float n100 = rift_hash13(cell + vec3(1.0, 0.0, 0.0));
    float n010 = rift_hash13(cell + vec3(0.0, 1.0, 0.0));
    float n110 = rift_hash13(cell + vec3(1.0, 1.0, 0.0));
    float n001 = rift_hash13(cell + vec3(0.0, 0.0, 1.0));
    float n101 = rift_hash13(cell + vec3(1.0, 0.0, 1.0));
    float n011 = rift_hash13(cell + vec3(0.0, 1.0, 1.0));
    float n111 = rift_hash13(cell + vec3(1.0, 1.0, 1.0));

    float lower = mix(
        mix(n000, n100, f.x), mix(n010, n110, f.x), f.y);
    float upper = mix(
        mix(n001, n101, f.x), mix(n011, n111, f.x), f.y);
    return mix(lower, upper, f.z);
}

float rift_fbm3(vec3 p) {
    float result = 0.0;
    float weight = 0.5;
    #if COSMIC_RIFT_QUALITY == 1
        const int cosmicOctaves = 3;
    #elif COSMIC_RIFT_QUALITY == 2
        const int cosmicOctaves = 4;
    #else
        const int cosmicOctaves = 5;
    #endif
    for (int octave = 0; octave < cosmicOctaves; octave++) {
        result += rift_noise3(p) * weight;
        p = mat3(
             0.00,  0.80,  0.60,
            -0.80,  0.36, -0.48,
            -0.60, -0.48,  0.64) * p * 1.91
            + vec3(11.7, 7.3, 19.1);
        weight *= 0.5;
    }
    return result;
}

// Continuous upper-hemisphere mapping used only by discrete particle cells.
// The nebula itself stays fully three-dimensional in normalized direction.
vec2 rift_octahedral_sky(vec3 direction) {
    vec3 dir = normalize(direction);
    vec2 uv = dir.xz / max(
        abs(dir.x) + abs(dir.y) + abs(dir.z), 0.0001);
    if (dir.y < 0.0) {
        uv = (1.0 - abs(uv.yx)) * sign(uv.xy);
    }
    return uv;
}

// Smooth circular weights let three different cosmic fields dissolve into
// one another without a cut at either the frame boundaries or cycle wrap.
float rift_cycle_weight(float phase, float center) {
    float distanceToCenter = abs(
        fract(phase - center + 0.5) - 0.5) * 3.0;
    return 1.0 - smoothstep(0.0, 1.0, distanceToCenter);
}

// Three 2048px seamless fields are packed horizontally with four-pixel
// gutters.  Addressing their pixel centers prevents linear filtering from
// leaking one scene into another while retaining a single sampler unit.
vec2 rift_cycle_atlas_uv(vec2 uv, float frameIndex) {
    vec2 localUv = fract(uv);
    float pixelX = frameIndex * 2056.0 + 4.5
        + localUv.x * 2047.0;
    float pixelY = 0.5 + localUv.y * 2047.0;
    return vec2(pixelX / 6168.0, pixelY / 2048.0);
}

// One shared, twice-smoothed weather envelope drives every catastrophic
// layer. Keeping it independent from all animation clocks prevents a weather
// change from seeking to a different point on the rift's timeline.
float rift_weather_transition() {
    float rain = smoothstep(0.015, 0.985,
        clamp(rainStrength, 0.0, 1.0));
    return rain * rain * (3.0 - 2.0 * rain);
}

float rift_visibility() {
    float night = day_blend_float(0.0, 0.0, 1.0);
    float rainVeil = rift_weather_transition();
    // Rain no longer erases the phenomenon. At night it pushes the rift into
    // its catastrophic state, while daylight rain still receives no rift.
    return night * mix(1.0, 0.86, rainVeil);
}

float rift_cataclysm_amount() {
    return day_blend_float(0.0, 0.0, 1.0)
        * rift_weather_transition();
}

// Continuous mirrored addressing hides the edge of the generated storm map
// without relying on driver-specific custom-texture wrapping modes.
vec2 rift_cataclysm_uv(vec2 uv) {
    vec2 folded = fract(uv * 0.5) * 2.0;
    return 1.0 - abs(folded - 1.0);
}

void rift_cataclysm_boundaries(
        vec2 p, float textureClock, float fractureClock,
        float cataclysmAmount, inout float lowerEdge,
        inout float upperEdge, out float catastrophePressure) {
    catastrophePressure = 0.5 + 0.5 * sin(
        textureClock * 0.47 + rift_noise(vec2(p.x * 0.34, 71.0)) * 5.0);
    float lowerStormTear = rift_fbm(vec2(
        p.x * 1.72 - fractureClock * 0.11,
        113.0 + sin(p.x * 0.43) * 0.7)) - 0.5;
    float upperStormTear = rift_fbm(vec2(
        p.x * 2.13 + fractureClock * 0.13,
        197.0 - cos(p.x * 0.37) * 0.8)) - 0.5;
    float violentLower = sin(p.x * 5.7 - fractureClock * 0.33)
        * (0.025 + rift_noise(vec2(p.x * 0.8, 29.0)) * 0.035);
    float violentUpper = sin(p.x * 6.4 + fractureClock * 0.29 + 1.7)
        * (0.025 + rift_noise(vec2(p.x * 0.7, 53.0)) * 0.035);
    lowerEdge -= cataclysmAmount
        * (0.075 + catastrophePressure * 0.095);
    upperEdge += cataclysmAmount
        * (0.075 + (1.0 - catastrophePressure) * 0.095);
    lowerEdge += cataclysmAmount
        * (lowerStormTear * 0.26 + violentLower);
    upperEdge += cataclysmAmount
        * (upperStormTear * 0.26 + violentUpper);
}

// This is the original Aurora Astral projection.  It maps an endless wave
// over the sky hemisphere, hiding its continuation behind the horizon rather
// than introducing visible start and end points.
vec2 rift_astral_projection(vec3 viewDir, float astralTime) {
    float rawDenom = viewDir.y + 0.15;
    #ifdef DISTANT_RENDER_MOD
        const float horizonFloor = 0.2;
    #else
        const float horizonFloor = 0.05;
    #endif
    float denom = horizonFloor
        + rawDenom * rawDenom / (abs(rawDenom) + horizonFloor);
    vec2 p = viewDir.xz / denom;

    // Aurora-like sky migration slowly changes the rift orientation so its
    // two spectral sides can enter through different screen edges instead of
    // remaining a vertical object locked to the middle of the sky.
    float orientation = sin(astralTime * 0.55) * 0.30
        + sin(astralTime * 0.21 + 1.7) * 0.14;
    float orientationCos = cos(orientation);
    float orientationSin = sin(orientation);
    p = mat2(
        orientationCos, -orientationSin,
        orientationSin,  orientationCos) * p;

    // Horizontal advection is inherited from the original aurora.  It moves
    // the band through the sky without changing its apparent altitude.
    p.x += astralTime * 0.24;
    return p;
}

float rift_independent_astral_path(
    vec2 p,
    float astralTime,
    float clock,
    float ribbonIndex
) {
    // Aurora Astral gives every ribbon its own deterministic seed and drift.
    // The two rift boundaries retain that independence at geometry level,
    // rather than merely changing their surface brightness independently.
    float seed = ribbonIndex * 17.0;
    float drift = seed * 437.0
        + astralTime * (1.10 + ribbonIndex * 0.37);
    float ribbonMix = ribbonIndex - 1.0;
    float primaryFrequency = mix(0.20, 0.255, ribbonMix);
    float counterFrequency = mix(0.77, 0.96, ribbonMix);
    float counterDrift = mix(0.41, 0.59, ribbonMix);

    float path = sin(p.x * primaryFrequency + drift) * 0.38;
    path += sin(p.x * counterFrequency
        - drift * counterDrift) * 0.125;

    // Opposite flow vectors ensure the large silhouette of one boundary does
    // not mirror the other, while keeping the movement slow and aurora-like.
    float flowDirection = ribbonIndex < 1.5 ? -1.0 : 1.0;
    float broadFlow = rift_fbm(vec2(
        p.x * mix(0.17, 0.22, ribbonMix)
            + clock * flowDirection * mix(0.013, 0.017, ribbonMix),
        seed * 0.37
            + clock * mix(0.0043, -0.0036, ribbonMix))) - 0.5;
    float turbulence = rift_noise(vec2(
        p.x * mix(0.53, 0.66, ribbonMix)
            + clock * flowDirection * mix(0.058, 0.071, ribbonMix),
        seed + clock * mix(0.021, -0.018, ribbonMix))) - 0.5;
    path += broadFlow * 0.32;
    path += turbulence * 0.15;
    path += sin(p.x * mix(2.25, 2.73, ribbonMix)
        + clock * flowDirection * mix(0.16, 0.13, ribbonMix)
        + seed) * 0.044;
    return path;
}

void rift_twin_boundaries(
    vec2 p,
    float astralTime,
    float clock,
    out float lowerEdge,
    out float upperEdge
) {
    float sharedSkyWander = (rift_fbm(vec2(
        p.x * 0.12 + clock * 0.0057,
        88.3 + clock * 0.0031)) - 0.5) * 0.22;
    float openingNoise = (rift_fbm(vec2(
        p.x * 0.18 - clock * 0.0079,
        143.7 + clock * 0.0047)) - 0.5) * 0.28;
    float migrationNoise = rift_fbm(vec2(
        p.x * 0.055 + clock * 0.0021,
        206.4 + clock * 0.0017)) - 0.5;
    // A wide, multi-period migration lets the opening travel through the
    // sides of the sky instead of orbiting one central anchor.  The spatial
    // waves bend different stretches by different amounts, so this is shape
    // travel rather than a rigid translation of the whole rift.
    float skyMigration = sin(
        clock * 0.018 + migrationNoise * 2.5) * 0.82
        + migrationNoise * 0.95;
    skyMigration += sin(
        clock * 0.0071 + migrationNoise * 1.3 + 0.8) * 0.32;
    float flexibleSweep = sin(
        p.x * 0.105 + clock * 0.0093
            + migrationNoise * 1.8) * 0.30;
    flexibleSweep += sin(
        p.x * 0.24 - clock * 0.0061 + 2.4) * 0.14;
    skyMigration += flexibleSweep;
    float blueMigration = (rift_fbm(vec2(
        p.x * 0.075 - clock * 0.0037,
        271.3 + clock * 0.0013)) - 0.5) * 0.40;
    float violetMigration = (rift_fbm(vec2(
        p.x * 0.069 + clock * 0.0032,
        347.9 - clock * 0.0015)) - 0.5) * 0.40;
    float halfOpening = (0.66 + openingNoise)
        * COSMIC_RIFT_SCALE;

    float blueRaw = rift_independent_astral_path(
        p, astralTime, clock, 1.0)
        + sharedSkyWander + skyMigration
        + blueMigration - halfOpening;
    float violetRaw = rift_independent_astral_path(
        p, astralTime, clock, 2.0)
        + sharedSkyWander + skyMigration
        + violetMigration + halfOpening;

    // Start from two genuinely independent target paths.  A soft repulsion is
    // applied only near contact so colors never swap sides or form a numerical
    // zero-width seam.  Away from a pinch both paths retain their own motion.
    float midpoint = (blueRaw + violetRaw) * 0.5;
    float rawGap = violetRaw - blueRaw;
    float minimumGap = 0.16 * COSMIC_RIFT_SCALE;
    float maximumGap = 1.94 * COSMIC_RIFT_SCALE;
    float softness = 0.065 * COSMIC_RIFT_SCALE;

    float floorDelta = rawGap - minimumGap;
    float stableGap = minimumGap + 0.5 * (
        floorDelta + sqrt(floorDelta * floorDelta
            + softness * softness));
    stableGap = 0.5 * (stableGap + maximumGap
        - sqrt((stableGap - maximumGap)
            * (stableGap - maximumGap) + softness * softness));
    stableGap = clamp(stableGap, minimumGap, maximumGap);

    lowerEdge = midpoint - stableGap * 0.5;
    upperEdge = midpoint + stableGap * 0.5;
}

// Piecewise-linear noise deliberately keeps its interpolation unsmoothed.
// It supplies the crisp angular direction changes expected from a fracture,
// while the ordinary smooth noise remains available for gas and atmosphere.
float rift_linear_noise1(float x, float seed) {
    float cell = floor(x);
    float local = fract(x);
    float first = rift_hash12(vec2(cell, seed));
    float second = rift_hash12(vec2(cell + 1.0, seed));
    return mix(first, second, local) * 2.0 - 1.0;
}

// Signed fracture geometry for a single polar boundary.  Positive values are
// part of the dimensional opening and negative values are ordinary sky.  The
// result is unioned with the main band before any mask, rim or glow is made,
// so a fracture is literally a deformation of the rift rather than a layer.
float rift_integrated_fracture(
    float alongCoordinate,
    float outwardDistance,
    float clock,
    float sideSeed
) {
    float outward = max(outwardDistance, 0.0);
    float cellSize = 3.15 + rift_hash12(
        vec2(sideSeed, 8.7)) * 0.70;
    // A coherent animated skeleton carries all roots belonging to this lane.
    // Its low-frequency bend and drift move the fracture family together,
    // while every branch still inherits the parent boundary's full motion.
    float skeletonThread = rift_linear_noise1(
        alongCoordinate * 0.23 + clock * 0.006
            + sideSeed * 0.11,
        sideSeed + 101.3);
    float skeletonOffset = skeletonThread * 0.22
        + sin(clock * 0.009 + sideSeed * 0.73) * 0.14;
    float cellCoordinate = (alongCoordinate + skeletonOffset)
        / cellSize;
    float cellId = floor(cellCoordinate);
    float localAlong = (fract(cellCoordinate) - 0.5) * cellSize;

    vec2 identity = rift_hash22(vec2(cellId, sideSeed));
    float forkIdentity = rift_hash12(vec2(
        cellId + 7.3, sideSeed * 2.17));
    float secondForkIdentity = rift_hash12(vec2(
        cellId - 11.9, sideSeed * 3.41));

    // Roots remain fixed in boundary space.  Motion changes extension,
    // angular stress and fork opening, so the crack grows from its attachment
    // instead of sliding across the rift as an independent object.
    float presence = smoothstep(0.50, 0.66, identity.x);
    float growthPhase = 0.5 + 0.5 * sin(
        clock * (0.018 + identity.y * 0.012)
            + cellId * 1.37 + sideSeed);
    float growth = mix(0.58, 1.0,
        smoothstep(0.06, 0.94, growthPhase));
    float branchLength = mix(0.56, 1.24, identity.y) * growth;
    float normalizedOutward = outward / max(branchLength, 0.001);

    float rootAnchor = (rift_hash12(vec2(
        cellId * 2.3, sideSeed + 19.1)) - 0.5) * 0.78;
    rootAnchor += sin(
        clock * 0.022 + cellId * 0.77 + sideSeed) * 0.055;
    rootAnchor += skeletonThread * 0.070;
    float stemSway = sin(
        clock * 0.037 + cellId * 0.91 + sideSeed) * 0.15;
    float stemCenter = rootAnchor + outward * stemSway;
    float angularStem = rift_linear_noise1(
        normalizedOutward * mix(4.2, 6.2, forkIdentity)
            + clock * 0.010
            + cellId * 1.31,
        sideSeed + 31.7);
    stemCenter += angularStem * normalizedOutward
        * mix(0.11, 0.19, forkIdentity);

    float mainDistance = abs(localAlong - stemCenter);
    float mainWidth = mix(
        0.105, 0.0025,
        pow(clamp(normalizedOutward, 0.0, 1.0), 0.68));
    float mainSigned = mainWidth - mainDistance;
    mainSigned = min(mainSigned,
        (1.0 - normalizedOutward) * 0.16);

    float forkStart = mix(0.30, 0.47, forkIdentity);
    float forkProgress = max(
        outward - branchLength * forkStart, 0.0);
    float forkSign = forkIdentity < 0.5 ? -1.0 : 1.0;
    float forkCenter = stemCenter + forkSign * forkProgress
        * mix(0.55, 0.92, secondForkIdentity)
        * mix(0.70, 1.0, growth);
    float angularFork = rift_linear_noise1(
        normalizedOutward * 5.7
            - clock * 0.008 + cellId * 2.17,
        sideSeed + 73.1);
    forkCenter += angularFork * forkProgress * 0.18;
    float forkDistance = abs(localAlong - forkCenter);
    float forkWidth = mainWidth * 0.78;
    float forkSigned = forkWidth - forkDistance;
    forkSigned = min(forkSigned,
        (normalizedOutward - forkStart) * 0.12);
    forkSigned = min(forkSigned,
        (0.96 - normalizedOutward) * 0.13);

    // A rarer, shorter counter-fork prevents every event from reading as the
    // same Y silhouette while preserving plenty of clean boundary space.
    float secondForkGate = smoothstep(
        0.50, 0.60, normalizedOutward)
        * (1.0 - smoothstep(0.77, 0.92, normalizedOutward))
        * smoothstep(0.67, 0.91, secondForkIdentity);
    float secondForkProgress = max(
        outward - branchLength * 0.52, 0.0);
    float secondForkCenter = stemCenter - forkSign
        * secondForkProgress * mix(0.42, 0.70, forkIdentity)
        * mix(0.72, 1.0, growth);
    float secondForkDistance = abs(
        localAlong - secondForkCenter);
    float secondForkWidth = mainWidth * 0.62;
    float secondForkSigned = secondForkWidth - secondForkDistance;
    secondForkSigned = min(secondForkSigned,
        (normalizedOutward - 0.52) * 0.11);
    secondForkSigned = min(secondForkSigned,
        (0.88 - normalizedOutward) * 0.12);
    secondForkSigned -= (1.0 - secondForkGate) * 0.45;

    float fractureSigned = max(
        mainSigned, max(forkSigned, secondForkSigned));
    // A deep overlap buries the old polar border beneath the dimensional
    // body.  No colored separator can remain between parent and branch.
    fractureSigned = min(fractureSigned,
        outwardDistance + 0.120);
    fractureSigned -= (1.0 - presence) * 0.75;
    return fractureSigned;
}

float rift_star_layer(
    vec2 p,
    float threshold,
    float clock,
    float seed
) {
    vec2 cell = floor(p);
    vec2 localPosition = fract(p) - 0.5;
    vec2 offset = (rift_hash22(cell + seed) - 0.5) * 0.72;
    float starSeed = rift_hash12(cell + seed * 3.71);
    float pointDistance = length(localPosition - offset);
    float starPresence = smoothstep(
        threshold, min(threshold + 0.040, 0.999), starSeed);
    float radius = mix(0.014, 0.030,
        smoothstep(threshold, 1.0, starSeed));
    float pointAA = clamp(fwidth(pointDistance), 0.004, 0.070);
    float core = 1.0 - smoothstep(
        radius - pointAA, radius + pointAA, pointDistance);
    float halo = 1.0 - smoothstep(
        radius * 1.4, radius * 4.2 + pointAA, pointDistance);
    float twinkle = 0.68 + 0.32 * sin(
        clock * (0.75 + starSeed * 1.6) + starSeed * 43.0);
    return (core + halo * 0.14) * starPresence * twinkle;
}

vec3 rift_nebula_palette(float phase, float colorAxis) {
    vec3 midnightBlue = vec3(0.018, 0.055, 0.22);
    vec3 astralBlue = vec3(0.055, 0.29, 0.67);
    vec3 violet = vec3(0.32, 0.075, 0.59);
    vec3 magenta = vec3(0.60, 0.075, 0.43);
    vec3 coolSide = mix(midnightBlue, astralBlue,
        smoothstep(0.08, 0.72, phase));
    vec3 warmSide = mix(violet, magenta,
        smoothstep(0.30, 0.90, phase));
    return mix(coolSide, warmSide,
        smoothstep(0.18, 0.88, colorAxis));
}

float rift_horizon_fade(vec3 viewDir) {
    #ifdef DISTANT_RENDER_MOD
        return smoothstep(-0.15, 0.018, viewDir.y);
    #else
        // Keep the joined mouth visible at the true horizon. Fine cosmic
        // detail is faded separately, so this wider visibility does not
        // reintroduce the projection fold that used to form there.
        return smoothstep(-0.025, 0.020, viewDir.y);
    #endif
}

// Evaluate the inexpensive two-tone field shared by the visible sky and the
// terrain-fog lookup. Keeping this separate from the rift interior prevents
// stars, atlas detail and bright boundary energy from being projected onto
// partially fogged terrain at short render distances.
void rift_global_sky_fields(
        vec3 viewDir,
        vec3 baseSkyColor,
        float clock,
        out vec2 celestialUv,
        out float skyPaletteFlow,
        out float skyAltitude,
        out vec3 dualSkyColor,
        out float skyTintAlpha,
        out vec3 twoToneSky,
        out vec3 riftPaletteInfluence) {
    vec3 cosmicDirection = normalize(viewDir);
    celestialUv = rift_octahedral_sky(cosmicDirection);

    // A continuous three-dimensional axis avoids the singular direction at
    // the zenith that previously formed a large dark gradient wedge.
    vec3 paletteDirection = normalize(vec3(
        cos(clock * 0.0063),
        0.24 + sin(clock * 0.0027) * 0.14,
        sin(clock * 0.0063)));
    skyPaletteFlow = 0.5 + 0.5 * dot(
        cosmicDirection, paletteDirection);
    skyPaletteFlow += (rift_noise(
        celestialUv * 1.85
            + vec2(clock * 0.0031, -clock * 0.0024)) - 0.5)
        * 0.14;
    skyPaletteFlow = smoothstep(
        0.02, 0.98, clamp(skyPaletteFlow, 0.0, 1.0));

    skyAltitude = smoothstep(-0.02, 0.82, viewDir.y);
    vec3 sapphireSky = mix(
        vec3(0.014, 0.050, 0.140),
        vec3(0.028, 0.125, 0.300), skyAltitude);
    vec3 violetSky = mix(
        vec3(0.062, 0.018, 0.110),
        vec3(0.150, 0.034, 0.230), skyAltitude);
    dualSkyColor = mix(sapphireSky, violetSky, skyPaletteFlow);

    float skyColorBreath = 0.92
        + sin(clock * 0.027 + skyPaletteFlow * 4.3) * 0.08;
    skyTintAlpha = clamp(
        (0.24 + skyAltitude * 0.10) * skyColorBreath,
        0.20, 0.38);

    vec3 sky = max(baseSkyColor, vec3(0.0));
    twoToneSky = sky * 0.72
        + dualSkyColor * (0.72 + skyAltitude * 0.20);
    // Sub-one-code-value screen dither breaks 8-bit contouring without
    // becoming visible noise or disturbing the slow celestial motion.
    float skyDither = (rift_hash12(
        floor(gl_FragCoord.xy) + vec2(41.7, 93.1)) - 0.5)
        / 255.0;
    twoToneSky = max(twoToneSky + vec3(skyDither), vec3(0.0));

    riftPaletteInfluence = mix(
        vec3(0.91, 1.04, 1.16),
        vec3(1.13, 0.91, 1.08), skyPaletteFlow);
}

// Fog-safe part of Aurora Astral. This is deliberately limited to the broad
// night-sky tint; the detailed opening remains exclusive to the real sky.
vec4 getCosmicSkyWash(vec3 viewDirection, vec3 baseSkyColor) {
    #ifndef COSMIC_RIFT
        return vec4(0.0);
    #else
        float visibility = rift_visibility();
        if (visibility < 0.002) return vec4(0.0);

        vec3 viewDir = normalize(viewDirection);
        float clock = frameTimeCounter * COSMIC_RIFT_MOTION;
        vec2 celestialUv;
        float skyPaletteFlow;
        float skyAltitude;
        vec3 dualSkyColor;
        float skyTintAlpha;
        vec3 twoToneSky;
        vec3 riftPaletteInfluence;
        rift_global_sky_fields(
            viewDir, baseSkyColor, clock,
            celestialUv, skyPaletteFlow, skyAltitude,
            dualSkyColor, skyTintAlpha, twoToneSky,
            riftPaletteInfluence);

        return vec4(
            twoToneSky * skyTintAlpha * visibility,
            skyTintAlpha * visibility);
    #endif
}

// A cheap post-cloud energy layer used only by rainy nights. It redraws no
// galaxy body; it exposes moving fracture energy and short rim-bound knots so
// dense storm clouds can be lit from the catastrophic rift behind them.
vec3 getCosmicCataclysmFront(vec3 viewDirection) {
    #ifndef COSMIC_RIFT
        return vec3(0.0);
    #else
        float amount = rift_cataclysm_amount();
        if (amount < 0.002) return vec3(0.0);

        vec3 viewDir = normalize(viewDirection);
        float horizonFade = smoothstep(0.025, 0.18, viewDir.y);
        if (horizonFade < 0.001) return vec3(0.0);

        float clock = frameTimeCounter * COSMIC_RIFT_MOTION;
        // These clocks must remain identical to the calm rift clocks. Weather
        // changes only the amplitude of the catastrophic deformation; it
        // never changes accumulated phase or teleports the projected shape.
        float boundaryClock = clock * 0.30;
        float fractureClock = clock * 0.78;
        float textureClock = clock * 0.80;
        float astralTime = boundaryClock * 0.05;
        vec2 p = rift_astral_projection(viewDir, astralTime);

        float lowerEdge;
        float upperEdge;
        rift_twin_boundaries(
            p, astralTime, boundaryClock, lowerEdge, upperEdge);
        float pressure;
        rift_cataclysm_boundaries(
            p, textureClock, fractureClock, amount,
            lowerEdge, upperEdge, pressure);

        float bandSigned = min(p.y - lowerEdge, upperEdge - p.y);
        float aa = max(fwidth(bandSigned), 0.0008);
        float band = smoothstep(-aa * 1.5, aa * 1.5, bandSigned);
        float edgeDistance = min(abs(p.y - lowerEdge), abs(p.y - upperEdge));
        float tornRim = exp(-edgeDistance * edgeDistance
            / max(0.028 * 0.028, aa * aa * 2.5));

        vec2 celestialUv = rift_octahedral_sky(viewDir);
        vec2 stormUvA = celestialUv * vec2(1.62, 1.28)
            + vec2(textureClock * 0.024, -textureClock * 0.034);
        vec2 stormUvB = mat2(0.62, -0.78, 0.78, 0.62) * celestialUv
            * vec2(1.18, 1.74)
            + vec2(-textureClock * 0.031, textureClock * 0.020);
        vec3 stormA = texture2D(astralRiftCataclysm,
            rift_cataclysm_uv(stormUvA)).rgb;
        #if PROFILE_QUALITY == 1
            vec3 stormMap = stormA;
        #else
            vec3 stormB = texture2D(astralRiftCataclysm,
                rift_cataclysm_uv(stormUvB)).rgb;
            vec3 stormMap = mix(stormA, stormB,
                0.46 + sin(textureClock * 0.23) * 0.16);
        #endif
        float energy = dot(stormMap, vec3(0.2126, 0.7152, 0.0722));
        float veins = smoothstep(0.20, 0.67, energy);

        float pulseClock = textureClock * 0.31;
        float pulse = pow(max(sin(fract(pulseClock) * 3.14159265), 0.0), 14.0)
            * smoothstep(0.58, 0.96,
                rift_hash12(vec2(floor(pulseClock), 409.0)));
        vec3 energyColor = mix(
            vec3(0.05, 0.30, 0.82),
            vec3(0.50, 0.10, 0.70),
            clamp(stormMap.r * 1.6, 0.0, 1.0));
        vec3 front = energyColor
            * (veins * band * (0.075 + pressure * 0.055)
                + tornRim * (0.060 + veins * 0.090)
                + pulse * veins * band * 0.18);
        return front * amount * horizonFade;
    #endif
}

vec4 getCosmicRift(vec3 viewDirection, vec3 baseSkyColor) {
    #ifndef COSMIC_RIFT
        return vec4(0.0);
    #else
        float visibility = rift_visibility();
        if (visibility < 0.002) return vec4(0.0);

        vec3 viewDir = normalize(viewDirection);
        float horizonFade = rift_horizon_fade(viewDir);
        // The broad two-tone sky must remain continuous across y=0. Only the
        // dimensional opening is clipped at the true horizon; fading the sky
        // tint there creates a dark horizontal seam above distant terrain.
        if (horizonFade < 0.001)
            return getCosmicSkyWash(viewDirection, baseSkyColor);

        float clock = frameTimeCounter * COSMIC_RIFT_MOTION;
        float cataclysmAmount = rift_cataclysm_amount();
        // Geometry, fractures and the cosmic volume now have related but
        // deliberately different clocks.  The main opening crosses the sky
        // slowly, while texture motion remains alive without racing its rim.
        // Never multiply accumulated time by live rain strength: doing so
        // changes phase during a weather transition and makes the storm rift
        // jump to a seemingly fixed coordinate. Both states now share one
        // continuous path; cataclysmAmount only blends deformation and light.
        float boundaryClock = clock * 0.30;
        float fractureClock = clock * 0.78;
        float textureClock = clock * 0.80;
        float astralTime = boundaryClock * 0.05;
        vec2 p = rift_astral_projection(viewDir, astralTime);

        float lowerEdge;
        float upperEdge;
        rift_twin_boundaries(
            p, astralTime, boundaryClock, lowerEdge, upperEdge);

        // During a rift storm the two sides no longer breathe coherently.
        // Several spatial scales tear the opening asymmetrically, while a
        // slow pressure pulse makes it surge rather than merely move faster.
        float catastrophePressure;
        rift_cataclysm_boundaries(
            p, textureClock, fractureClock, cataclysmAmount,
            lowerEdge, upperEdge, catastrophePressure);

        // Both independently moving sides converge into one dimensional
        // mouth near the horizon.  Sharing only this final midpoint preserves
        // their freedom overhead while eliminating two unrelated loose ends.
        float endpointMidpoint = (lowerEdge + upperEdge) * 0.5;
        float endpointGapScale = mix(0.055, 1.0,
            smoothstep(0.010, 0.145, viewDir.y));
        lowerEdge = mix(endpointMidpoint, lowerEdge, endpointGapScale);
        upperEdge = mix(endpointMidpoint, upperEdge, endpointGapScale);
        float endpointJoin = 1.0
            - smoothstep(0.018, 0.145, viewDir.y);

        // Fracture signed fields are welded into the two polar half-spaces
        // before the opening mask is evaluated.  They grow and flex in local
        // boundary space, inheriting all motion of their parent edge.
        float fractureHorizonFade = smoothstep(
            0.045, 0.175, viewDir.y);
        // Two decorrelated skeleton lanes per pole provide an even population
        // of roots without turning the entire boundary into constant noise.
        float lowerFractureSigned = max(
            rift_integrated_fracture(
                p.x, lowerEdge - p.y, fractureClock, 13.0),
            rift_integrated_fracture(
                p.x + 1.57, lowerEdge - p.y, fractureClock, 29.0))
            - (1.0 - fractureHorizonFade) * 0.85;
        float upperFractureSigned = max(
            rift_integrated_fracture(
                p.x, p.y - upperEdge, fractureClock, 47.0),
            rift_integrated_fracture(
                p.x - 1.31, p.y - upperEdge, fractureClock, 71.0))
            - (1.0 - fractureHorizonFade) * 0.85;

        float lowerSideSigned = p.y - lowerEdge;
        float upperSideSigned = upperEdge - p.y;
        float mainBandSigned = min(
            lowerSideSigned, upperSideSigned);
        float signedInside = max(mainBandSigned,
            max(lowerFractureSigned, upperFractureSigned));

        // Distances are taken from the unioned side fields.  Consequently the
        // original polar ribbon bends around every branch and no internal
        // seam or separately drawn fracture border can survive at the root.
        float lowerDistance = abs(max(
            lowerSideSigned, lowerFractureSigned));
        float upperDistance = abs(max(
            upperSideSigned, upperFractureSigned));
        float nearestEdge = min(lowerDistance, upperDistance);
        float edgeAA = max(fwidth(signedInside), 0.00065);
        float bandMask = smoothstep(
            -edgeAA * 1.20, edgeAA * 1.20, signedInside);
        // One signed opening now contains both the main rift and its polar
        // deformations.  Cosmos, rim and glow all derive from this geometry.
        float interiorFeather = max(
            0.055 * COSMIC_RIFT_SCALE, edgeAA * 1.8);
        float interiorMask = smoothstep(
            0.0, interiorFeather, signedInside);

        // Broad aurora-like edge energy is the only component permitted just
        // outside the band.  Skip all expensive galaxy work farther away.
        float haloWidth = max(0.18 * COSMIC_RIFT_SCALE,
            edgeAA * 1.65);
        float broadHalo = exp(-nearestEdge * nearestEdge
            / max(haloWidth * haloWidth, 0.0001));
        float atmosphereWidth = max(0.66 * COSMIC_RIFT_SCALE,
            edgeAA * 4.0);
        float atmosphereWash = exp(-nearestEdge * nearestEdge
            / max(atmosphereWidth * atmosphereWidth, 0.0001));

        // Aurora Astral's color travel is reinterpreted as a two-tone night
        // field. The same lightweight evaluation is reused by terrain fog so
        // a short render distance cannot expose a differently colored band.
        vec3 cosmicDirection = normalize(viewDir);
        vec2 celestialUv;
        float skyPaletteFlow;
        float skyAltitude;
        vec3 dualSkyColor;
        float skyTintAlpha;
        vec3 twoToneSky;
        vec3 riftPaletteInfluence;
        rift_global_sky_fields(
            viewDir, baseSkyColor, clock,
            celestialUv, skyPaletteFlow, skyAltitude,
            dualSkyColor, skyTintAlpha, twoToneSky,
            riftPaletteInfluence);
        // Even pixels far from the opening receive the inexpensive two-color
        // night blend; expensive galaxy work remains restricted to the rift.
        if (bandMask < 0.001 && broadHalo < 0.002
                && atmosphereWash < 0.002)
            return vec4(
                twoToneSky * skyTintAlpha * visibility,
                skyTintAlpha * visibility);

        // The field below depends only on normalized sky direction and time.
        // Boundary coordinates never enter these equations.  The unified
        // opening mask reveals this independent sky volume.

        // Octahedral projection naturally compresses at the horizon seam.
        // Remove fine atlas, particle and noise detail before that seam while
        // leaving the two colored boundary ribbons visible.  The endpoint
        // therefore closes as a clean glow instead of a folded texture blob.
        float horizonDetailFade = smoothstep(
            0.025, 0.160, viewDir.y);

        // Three exposure-matched cosmic veils share one atlas.  None contains
        // an axial dust lane, so the interior can no longer read as a second
        // crack.  A forty-minute circular dissolve changes the celestial
        // scene too slowly for its transition point to become perceptible.
        vec2 atlasDirection = mat2(
            0.93969, -0.34202,
            0.34202,  0.93969) * celestialUv;

        // Two evolving flow fields bend local detail independently.  Their
        // displacement remains within a few percent of the texture, enough
        // to read clearly without liquefying the artwork.  Their higher
        // temporal rate makes the gas visibly alive instead of merely
        // translating a static photograph.
        float flowNoiseX = rift_fbm(
            celestialUv * 3.10
            + vec2(textureClock * 0.028, -textureClock * 0.021));
        float flowNoiseY = rift_fbm(
            celestialUv.yx * 2.85 + vec2(17.4, 8.1)
            + vec2(-textureClock * 0.024, textureClock * 0.034));
        vec2 livingFlow = (vec2(flowNoiseX, flowNoiseY) - 0.5)
            * 0.050;
        vec2 microShear = vec2(
            sin(celestialUv.y * 10.7 + textureClock * 0.19),
            cos(celestialUv.x * 9.3 - textureClock * 0.16)) * 0.012;

        // Anisotropic stretching is intentionally stronger than a camera
        // drift: clouds alternately lengthen and relax in different axes.
        float stretchX = sin(textureClock * 0.17 + flowNoiseY * 2.0) * 0.040
            + sin(textureClock * 0.061) * 0.016;
        float stretchY = cos(textureClock * 0.137 + flowNoiseX * 2.4) * 0.035
            + sin(textureClock * 0.047) * 0.014;
        vec2 atlasScale = vec2(0.80 + stretchX, 0.76 + stretchY);
        vec2 atlasUvA = atlasDirection
            * atlasScale
            + vec2(0.500 + textureClock * 0.0017,
                0.500 - textureClock * 0.0013)
            + livingFlow + microShear;
        vec2 atlasUvB = atlasDirection
            * vec2(atlasScale.x * 0.992, atlasScale.y * 0.984)
            + vec2(0.506 + textureClock * 0.0014,
                0.494 - textureClock * 0.0011)
            + vec2(-livingFlow.y, livingFlow.x) * 0.82
            + microShear.yx * 0.70;
        vec2 atlasUvC = atlasDirection
            * vec2(atlasScale.x * 1.012, atlasScale.y * 1.008)
            + vec2(0.494 + textureClock * 0.0019,
                0.506 - textureClock * 0.0015)
            + vec2(livingFlow.y, -livingFlow.x) * 0.70
            - microShear * 0.55;

        float textureCycle = fract(textureClock * 0.00042);
        vec3 textureWeights = vec3(
            rift_cycle_weight(textureCycle, 0.0),
            rift_cycle_weight(textureCycle, 1.0 / 3.0),
            rift_cycle_weight(textureCycle, 2.0 / 3.0));
        textureWeights /= max(dot(textureWeights, vec3(1.0)), 0.0001);

        vec3 atlasLayerA = texture2D(
            astralCosmosCycleAtlas,
            rift_cycle_atlas_uv(atlasUvA, 0.0)).rgb;
        vec3 atlasLayerB = texture2D(
            astralCosmosCycleAtlas,
            rift_cycle_atlas_uv(atlasUvB, 1.0)).rgb;
        vec3 atlasLayerC = texture2D(
            astralCosmosCycleAtlas,
            rift_cycle_atlas_uv(atlasUvC, 2.0)).rgb;
        vec3 generatedCosmos = pow(clamp(
            atlasLayerA * textureWeights.x
                + atlasLayerB * textureWeights.y
                + atlasLayerC * textureWeights.z,
            vec3(0.0), vec3(1.0)), vec3(1.12)) * 0.66;

        // Independent exposure, contrast and spectral breaths keep every
        // frame changing even during the very long cross-texture dissolve.
        float irregularBreath =
            sin(textureClock * 0.21 + flowNoiseX * 5.2) * 0.52
            + sin(textureClock * 0.083 + flowNoiseY * 7.1) * 0.31
            + sin(textureClock * 0.037
                + dot(cosmicDirection, vec3(1.7, -2.3, 1.1))) * 0.17;
        float textureBreath = 0.92 + irregularBreath * 0.20;
        float contrastBreath = 1.03
            + sin(textureClock * 0.137 + flowNoiseX * 3.3) * 0.12;
        float spectralBreath = 0.5 + 0.5 * sin(
            textureClock * 0.057 + flowNoiseY * 4.0
            + dot(cosmicDirection, vec3(2.1, -1.4, 1.7)));
        vec3 spectralTint = mix(
            vec3(0.92, 0.98, 1.08),
            vec3(1.07, 0.93, 1.04), spectralBreath);
        generatedCosmos = max(
            (generatedCosmos - vec3(0.045)) * contrastBreath
                + vec3(0.045),
            vec3(0.0)) * spectralTint * riftPaletteInfluence
                * textureBreath;
        generatedCosmos *= mix(0.08, 1.0, horizonDetailFade);

        // A dedicated, original catastrophe texture replaces the calm atlas
        // only as unstable energy. Counter-moving mirrored samples keep its
        // fractures alive and prevent a recognizable static image projection.
        mat2 catastropheRotation = mat2(
            0.81915, -0.57358,
            0.57358,  0.81915);
        vec2 catastropheUvA = atlasDirection * vec2(1.42, 1.16)
            + livingFlow * 1.8
            + vec2(textureClock * 0.018, -textureClock * 0.027);
        vec2 catastropheUvB = catastropheRotation * atlasDirection
            * vec2(1.08, 1.57) - livingFlow.yx * 1.35
            + vec2(-textureClock * 0.023, textureClock * 0.015);
        vec3 catastropheMapA = texture2D(astralRiftCataclysm,
            rift_cataclysm_uv(catastropheUvA)).rgb;
        #if PROFILE_QUALITY == 1
            vec3 catastropheMap = catastropheMapA;
        #else
            vec3 catastropheMapB = texture2D(astralRiftCataclysm,
                rift_cataclysm_uv(catastropheUvB)).rgb;
            vec3 catastropheMap = mix(
                catastropheMapA, catastropheMapB,
                0.38 + 0.22 * sin(textureClock * 0.19 + flowNoiseX * 4.0));
        #endif
        float catastropheLuma = dot(catastropheMap,
            vec3(0.2126, 0.7152, 0.0722));
        float catastropheVeins = smoothstep(0.16, 0.62, catastropheLuma);
        float catastropheCore = smoothstep(0.38, 0.82, catastropheLuma);
        catastropheMap = pow(max(catastropheMap, vec3(0.0)), vec3(1.28));
        catastropheMap *= mix(
            vec3(0.76, 1.08, 1.28),
            vec3(1.12, 0.76, 1.22),
            catastropheMap.r / max(catastropheMap.r + catastropheMap.b, 0.001));
        float generatedLuma = dot(generatedCosmos,
            vec3(0.2126, 0.7152, 0.0722));
        float generatedStructure = smoothstep(
            0.020, 0.20, generatedLuma);

        #if COSMIC_RIFT_QUALITY == 1
            // Balanced derives the secondary volume masks from the already
            // sampled painted cosmos and its two animated flow fields. The
            // atlas still carries all fine nebula detail, while this removes
            // five independent 3D FBM evaluations (the dominant ALU cost on
            // GTX 16-class GPUs) without removing a visible cosmic layer.
            float atlasEnergy = clamp(generatedLuma * 2.35, 0.0, 1.0);
            float atlasColorStructure = clamp(
                abs(atlasLayerA.b - atlasLayerC.r) * 1.85,
                0.0, 1.0);
            float deepGas = clamp(
                flowNoiseX * 0.48 + flowNoiseY * 0.27
                    + atlasEnergy * 0.25,
                0.0, 1.0);
            float cloudGas = clamp(
                flowNoiseY * 0.46 + generatedStructure * 0.34
                    + atlasColorStructure * 0.20,
                0.0, 1.0);
            float dustGas = clamp(
                (1.0 - atlasEnergy) * 0.45
                    + flowNoiseX * 0.31
                    + atlasColorStructure * 0.24,
                0.0, 1.0);
            float opacityFieldA = flowNoiseX;
            float opacityFieldB = clamp(
                flowNoiseY * 0.72 + atlasEnergy * 0.28,
                0.0, 1.0);
        #else
            vec3 deepPosition = cosmicDirection * 3.65
                + vec3(textureClock * 0.013, -textureClock * 0.010,
                    textureClock * 0.011);
            float deepGas = rift_fbm3(deepPosition);

            vec3 cloudPosition = cosmicDirection * 7.20
                + vec3(-textureClock * 0.020, textureClock * 0.017,
                    textureClock * 0.014)
                + vec3(deepGas - 0.5) * vec3(0.78, -0.46, 0.63);
            float cloudGas = rift_fbm3(cloudPosition);

            vec3 dustPosition = cosmicDirection * 13.4
                + vec3(textureClock * 0.031, textureClock * 0.022,
                    -textureClock * 0.027)
                + vec3(cloudGas - 0.5) * vec3(-0.52, 0.67, 0.41);
            float dustGas = rift_fbm3(dustPosition);

            // Broad opacity continents cross the celestial volume independently
            // of both polar boundaries.  Some long sections nearly dissolve into
            // the night while others remain legible; two counter-moving fields
            // continuously exchange those roles without a visible transition.
            vec3 opacityPositionA = cosmicDirection * 1.82
                + vec3(textureClock * 0.0024, -textureClock * 0.0018,
                    textureClock * 0.0020);
            vec3 opacityPositionB = cosmicDirection * 2.64
                + vec3(13.7, 29.1, 7.3)
                + vec3(-textureClock * 0.0017, textureClock * 0.0022,
                    textureClock * 0.0014);
            float opacityFieldA = rift_fbm3(opacityPositionA);
            float opacityFieldB = rift_fbm3(opacityPositionB);
        #endif
        float opacityExchange = 0.5 + 0.5 * sin(
            dot(cosmicDirection, vec3(3.1, -2.0, 2.6))
            + textureClock * 0.010
            + (opacityFieldA - opacityFieldB) * 2.2);
        float opacitySignal = opacityFieldA * 0.58
            + opacityFieldB * 0.24 + opacityExchange * 0.18;
        float opacityRegion = smoothstep(0.35, 0.66, opacitySignal);
        float regionalOpacity = mix(0.070, 1.0, opacityRegion);
        float regionalBreath = 0.90 + 0.10 * sin(
            textureClock * 0.026 + opacityFieldA * 5.1
            - opacityFieldB * 3.2);
        regionalOpacity = clamp(
            regionalOpacity * regionalBreath, 0.055, 1.0);
        float regionalEdgePresence = mix(0.16, 1.0,
            smoothstep(0.06, 0.94, regionalOpacity));

        float nebulaBody = smoothstep(0.28, 0.73,
            deepGas * 0.62 + cloudGas * 0.38);
        float luminousClouds = smoothstep(0.43, 0.80, cloudGas)
            * (0.34 + deepGas * 0.66);
        float filamentField = 1.0 - abs(cloudGas * 2.0 - 1.0);
        float luminousFilaments = pow(smoothstep(
            0.53, 0.91, filamentField), 2.0) * nebulaBody;
        // Broad shadow pockets add depth without drawing another linear seam.
        float dustPockets = smoothstep(0.54, 0.86, dustGas)
            * smoothstep(0.34, 0.76, deepGas);
        float stellarCloud = smoothstep(0.34, 0.76,
            deepGas * 0.52 + cloudGas * 0.48);

        // Sparse event envelopes produce localized, irregular cosmic glints.
        // They brighten existing gas for a moment instead of flashing the
        // entire rift or introducing a repeating global strobe.
        float flashClockA = textureClock * 0.085;
        float flashLocalA = fract(flashClockA);
        float flashEnvelopeA = pow(max(sin(
            flashLocalA * 3.14159265), 0.0), 10.0);
        float flashChanceA = smoothstep(0.55, 0.96,
            rift_hash12(vec2(floor(flashClockA), 19.7)));
        float flashClockB = textureClock * 0.052 + 7.3;
        float flashLocalB = fract(flashClockB);
        float flashEnvelopeB = pow(max(sin(
            flashLocalB * 3.14159265), 0.0), 12.0);
        float flashChanceB = smoothstep(0.62, 0.97,
            rift_hash12(vec2(floor(flashClockB), 47.1)));
        float cosmicFlash = clamp(
            flashEnvelopeA * flashChanceA
            + flashEnvelopeB * flashChanceB * 0.72,
            0.0, 1.0);
        float flashMask = smoothstep(0.40, 0.78,
            cloudGas * 0.58 + flowNoiseY * 0.42)
            * (0.36 + stellarCloud * 0.64);

        float colorAxis = clamp(0.5 + 0.5 * dot(
            cosmicDirection, normalize(vec3(0.54, 0.28, -0.79)))
            + (deepGas - 0.5) * 0.20, 0.0, 1.0);
        float colorPhase = fract(deepGas * 0.51
            + cloudGas * 0.31 + dustGas * 0.13
            + dot(cosmicDirection, vec3(0.19, -0.33, 0.27))
            - textureClock * 0.0085);
        vec3 nebulaColor = rift_nebula_palette(
            colorPhase, colorAxis);
        vec3 cloudColor = mix(
            vec3(0.035, 0.24, 0.62),
            vec3(0.56, 0.075, 0.43), colorAxis);
        vec3 filamentColor = mix(
            vec3(0.08, 0.46, 0.82),
            vec3(0.63, 0.14, 0.67),
            smoothstep(0.24, 0.86, colorAxis));

        vec3 interior = mix(
            vec3(0.0007, 0.0018, 0.0090),
            vec3(0.0040, 0.0100, 0.0380), deepGas);
        interior += max(baseSkyColor, vec3(0.0)) * 0.028;
        interior += nebulaColor * nebulaBody * 0.36;
        interior += cloudColor * luminousClouds * 0.28;
        interior += filamentColor * luminousFilaments * 0.17;
        interior += mix(
            vec3(0.045, 0.14, 0.40),
            vec3(0.22, 0.07, 0.36), colorAxis)
            * stellarCloud * 0.10;
        interior *= 1.0 - dustPockets * 0.20;
        interior *= mix(
            vec3(1.0), riftPaletteInfluence, 0.24);
        interior = interior * 0.22
            + generatedCosmos * (0.78 + stellarCloud * 0.040);
        float stormSurge = 0.70 + catastrophePressure * 0.30
            + cosmicFlash * 0.24;
        interior = mix(interior, interior * vec3(0.70, 0.82, 1.02),
            cataclysmAmount * 0.42);
        interior += catastropheMap
            * (catastropheVeins * 0.34 + catastropheCore * 0.24)
            * cataclysmAmount * stormSurge * horizonDetailFade;
        vec3 flashColor = mix(
            vec3(0.12, 0.38, 0.72),
            vec3(0.42, 0.16, 0.62), spectralBreath);
        interior += (generatedCosmos * 0.16 + flashColor * 0.060)
            * cosmicFlash * flashMask * horizonDetailFade;

        // A fourth, fully procedural veil travels independently of all three
        // images.  Its low amplitude supplies living depth rather than a new
        // visible shape or another photographic layer.
        float movingVeil = smoothstep(0.42, 0.77,
            flowNoiseX * 0.56 + flowNoiseY * 0.44);
        vec3 movingVeilColor = mix(
            vec3(0.025, 0.12, 0.30),
            vec3(0.17, 0.045, 0.26), flowNoiseY);
        interior += movingVeilColor * movingVeil
            * (0.026 + luminousClouds * 0.022)
            * horizonDetailFade;

        // Broad counter-traveling light pulses pass through the volume.  They
        // change illumination rather than drawing sharp new filaments.
        float travelingLightA = 0.5 + 0.5 * sin(
            dot(cosmicDirection, vec3(7.2, -3.8, 5.1))
            + textureClock * 0.19 + flowNoiseX * 3.2);
        float travelingLightB = 0.5 + 0.5 * sin(
            dot(cosmicDirection, vec3(-4.6, 6.7, 3.3))
            - textureClock * 0.13 + flowNoiseY * 3.7);
        float livingLight = smoothstep(0.18, 0.88,
            travelingLightA * 0.58 + travelingLightB * 0.42);
        vec3 livingLightColor = mix(
            vec3(0.035, 0.15, 0.34),
            vec3(0.21, 0.055, 0.29), spectralBreath);
        interior += livingLightColor * livingLight
            * (0.018 + nebulaBody * 0.030)
            * horizonDetailFade;

        float volumeBreath = 0.90
            + sin(textureClock * 0.143 + deepGas * 5.2) * 0.10
            + sin(textureClock * 0.053 + cloudGas * 6.1) * 0.055;
        interior *= volumeBreath;

        // Close to the projection seam, replace all fine structure with one
        // dim breathing color field.  This is deliberately textureless.
        vec3 quietHorizon = mix(
            vec3(0.006, 0.022, 0.070),
            vec3(0.030, 0.010, 0.060), colorAxis);
        float quietHorizonBreath = 0.84
            + sin(textureClock * 0.052 + colorAxis * 4.0) * 0.08
            + sin(textureClock * 0.113 + flowNoiseY * 5.0) * 0.04;
        quietHorizon *= quietHorizonBreath;
        quietHorizon += max(baseSkyColor, vec3(0.0)) * 0.018;
        quietHorizon = mix(
            quietHorizon, dualSkyColor * 0.58, 0.30);
        interior = mix(quietHorizon, interior, horizonDetailFade);

        // Octahedral sky cells create resolution-independent particle depth.
        // Their coordinates also ignore the opening, so moving boundaries
        // reveal and hide stars rather than dragging or stretching them.
        vec2 farStarUv = celestialUv * vec2(168.0, 154.0)
            + vec2(textureClock * 0.018, -textureClock * 0.014);
        vec2 nearStarUv = celestialUv * vec2(286.0, 264.0)
            + vec2(-textureClock * 0.031, textureClock * 0.024);
        float farStars = rift_star_layer(
            farStarUv, 0.865, textureClock, 7.0);
        float nearStars = rift_star_layer(
            nearStarUv, 0.948, textureClock, 23.0);
        float starOcclusion = 1.0 - dustPockets * 0.32;
        farStars *= (0.32 + stellarCloud * 0.92) * starOcclusion;
        nearStars *= (0.58 + stellarCloud * 0.72) * starOcclusion;

        float farHue = rift_hash12(floor(farStarUv) + 5.3);
        float nearHue = rift_hash12(floor(nearStarUv) + 17.7);
        vec3 farStarColor = mix(
            vec3(0.20, 0.44, 0.92),
            vec3(0.66, 0.38, 0.88), farHue);
        vec3 nearStarColor = mix(
            vec3(0.52, 0.76, 1.0),
            vec3(0.96, 0.58, 0.83), nearHue);
        vec3 starRadiance = farStarColor * farStars * 0.12;
        starRadiance += nearStarColor * nearStars * 0.25;
        starRadiance *= horizonDetailFade;

        // Colored, broken aurora boundaries.  Geometry, curtain texture and
        // traveling energy packets all evolve independently for each side.
        float lowerTexture = rift_fbm(vec2(
            p.x * 0.72 - boundaryClock * 0.048,
            17.0 + boundaryClock * 0.019));
        float upperTexture = rift_fbm(vec2(
            p.x * 0.68 + boundaryClock * 0.043,
            53.0 - boundaryClock * 0.017));
        float lowerPackets = 0.5 + 0.5 * sin(
            p.x * 3.7 - boundaryClock * 0.20 + lowerTexture * 5.0);
        float upperPackets = 0.5 + 0.5 * sin(
            p.x * 3.3 + boundaryClock * 0.17 + upperTexture * 5.0);
        lowerPackets = 0.34 + pow(max(lowerPackets, 0.0), 3.0) * 0.66;
        upperPackets = 0.34 + pow(max(upperPackets, 0.0), 3.0) * 0.66;

        // Endpoint ribbons keep their color but lose noisy texture and packet
        // breakup as they converge at the horizon.
        float endpointDetail = mix(0.18, 1.0, horizonDetailFade);
        lowerTexture = mix(0.50, lowerTexture, endpointDetail);
        upperTexture = mix(0.50, upperTexture, endpointDetail);
        lowerPackets = mix(0.55, lowerPackets, endpointDetail);
        upperPackets = mix(0.55, upperPackets, endpointDetail);

        // The boundary is a layered spectral curtain rather than a laser:
        // one smooth color band dissolving into a wider atmospheric aura.
        float ribbonWidth = max(0.060 * COSMIC_RIFT_SCALE,
            edgeAA * 1.15);
        float lowerRibbon = exp(-lowerDistance * lowerDistance
            / max(ribbonWidth * ribbonWidth, 0.00002));
        float upperRibbon = exp(-upperDistance * upperDistance
            / max(ribbonWidth * ribbonWidth, 0.00002));
        lowerRibbon *= (0.38 + lowerTexture * 0.62) * lowerPackets;
        upperRibbon *= (0.38 + upperTexture * 0.62) * upperPackets;
        lowerRibbon *= mix(1.0,
            0.62 + catastrophePressure * 1.15 + catastropheVeins * 0.48,
            cataclysmAmount);
        upperRibbon *= mix(1.0,
            1.74 - catastrophePressure * 1.04 + catastropheCore * 0.52,
            cataclysmAmount);

        float lowerAura = exp(-lowerDistance * lowerDistance
            / max(haloWidth * haloWidth, 0.0001));
        float upperAura = exp(-upperDistance * upperDistance
            / max(haloWidth * haloWidth, 0.0001));
        lowerAura *= 0.30 + lowerTexture * 0.70;
        upperAura *= 0.30 + upperTexture * 0.70;
        // The same regional field attenuates both the opening and its welded
        // rim, preventing a permanently visible outline around faded areas.
        lowerRibbon *= regionalEdgePresence;
        upperRibbon *= regionalEdgePresence;
        lowerAura *= regionalEdgePresence;
        upperAura *= regionalEdgePresence;

        vec3 lowerColor = mix(
            vec3(0.020, 0.105, 0.34),
            vec3(0.090, 0.36, 0.66), lowerTexture);
        vec3 upperColor = mix(
            vec3(0.105, 0.030, 0.32),
            vec3(0.36, 0.090, 0.55), upperTexture);
        lowerColor *= mix(
            vec3(1.0), riftPaletteInfluence, 0.18);
        upperColor *= mix(
            vec3(1.0), riftPaletteInfluence, 0.18);

        float endpointBridgeWidth = max(
            0.14 * COSMIC_RIFT_SCALE, edgeAA * 2.1);
        float endpointDistance = abs(p.y - endpointMidpoint);
        float endpointBridge = exp(
            -endpointDistance * endpointDistance
                / max(endpointBridgeWidth * endpointBridgeWidth, 0.0001))
            * endpointJoin * regionalEdgePresence;
        vec3 endpointColor = mix(
            lowerColor, upperColor, 0.5);

        // Every dimensional contribution uses the unified opening mask, so
        // polar fractures share the exact opacity and life of the parent rift.
        float opacityLife = clamp(0.90
            + sin(textureClock * 0.173 + deepGas * 3.7) * 0.12
            + sin(textureClock * 0.071 + cloudGas * 4.9) * 0.065
            + sin(textureClock * 0.311 + flowNoiseX * 5.7) * 0.035,
            0.70, 1.08);
        float endpointOpacity = mix(0.72, 1.0, horizonDetailFade);
        float interiorAlpha = interiorMask
            * (0.50 + nebulaBody * 0.030
                + luminousFilaments * 0.012
                + generatedStructure * 0.075)
            * opacityLife * endpointOpacity * regionalOpacity;
        float edgeVeil = clamp(
            (lowerAura + upperAura) * 0.021
            + (lowerRibbon + upperRibbon) * 0.030
            + endpointBridge * 0.050,
            0.0, 0.115);

        // A broad, low-opacity sky wash shares the rift palette and altitude
        // gradient.  It changes the surrounding sky rather than drawing a
        // colored rectangle, making the cosmic window belong to the night.
        float atmospherePulse = 0.72
            + sin(textureClock * 0.041 + flowNoiseY * 3.8) * 0.18
            + sin(textureClock * 0.017 + deepGas * 4.2) * 0.10;
        float atmosphereAlpha = clamp(atmosphereWash
            * (0.025 + atmospherePulse * 0.025)
            * mix(0.32, 1.0, regionalOpacity), 0.0, 0.060);
        // The local veil is always at least as luminous as the global sky;
        // it can no longer carve a dark band around the moving rift.
        vec3 atmosphericColor = twoToneSky
            + dualSkyColor * (0.055 + movingVeil * 0.035);

        // The weak boundary veil carries the live sky through unchanged, so
        // its outer falloff can add colored atmosphere without creating dark
        // smoke around the band.
        vec3 radiance = twoToneSky * skyTintAlpha;
        radiance = atmosphericColor * atmosphereAlpha
            + radiance * (1.0 - atmosphereAlpha);
        radiance = twoToneSky * edgeVeil
            + radiance * (1.0 - edgeVeil);
        radiance = interior * interiorAlpha
            + radiance * (1.0 - interiorAlpha);
        radiance += starRadiance * interiorMask * regionalOpacity;
        radiance += lowerColor
            * (lowerAura * 0.065 + lowerRibbon * 0.108);
        radiance += upperColor
            * (upperAura * 0.065 + upperRibbon * 0.108);
        radiance += endpointColor * endpointBridge * 0.16;

        float totalAlpha = 1.0
            - (1.0 - interiorAlpha)
                * (1.0 - atmosphereAlpha) * (1.0 - edgeVeil)
                * (1.0 - skyTintAlpha);
        // Preserve the global sky layer through the horizon while fading only
        // the rift-specific interior, rim and atmospheric energy. This avoids
        // multiplying the entire sky by the narrow local horizon mask.
        vec3 globalSkyRadiance = twoToneSky * skyTintAlpha;
        radiance = mix(
            globalSkyRadiance,
            radiance * COSMIC_RIFT_INTENSITY,
            horizonFade);
        totalAlpha = mix(skyTintAlpha, totalAlpha, horizonFade);
        radiance *= visibility;
        totalAlpha *= visibility;

        return vec4(max(radiance, vec3(0.0)),
            clamp(totalAlpha, 0.0, 0.88));
    #endif
}

#endif // ASTRAL_COSMIC_RIFT_GLSL
