/*
 * =======================================================================================
 *
 * CELESTIAL EFFECTS MODULE
 * Lunara Shader - Sky & Space Visual Effects
 *
 * =======================================================================================
 *
 * Copyright (c) 2025 Lunara Development
 * All rights reserved. Unauthorized reproduction or distribution is prohibited.
 *
 * This consolidated module handles all celestial and space-based atmospheric effects
 * including aurora borealis, night nebula, stars systems, and dimensional star effects.
 * Optimized for performance with quality-adaptive rendering.
 *
 * =======================================================================================
 */

#ifndef LU_INCLUDE_CELESTIAL_EFFECTS
#define LU_INCLUDE_CELESTIAL_EFFECTS

#include "/lib/colors/core.glsl"

//=======================================================================================
// STARS SYSTEM MODULE
//=======================================================================================

/**
 * Generates pseudo-random noise for star placement
 * @param pos: 2D position coordinate
 * @return: Random value [0.0, 1.0] for star generation
 */
float lu_GetStarNoise(vec2 pos) {
    return fract(sin(dot(pos, vec2(12.9898, 4.1414))) * 43758.54953);
}

/**
 * Calculates star coordinate mapping from view position
 * @param viewPos: View space position vector
 * @param sphereness: Spherical projection factor [0.0, 1.0]
 * @return: 2D star coordinate for texture sampling
 */
vec2 lu_GetStarCoord(vec3 viewPos, float sphereness) {
    vec3 wpos = normalize((gbufferModelViewInverse * vec4(viewPos * 1000.0, 1.0)).xyz);
    vec3 starCoord = wpos / (wpos.y + length(wpos.xz) * sphereness);
    starCoord.x += 0.006 * syncedTime;  // Subtle star movement
    return starCoord.xz;
}

/**
 * Generates overworld star field with configurable density
 * @param starCoord: 2D star coordinate from GetStarCoord
 * @param VdotU: View direction dot up vector
 * @param VdotS: View direction dot sun vector
 * @return: RGB star color contribution
 */
vec3 lu_GetStars(vec2 starCoord, float VdotU, float VdotS) {
    if (VdotU < 0.0) return vec3(0.0);  // No stars below horizon

    starCoord *= 0.2;
    float starFactor = 1024.0;
    vec2 rawCoord = starCoord; // continuous coord for galaxy/meteor
    starCoord = floor(starCoord * starFactor) / starFactor;

    // Multi-layer star noise for realistic distribution
    float star = 1.0;
    float starSeed = lu_GetStarNoise(starCoord.xy);
    star *= starSeed;
    star *= lu_GetStarNoise(starCoord.xy + 0.1);
    star *= lu_GetStarNoise(starCoord.xy + 0.23);

    #if LU_NIGHT_STAR_AMOUNT == 2
        star -= 0.7;  // High star density
    #else
        star -= 0.6;  // Standard star density
        star *= 0.65;
    #endif
    star = lu_max0(star);
    star *= star;  // Square for more dramatic falloff

    // === ADVANCED NIGHT SKY: per-star twinkle ===
    #ifdef LU_STAR_TWINKLE
        float twinklePhase = starSeed * 6.2831 + syncedTime * (0.8 + starSeed * 1.6);
        float twinkle = 0.55 + 0.45 * sin(twinklePhase);
        star *= twinkle;
    #endif

    // Apply atmospheric and lighting conditions
    star *= lu_min1(VdotU * 3.0) * lu_max0(1.0 - pow(abs(VdotS) * 1.002, 100.0));
    star *= invRainFactor * lu_pow2(lu_pow2(invNoonFactor2)) * (1.0 - 0.5 * sunVisibility);

    // === ADVANCED NIGHT SKY: per-star stellar color classes ===
    vec3 starColor = vec3(0.38, 0.4, 0.5);
    #ifdef LU_STAR_COLORS
        float hueSeed = lu_GetStarNoise(starCoord.xy + 0.51);
        vec3 blueStar   = vec3(0.30, 0.42, 0.65);
        vec3 whiteStar  = vec3(0.50, 0.50, 0.55);
        vec3 orangeStar = vec3(0.65, 0.42, 0.28);
        vec3 stellar = mix(blueStar, whiteStar, lu_clamp01(hueSeed * 3.0));
        stellar = mix(stellar, orangeStar, lu_clamp01(hueSeed * 3.0 - 2.0));
        starColor = stellar;
    #endif

    vec3 result = 40.0 * star * starColor;

    // Common sky-visibility factor for galaxy & meteors
    float skyVis = invRainFactor * lu_pow2(lu_pow2(invNoonFactor2)) * (1.0 - 0.5 * sunVisibility);

    // === ADVANCED NIGHT SKY: Milky Way galaxy band ===
    #ifdef LU_GALAXY_BAND
        vec2 gCoord = rawCoord;
        float gAngle = 0.55;
        float cs = cos(gAngle), sn = sin(gAngle);
        vec2 gRot = vec2(cs * gCoord.x - sn * gCoord.y, sn * gCoord.x + cs * gCoord.y);
        float bandDist = abs(gRot.y);
        float band = exp(-bandDist * bandDist * 18.0);
        float gNoise = lu_GetStarNoise(floor(rawCoord * 64.0) / 64.0) * 0.6
                     + 0.4 * lu_GetStarNoise(floor(rawCoord * 16.0) / 16.0);
        gNoise = lu_smoothstep1(gNoise);
        float along = 0.7 + 0.3 * sin(gRot.x * 2.4);
        vec3 galaxyTint = mix(vec3(0.18, 0.10, 0.32), vec3(0.10, 0.18, 0.42), gNoise);
        float galaxyVis = lu_min1(VdotU * 2.5) * skyVis;
        result += band * gNoise * along * galaxyTint * (float(LU_GALAXY_I) * 0.01) * galaxyVis;
    #endif

    // === ADVANCED NIGHT SKY: shooting stars / meteors ===
    #ifdef LU_SHOOTING_STARS
        const int meteorCount = 2;
        for (int m = 0; m < meteorCount; m++) {
            float mSeed = float(m) * 7.137;
            float gap = 16.0 / float(LU_SHOOTING_STARS_RATE);
            float cycle = floor(syncedTime / gap + mSeed);
            float lifeT = fract(syncedTime / gap + mSeed);
            vec2 spawn = vec2(
                lu_GetStarNoise(vec2(cycle, mSeed)) * 4.0 - 2.0,
                lu_GetStarNoise(vec2(cycle + 0.7, mSeed)) * 2.0 + 0.3
            );
            float dirA = lu_GetStarNoise(vec2(cycle * 1.3, mSeed)) * 6.2831;
            vec2 dir = vec2(cos(dirA), sin(dirA));
            float visT = lu_clamp01(lifeT * 4.0);
            float fade = 1.0 - visT;
            vec2 head = spawn + dir * lifeT * 1.5;
            vec2 toPix = rawCoord - head;
            float perp = abs(dot(toPix, vec2(-dir.y, dir.x)));
            float along2 = dot(toPix, -dir);
            float trail = lu_clamp01(along2 / 0.35);
            float core  = exp(-perp * perp * 8000.0);
            float tail  = exp(-perp * perp * 1600.0) * trail;
            float meteor = (core + 0.6 * tail) * fade;
            meteor *= step(lifeT, 0.25) * lu_min1(VdotU * 3.0)
                    * invRainFactor * lu_pow2(invNoonFactor2) * (1.0 - 0.7 * sunVisibility);
            result += meteor * vec3(1.2, 1.0, 0.7) * 6.0;
        }
    #endif

    return result;
}

//=======================================================================================
// GIANT PLANET MODULE
//=======================================================================================

/**
 * Smooth 2D value noise for planet surface texture
 */
float lu_planetNoise(vec2 p) {
    vec2 i = floor(p);
    vec2 f = fract(p);
    f = f * f * (3.0 - 2.0 * f);
    float a = lu_GetStarNoise(i);
    float b = lu_GetStarNoise(i + vec2(1.0, 0.0));
    float c = lu_GetStarNoise(i + vec2(0.0, 1.0));
    float d = lu_GetStarNoise(i + vec2(1.0, 1.0));
    return mix(mix(a, b, f.x), mix(c, d, f.x), f.y);
}

float lu_planetFbm(vec2 p) {
    float v = 0.0;
    float a = 0.5;
    for (int i = 0; i < 4; i++) {
        v += a * lu_planetNoise(p);
        p *= 2.07;
        a *= 0.5;
    }
    return v;
}

/**
 * Renders a giant procedural planet hanging in the night sky.
 * The planet has:
 *  - sphere shading via fake normal from screen-space disk
 *  - terminator (day/night) lit by the moon direction
 *  - thin atmospheric rim glow
 *  - style-dependent surface (gas bands / craters / rings)
 *  - sky-visibility gating identical to stars (hidden in rain & day)
 *
 * @param starCoord: continuous sky-direction coords (same as lu_GetStars input)
 * @param VdotU:     view dot up
 * @param VdotS:     view dot sun (used for daylight fade only)
 * @return RGB planet contribution
 */
vec3 lu_GetGiantPlanet(vec2 starCoord, float VdotU, float VdotS) {
    if (VdotU < 0.05) return vec3(0.0);

    // Fixed position in the sky (slowly rotates with star coord)
    vec2 planetCenter = vec2(0.85, 0.55);
    float radius = float(LU_PLANET_SIZE) * 0.012;

    vec2 d = starCoord - planetCenter;
    float r = length(d);
    if (r > radius * 2.2) return vec3(0.0); // early out

    // Disk mask + soft anti-aliased edge
    float edgeAA = radius * 0.04;
    float diskMask = lu_clamp01((radius - r) / edgeAA);

    // Atmospheric rim glow (just outside disk)
    float rim = exp(-(r - radius) * (r - radius) / (radius * radius * 0.06));
    rim *= lu_clamp01((r - radius * 0.95) / (radius * 0.6));

    // Fake sphere normal from 2D disk (z = sqrt(R^2 - r^2))
    vec2 nUV = d / max(radius, 1e-4);
    float nz = sqrt(max(1.0 - dot(nUV, nUV), 0.0));
    vec3 sphereN = vec3(nUV, nz);

    // Light direction = opposite to the sun (it's night-time), in screen-projected sense
    // Use a stylized fixed direction so the lit limb feels natural.
    vec3 lightDir = normalize(vec3(-0.55, 0.35, 0.75));
    float NdotL = lu_clamp01(dot(sphereN, lightDir));
    float terminator = pow(NdotL, 0.85);

    // Surface UV in spherical-ish projection
    vec2 surfUV = vec2(
        atan(sphereN.x, sphereN.z) * 0.5,
        asin(sphereN.y * 0.999)
    );

    vec3 surface;

    #if LU_PLANET_STYLE == 1
        // Banded gas giant (Jupiter-like) -------------------------------
        float bands = sin(surfUV.y * 14.0 + lu_planetFbm(surfUV * 4.0) * 2.5);
        float swirl = lu_planetFbm(surfUV * 6.0 + vec2(syncedTime * 0.005, 0.0));
        bands = bands * 0.5 + 0.5;
        vec3 cBand1 = vec3(0.85, 0.62, 0.40);
        vec3 cBand2 = vec3(0.55, 0.32, 0.20);
        vec3 cBand3 = vec3(0.95, 0.80, 0.55);
        surface = mix(cBand2, cBand1, bands);
        surface = mix(surface, cBand3, lu_clamp01(swirl * 1.2 - 0.3) * 0.6);
    #elif LU_PLANET_STYLE == 2
        // Rocky cracked moon ---------------------------------------------
        float craters = lu_planetFbm(surfUV * 18.0);
        craters = pow(craters, 1.6);
        vec3 cRock1 = vec3(0.55, 0.45, 0.40);
        vec3 cRock2 = vec3(0.25, 0.18, 0.15);
        surface = mix(cRock2, cRock1, craters);
        // dark cracks
        float cracks = smoothstep(0.45, 0.5, lu_planetFbm(surfUV * 8.0));
        surface *= 1.0 - 0.6 * cracks;
    #else
        // Ringed world (Saturn-like, default) ----------------------------
        float bands = sin(surfUV.y * 11.0 + lu_planetFbm(surfUV * 3.0) * 1.8);
        bands = bands * 0.5 + 0.5;
        vec3 cBand1 = vec3(0.92, 0.82, 0.55);
        vec3 cBand2 = vec3(0.70, 0.55, 0.30);
        surface = mix(cBand2, cBand1, bands);
    #endif

    // Apply lighting: lit side full, dark side dim blue ambient
    vec3 nightAmbient = vec3(0.05, 0.07, 0.12);
    vec3 lit = surface * (0.15 + 0.85 * terminator);
    vec3 disk = mix(nightAmbient, lit, terminator * 0.85 + 0.15);

    // Soft limb darkening
    disk *= 0.55 + 0.45 * nz;

    // Atmospheric rim color (cool blue glow on lit side, fades on dark side)
    vec3 rimColor = vec3(0.30, 0.55, 0.95) * rim * (0.25 + 0.75 * NdotL);

    vec3 result = disk * diskMask + rimColor;

    // Saturn-style rings ------------------------------------------------
    #if LU_PLANET_STYLE == 3
        // Ring plane is tilted; compute distance to ring ellipse in disk space
        vec2 ringUV = d;
        // tilt: scale Y to flatten ring into ellipse
        ringUV.y /= 0.32;
        float ringR = length(ringUV);
        float ringInner = radius * 1.25;
        float ringOuter = radius * 2.05;
        if (ringR > ringInner && ringR < ringOuter && r > radius * 0.95) {
            float t = (ringR - ringInner) / (ringOuter - ringInner);
            // banding inside ring
            float ringBands = 0.4 + 0.6 * lu_planetFbm(vec2(t * 12.0, 0.0));
            // Cassini gap
            float gap = smoothstep(0.40, 0.45, t) * (1.0 - smoothstep(0.50, 0.55, t));
            ringBands *= 1.0 - 0.85 * gap;
            // edge fade
            float ringMask = smoothstep(0.0, 0.1, t) * (1.0 - smoothstep(0.85, 1.0, t));
            // rings only visible in front of (not behind) planet's far side
            // approximation: rings dim when behind the lit hemisphere
            float ringLit = 0.4 + 0.6 * lu_clamp01(dot(normalize(vec3(ringUV.x, 0.0, sqrt(max(1.0 - ringUV.x * ringUV.x / (ringOuter * ringOuter), 0.0)))), lightDir));
            vec3 ringColor = vec3(0.85, 0.78, 0.60) * ringBands * ringMask * ringLit;
            // Ring should not overdraw the planet disk
            ringColor *= 1.0 - diskMask;
            result += ringColor;
        }
    #endif

    // Sky visibility (same gating as stars)
    float skyVis = lu_min1(VdotU * 3.0)
                 * invRainFactor
                 * lu_pow2(lu_pow2(invNoonFactor2))
                 * (1.0 - 0.85 * sunVisibility);
    result *= skyVis * (float(LU_PLANET_INTENSITY) * 0.01);

    return result;
}

//=======================================================================================
// END DIMENSION STARS MODULE
//=======================================================================================

/**
 * Generates End dimension stars with magical coloring
 * @param viewPos: View space position vector
 * @param VdotU: View direction dot up vector
 * @return: RGB End star color with magical enhancement
 */
vec3 lu_GetEnderStars(vec3 viewPos, float VdotU) {
    vec3 wpos = normalize((gbufferModelViewInverse * vec4(viewPos * 1000.0, 1.0)).xyz);

    vec3 starCoord = 0.65 * wpos / (abs(wpos.y) + length(wpos.xz));
    vec2 starCoord2 = starCoord.xz * 0.5;
    if (VdotU < 0.0) starCoord2 += 100.0;
    float starFactor = 1024.0;
    starCoord2 = floor(starCoord2 * starFactor) / starFactor;

    // Triple-layer End star generation
    float star = 1.0;
    star *= lu_GetStarNoise(starCoord2.xy);
    star *= lu_GetStarNoise(starCoord2.xy + 0.1);
    star *= lu_GetStarNoise(starCoord2.xy + 0.23);
    star = max(star - 0.6, 0.0);  // More visible End stars
    star *= star;

    // magical End star colors
    vec3 baseStarColor = endSkyColor * mix(
        vec3(3500.0, 2800.0, 4500.0),
        vec3(4200.0, 3500.0, 2200.0),
        lu_GetStarNoise(starCoord2.xy + 0.5)
    );
    vec3 enderStars = star * baseStarColor;

    // Apply directional visibility factors
    float VdotUM1 = abs(VdotU);
    float VdotUM2 = lu_pow2(1.0 - VdotUM1);
    enderStars *= VdotUM1 * VdotUM1 * (VdotUM2 + 0.015) + 0.015;

    return enderStars;
}

//=======================================================================================
// AURORA BOREALIS MODULE
//=======================================================================================

/**
 * Generates aurora borealis effects with configurable conditions
 * @param viewPos: View space position vector
 * @param VdotU: View direction dot up vector
 * @param dither: Temporal dithering value
 * @return: RGB aurora color contribution
 */
vec3 lu_GetAuroraBorealis(vec3 viewPos, float VdotU, float dither) {
    float visibility = lu_sqrt1(lu_clamp01(VdotU * 1.5 - 0.225)) - sunVisibility - rainFactor - maxBlindnessDarkness;
    visibility *= 1.0 - VdotU * 0.9;

    // Apply aurora visibility conditions
    #if LU_AURORA_CONDITION == 1 || LU_AURORA_CONDITION == 3
        visibility -= moonPhase;
    #endif
    #if LU_AURORA_CONDITION == 2 || LU_AURORA_CONDITION == 3
        visibility *= inSnowy;
    #endif
    #if LU_AURORA_CONDITION == 4
        visibility = max(visibility * inSnowy, visibility - moonPhase);
    #endif

    if (visibility > 0.0) {
        vec3 aurora = vec3(0.0);

        vec3 wpos = mat3(gbufferModelViewInverse) * viewPos;
             wpos.xz /= wpos.y;
        vec2 cameraPositionM = cameraPosition.xz * 0.0075;
             cameraPositionM.x += syncedTime * 0.04;

        // Quality-adaptive sample count
        #ifdef LU_DEFERRED1
            int sampleCount = 25;
            int sampleCountP = sampleCount + 5;
        #else
            int sampleCount = 10;
            int sampleCountP = sampleCount + 10;
        #endif

        float ditherM = dither + 5.0;
        float auroraAnimate = frameTimeCounter * 0.001;

        // Multi-layer aurora generation
        for (int i = 0; i < sampleCount; i++) {
            float current = lu_pow2((i + ditherM) / sampleCountP);

            vec2 planePos = wpos.xz * (0.8 + current) * 11.0 + cameraPositionM;

            #if LU_AURORA_STYLE == 1
                // Pixelated aurora style
                planePos = floor(planePos) * 0.0007;
                float noise = texture2D(noisetex, planePos).b;
                noise = lu_pow2(lu_pow2(lu_pow2(lu_pow2(1.0 - 2.0 * abs(noise - 0.5)))));
                noise *= lu_pow1_5(texture2D(noisetex, planePos * 100.0 + auroraAnimate).b);
            #else
                // Smooth aurora style
                planePos *= 0.0007;
                float noise = texture2D(noisetex, planePos).r;
                noise = lu_pow2(lu_pow2(lu_pow2(lu_pow2(1.0 - 2.0 * abs(noise - 0.5)))));
                noise *= texture2D(noisetex, planePos * 3.0 + auroraAnimate).b;
                noise *= texture2D(noisetex, planePos * 5.0 - auroraAnimate).b;
            #endif

            float currentM = 1.0 - current;
            aurora += noise * currentM * mix(
                vec3(7.0, 2.2, 12.0),
                vec3(6.0, 16.0, 12.0),
                lu_pow2(lu_pow2(currentM))
            );
        }

        #if LU_AURORA_STYLE == 1
            aurora *= 1.3;
        #else
            aurora *= 1.8;
        #endif

        #ifdef ATM_COLOR_MULTS
            aurora *= sqrtAtmColorMult;
        #endif

        return aurora * visibility / sampleCount;
    }

    return vec3(0.0);
}

//=======================================================================================
// NIGHT NEBULA MODULE
//=======================================================================================

#ifndef LU_HQ_NIGHT_NEBULA
    const int OCTAVE = 5;
#else
    const int OCTAVE = 8;
#endif
const float timescale = 5.0;
const float zoomScale = 3.5;
const vec4 CLOUD1_COL = vec4(0.41, 0.64, 0.97, 0.4);
const vec4 CLOUD2_COL = vec4(0.81, 0.55, 0.21, 0.2);
const vec4 CLOUD3_COL = vec4(0.51, 0.81, 0.98, 1.0);

float lu_sinM(float x) {
    return sin(mod(x, 2.0 * pi));
}

float lu_cosM(float x) {
    return cos(mod(x, 2.0 * pi));
}

float lu_rand(vec2 inCoord){
    return fract(lu_sinM(dot(inCoord, vec2(23.53, 44.0))) * 42350.45);
}

float lu_perlin(vec2 inCoord){
    vec2 i = floor(inCoord);
    vec2 j = fract(inCoord);
    vec2 coord = smoothstep(0.0, 1.0, j);

    float a = lu_rand(i);
    float b = lu_rand(i + vec2(1.0, 0.0));
    float c = lu_rand(i + vec2(0.0, 1.0));
    float d = lu_rand(i + vec2(1.0, 1.0));

    return mix(mix(a, b, coord.x), mix(c, d, coord.x), coord.y);
}

float lu_fbmCloud(vec2 inCoord, float minimum){
    float value = 0.0;
    float scale = 0.5;

    for (int i = 0; i < OCTAVE; i++){
        value += lu_perlin(inCoord) * scale;
        inCoord *= 2.0;
        scale *= 0.5;
    }

    return smoothstep(0.0, 1.0, (smoothstep(minimum, 1.0, value) - minimum) / (1.0 - minimum));
}

float lu_fbmCloud2(vec2 inCoord, float minimum){
    float value = 0.0;
    float scale = 0.5;

    for (int i = 0; i < OCTAVE; i++){
        value += lu_perlin(inCoord) * scale;
        inCoord *= 2.0;
        scale *= 0.5;
    }

    return (smoothstep(minimum, 1.0, value) - minimum) / (1.0 - minimum);
}

/**
 * Generates night nebula effects with fractal noise
 * @param viewPos: View space position vector
 * @param VdotU: View direction dot up vector
 * @param VdotS: View direction dot sun vector
 * @return: RGB nebula color contribution
 */
vec3 lu_GetNightNebula(vec3 viewPos, float VdotU, float VdotS) {
    float nebulaFactor = lu_pow2(lu_max0(VdotU) * lu_min1(nightFactor * 2.0)) * invRainFactor - maxBlindnessDarkness;
    if (nebulaFactor < 0.001) return vec3(0.0);

    vec2 UV = lu_GetStarCoord(viewPos, 0.75);
    float TIME = syncedTime * 0.003 + 15.0;

    float timescaled = TIME * timescale;
    vec2 zoomUV2 = vec2(
        zoomScale * UV.x + 0.03  * timescaled * lu_sinM(0.07 * timescaled),
        zoomScale * UV.y + 0.03  * timescaled * lu_cosM(0.06 * timescaled)
    );
    vec2 zoomUV3 = vec2(
        zoomScale * UV.x + 0.027 * timescaled * lu_sinM(0.07 * timescaled),
        zoomScale * UV.y + 0.025 * timescaled * lu_cosM(0.06 * timescaled)
    );
    vec2 zoomUV4 = vec2(
        zoomScale * UV.x + 0.021 * timescaled * lu_sinM(0.07 * timescaled),
        zoomScale * UV.y + 0.021 * timescaled * lu_cosM(0.07 * timescaled)
    );

    float tide = 0.05 * lu_sinM(TIME);
    float tide2 = 0.06 * lu_cosM(0.3 * TIME);

    // Multi-layer nebula texture generation
    vec4 nebulaTexture = vec4(vec3(0.0), 0.5 + 0.2 * lu_sinM(0.23 * TIME + UV.x - UV.y));
    nebulaTexture += lu_fbmCloud2(zoomUV3, 0.24 + tide) * CLOUD1_COL;
    nebulaTexture += lu_fbmCloud(zoomUV2 * 0.9, 0.33 - tide) * CLOUD2_COL;
    nebulaTexture = mix(nebulaTexture, CLOUD3_COL, lu_fbmCloud(vec2(0.9 * zoomUV4.x, 0.9 * zoomUV4.y), 0.25 + tide2));

    nebulaFactor *= 1.0 - lu_pow2(lu_pow2(lu_pow2(abs(VdotS))));
    nebulaTexture.a *= lu_min1(lu_pow2(lu_pow2(nebulaTexture.a))) * nebulaFactor;

    // Add star integration
    float starFactor = 1024.0;
    vec2 starCoord = floor(UV * 0.25 * starFactor) / starFactor;
    nebulaTexture.rgb *= 1.5 + 10.0 * lu_pow2(lu_max0(lu_GetStarNoise(starCoord) * lu_GetStarNoise(starCoord + 0.1) - 0.6));

    #if LU_NIGHT_NEBULA_I != 100
        #define LU_NIGHT_NEBULA_IM LU_NIGHT_NEBULA_I * 0.01
        nebulaTexture.a *= LU_NIGHT_NEBULA_IM;
    #endif

    #ifdef ATM_COLOR_MULTS
        nebulaTexture.rgb *= sqrtAtmColorMult;
    #endif

    return max(nebulaTexture.rgb * nebulaTexture.a, vec3(0.0));
}

//=======================================================================================
// LEGACY COMPATIBILITY
//=======================================================================================

// Maintain backward compatibility with existing celestial systems
#define LU_GetNightStars GetStars
#define LU_GetCelestialEffects GetAuroraBorealis

#endif // LU_INCLUDE_CELESTIAL_EFFECTS
