/* Prisma - end_space.glsl
End sky split into two layers:
- prismaEndFogSky: fog/background only
- prismaEndSpaceSky: visible sky with full-sphere stars and a stylized black hole
Inspired by the End atmosphere work in Solas and Iteration T.
*/

const float PRISMA_PI = 3.14159265358979323846;

float prismaHash12(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);
}

float prismaHash13(vec3 p) {
    p = fract(p * 0.1031);
    p += dot(p, p.zyx + 31.32);
    return fract((p.x + p.y) * p.z);
}

vec2 prismaRotate(vec2 p, float a) {
    float s = sin(a);
    float c = cos(a);
    return mat2(c, -s, s, c) * p;
}

float prismaNoise2(vec2 p) {
    vec2 i = floor(p);
    vec2 f = fract(p);
    f = f * f * (3.0 - 2.0 * f);

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

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

float prismaNoise3(vec3 p) {
    vec3 i = floor(p);
    vec3 f = fract(p);
    f = f * f * (3.0 - 2.0 * f);

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

    float x00 = mix(n000, n100, f.x);
    float x10 = mix(n010, n110, f.x);
    float x01 = mix(n001, n101, f.x);
    float x11 = mix(n011, n111, f.x);

    float y0 = mix(x00, x10, f.y);
    float y1 = mix(x01, x11, f.y);

    return mix(y0, y1, f.z);
}

float prismaFbm2(vec2 p) {
    float f = 0.0;
    float a = 0.5;
    for (int i = 0; i < 4; i++) {
        f += prismaNoise2(p) * a;
        p = p * 2.03 + vec2(7.1, 3.7);
        a *= 0.5;
    }
    return f;
}

float prismaFbm3(vec3 p) {
    float f = 0.0;
    float a = 0.5;
    for (int i = 0; i < 4; i++) {
        f += prismaNoise3(p) * a;
        p = p * 2.01 + vec3(6.2, 3.7, 8.1);
        a *= 0.5;
    }
    return f;
}

vec2 prismaSphereUv(vec3 dir) {
    dir = normalize(dir);
    return vec2(
        atan(dir.z, dir.x) / (2.0 * PRISMA_PI) + 0.5,
        asin(clamp(dir.y, -1.0, 1.0)) / PRISMA_PI + 0.5
    );
}

vec3 prismaHoleDir() {
    return normalize(vec3(-0.72, 0.23, 0.65));
}

vec3 prismaHoleBasisX(vec3 centerDir) {
    vec3 upRef = abs(centerDir.y) > 0.85 ? vec3(1.0, 0.0, 0.0) : vec3(0.0, 1.0, 0.0);
    return normalize(cross(upRef, centerDir));
}

vec3 prismaHoleBasisY(vec3 centerDir, vec3 basisX) {
    return normalize(cross(centerDir, basisX));
}

vec2 prismaHoleCoord(vec3 dir) {
    vec3 centerDir = prismaHoleDir();
    vec3 basisX = prismaHoleBasisX(centerDir);
    vec3 basisY = prismaHoleBasisY(centerDir, basisX);
    return vec2(dot(dir, basisX), dot(dir, basisY));
}

vec3 prismaWarpStarsDir(vec3 worldDir) {
    vec3 centerDir = prismaHoleDir();
    vec3 basisX = prismaHoleBasisX(centerDir);
    vec3 basisY = prismaHoleBasisY(centerDir, basisX);

    float facing = clamp(dot(worldDir, centerDir), 0.0, 1.0);
    vec2 p = vec2(dot(worldDir, basisX), dot(worldDir, basisY));
    float dist = length(p) + 1e-4;
    float lens = facing * facing * 0.11 / (1.0 + dist * 26.0);
    vec2 warped = p * (1.0 + lens);

    float axial = dot(worldDir, centerDir);
    return normalize(centerDir * axial + basisX * warped.x + basisY * warped.y);
}

float prismaStarLayer(vec2 uv, float scale, float threshold) {
    vec2 grid = uv * scale;
    vec2 cell = floor(grid);
    vec2 local = fract(grid) - 0.5;

    float rnd = prismaHash12(cell);
    float spawn = smoothstep(threshold, 1.0, rnd);

    vec2 jitter = vec2(
        prismaHash12(cell + vec2(2.31, 7.19)),
        prismaHash12(cell + vec2(8.11, 1.73))
    ) - 0.5;

    float size = mix(0.016, 0.090, prismaHash12(cell + vec2(6.7, 4.2)));
    float d = length(local - jitter * 0.55);
    float core = 1.0 - smoothstep(size, size * 2.9, d);
    core *= core;

    float sparkle = mix(0.74, 1.30, prismaHash12(cell + vec2(9.2, 3.4)));
    return spawn * core * sparkle;
}

vec3 prismaEndFogSky(vec3 worldDir) {
    worldDir = normalize(worldDir);

    float horizon = pow(clamp(1.0 - abs(worldDir.y), 0.0, 1.0), 1.20);
    float vertical = clamp(worldDir.y * 0.5 + 0.5, 0.0, 1.0);

    vec3 upperColor = vec3(0.015, 0.019, 0.034);
    vec3 lowerColor = vec3(0.013, 0.015, 0.028);
    vec3 horizonColor = vec3(0.037, 0.043, 0.062);

    vec3 color = mix(lowerColor, upperColor, vertical);
    color = mix(color, horizonColor, horizon * 0.82);

    return clamp(color, vec3(0.0), vec3(20.0));
}

vec3 prismaEndNebula(vec3 worldDir) {
    worldDir = normalize(worldDir);

    vec3 pA = worldDir;
    pA.xz = prismaRotate(pA.xz, 0.68);
    pA.xy = prismaRotate(pA.xy, -0.34);

    vec3 pB = worldDir.zxy;
    pB.xy = prismaRotate(pB.xy, -0.92);
    pB.yz = prismaRotate(pB.yz, 0.57);

    vec3 pC = worldDir.yzx;
    pC.xz = prismaRotate(pC.xz, 1.21);
    pC.xy = prismaRotate(pC.xy, 0.44);

    vec3 pD = worldDir;
    pD.yz = prismaRotate(pD.yz, -1.08);
    pD.xy = prismaRotate(pD.xy, 0.93);

    vec3 pSoft = worldDir;
    pSoft.xy = prismaRotate(pSoft.xy, 0.27);
    pSoft.yz = prismaRotate(pSoft.yz, -0.51);

    float nA = prismaFbm3(pA * 2.5 + vec3(0.0, frameTimeCounter * 0.003, -frameTimeCounter * 0.002));
    float nB = prismaFbm3(pB * 3.9 + vec3(3.1, -2.4, frameTimeCounter * 0.0025));
    float nC = prismaFbm3(pC * 1.7 - vec3(frameTimeCounter * 0.0015, 4.2, 1.7));
    float nD = prismaFbm3(pD * 2.9 + vec3(-2.1, 1.6, frameTimeCounter * 0.0018));
    float nSoft = prismaFbm3(pSoft * 1.15 + vec3(2.7, -1.5, frameTimeCounter * 0.0012));

    float bandA = pow(clamp(1.0 - abs(dot(worldDir, normalize(vec3(0.24, 0.83, -0.50)))), 0.0, 1.0), 0.72);
    float bandB = pow(clamp(1.0 - abs(dot(worldDir, normalize(vec3(-0.62, 0.18, 0.76)))), 0.0, 1.0), 0.86);
    float bandC = pow(clamp(1.0 - abs(dot(worldDir, normalize(vec3(0.71, -0.34, 0.61)))), 0.0, 1.0), 0.78);
    float horizon = pow(clamp(1.0 - abs(worldDir.y) * 0.55, 0.0, 1.0), 0.60);
    float fullSphere = max(0.46 + 0.54 * bandA, 0.40 + 0.60 * bandB);

    float cloudMask = smoothstep(0.42, 0.74, nA * 0.54 + nB * 0.46);
    float veilMask  = smoothstep(0.38, 0.76, nC * 0.46 + nD * 0.54);
    float knotMask  = smoothstep(0.53, 0.82, nB * 0.36 + nD * 0.64);
    float dustMask  = smoothstep(0.34, 0.72, nC * 0.60 + nSoft * 0.40);
    float broadMask = smoothstep(0.30, 0.66, nSoft);

    float blueMask   = max(cloudMask * (0.42 + 0.58 * fullSphere), broadMask * (0.18 + 0.30 * horizon));
    float violetMask = max(smoothstep(0.42, 0.76, nB * 0.64 + nD * 0.36) * (0.42 + 0.58 * bandB), broadMask * bandB * 0.34);
    float roseMask   = max(smoothstep(0.36, 0.72, nC * 0.58 + nB * 0.42) * (0.34 + 0.66 * bandB) * (0.46 + 0.54 * horizon), broadMask * (0.24 + 0.38 * bandC));
    float goldMask   = smoothstep(0.52, 0.82, nD * 0.54 + nA * 0.46) * (0.12 + 0.44 * horizon) * (0.18 + 0.82 * bandA);

    vec3 blueNebula   = mix(vec3(0.07, 0.12, 0.28), vec3(0.10, 0.24, 0.46), clamp(nA * 0.82 + nC * 0.18, 0.0, 1.0));
    vec3 violetNebula = mix(vec3(0.14, 0.08, 0.28), vec3(0.31, 0.11, 0.42), clamp(nB * 0.72 + nD * 0.28, 0.0, 1.0));
    vec3 roseNebula   = mix(vec3(0.24, 0.10, 0.22), vec3(0.42, 0.15, 0.30), clamp(nC * 0.66 + nB * 0.34, 0.0, 1.0));
    vec3 goldNebula   = mix(vec3(0.18, 0.15, 0.07), vec3(0.32, 0.23, 0.09), clamp(nD * 0.70 + nA * 0.30, 0.0, 1.0));

    vec3 nebula = vec3(0.0);
    nebula += blueNebula   * blueMask   * (0.42 + 0.28 * horizon);
    nebula += violetNebula * violetMask * (0.34 + 0.24 * fullSphere);
    nebula += roseNebula   * roseMask   * (0.20 + 0.24 * horizon);
    nebula += goldNebula   * goldMask   * 0.16;

    vec3 violetRoseWisps = mix(violetNebula, roseNebula, clamp(nSoft * 0.72 + 0.08, 0.0, 1.0));
    vec3 blueVioletWisps = mix(blueNebula, violetNebula, clamp(nD * 0.58 + 0.12, 0.0, 1.0));
    nebula += blueVioletWisps * veilMask * (0.14 + 0.12 * fullSphere);
    nebula += violetRoseWisps * broadMask * bandB * (0.18 + 0.10 * horizon);
    nebula += vec3(0.06, 0.08, 0.11) * dustMask * (0.18 + 0.18 * fullSphere);
    nebula += vec3(0.03, 0.02, 0.04) * knotMask * 0.10;

    return nebula;
}

vec3 prismaEndStars(vec3 worldDir, vec2 uv) {
    float starsA = prismaStarLayer(uv + vec2(0.13, 0.07), 320.0, 0.9882);
    float starsB = prismaStarLayer(uv * vec2(1.13, 0.91) + vec2(0.41, 0.23), 540.0, 0.9925);
    float starsC = prismaStarLayer(uv * vec2(0.77, 1.19) + vec2(0.73, 0.61), 180.0, 0.9852) * 0.58;

    float sphereFade = 0.86 + 0.14 * smoothstep(0.0, 0.05, abs(worldDir.y));

    vec3 starTintA = vec3(0.72, 0.82, 1.00);
    vec3 starTintB = vec3(1.00, 0.92, 0.82);
    vec3 starTintC = vec3(0.82, 0.92, 1.00);

    return sphereFade * (starsA * starTintA * 2.10 + starsB * starTintB * 2.95 + starsC * starTintC * 1.60);
}

float prismaBlackHoleFacing(vec3 worldDir) {
    return smoothstep(0.14, 0.34, dot(normalize(worldDir), prismaHoleDir()));
}

float prismaBlackHoleCoreMask(vec3 worldDir) {
    float dist = length(prismaHoleCoord(worldDir));
    float coreRadius = 0.168;
    float core = 1.0 - smoothstep(coreRadius - 0.012, coreRadius + 0.022, dist);
    return clamp(core * prismaBlackHoleFacing(worldDir), 0.0, 1.0);
}

float prismaBlackHoleShadowMask(vec3 worldDir) {
    float dist = length(prismaHoleCoord(worldDir));
    float shadow = 1.0 - smoothstep(0.17, 0.42, dist);
    return clamp(shadow * prismaBlackHoleFacing(worldDir), 0.0, 1.0);
}

vec3 prismaBlackHoleGlow(vec3 worldDir) {
    float facing = prismaBlackHoleFacing(worldDir);
    if (facing <= 0.0) return vec3(0.0);

    vec2 p = prismaRotate(prismaHoleCoord(worldDir), -0.58);
    float skyDist = length(p);
    float coreRadius = 0.168;

    float angle = atan(p.y, p.x);
    float outerMask = smoothstep(coreRadius + 0.008, coreRadius + 0.060, skyDist) * (1.0 - smoothstep(0.44, 0.62, skyDist));
    float sideFade = exp(-8.5 * abs(p.x));

    float arcCurve = coreRadius + 0.016 + 0.018 * exp(-abs(p.x) * 13.0);
    float lensedArc = exp(-540.0 * abs(abs(p.y) - arcCurve)) * sideFade;

    vec2 diskCoord = vec2(p.x * 1.15, p.y * 6.1);
    float diskBandA = exp(-72.0 * abs(length(diskCoord) - (coreRadius + 0.090)));
    float diskBandB = exp(-50.0 * abs(length(diskCoord) - (coreRadius + 0.135)));
    float equator = exp(-250.0 * abs(p.y)) * (diskBandA * 0.95 + diskBandB * 0.55) * sideFade;

    float detailA = prismaFbm2(vec2(angle * 3.5 + frameTimeCounter * 0.022, skyDist * 28.0 - frameTimeCounter * 0.018));
    float detailB = prismaFbm2(vec2(angle * 7.2 - frameTimeCounter * 0.015, skyDist * 48.0 + 3.7));
    float detail = clamp(detailA * 0.60 + detailB * 0.40, 0.0, 1.0);
    detail = 0.58 + detail * 0.70;

    float innerOcclusion = smoothstep(coreRadius + 0.008, coreRadius + 0.032, skyDist);
    float diskShape = max(lensedArc * 1.12, equator * 0.68);
    diskShape *= outerMask * innerOcclusion;

    vec3 diskCool = vec3(0.78, 0.82, 0.88);
    vec3 diskWarm = vec3(0.96, 0.92, 0.86);
    vec3 diskColor = mix(diskCool, diskWarm, clamp(detailA * 0.55 + 0.12, 0.0, 1.0));

    vec3 accretion = diskColor * diskShape * detail * 1.00;

    float photonRing = exp(-420.0 * abs(skyDist - (coreRadius + 0.010)));
    vec3 photonColor = vec3(0.90, 0.95, 1.00) * photonRing * 0.10;

    float softHalo = (1.0 - smoothstep(0.17, 0.36, skyDist)) * (1.0 - smoothstep(coreRadius + 0.006, coreRadius + 0.038, skyDist));
    vec3 halo = vec3(0.015, 0.019, 0.026) * softHalo * 0.42;

    return (accretion + photonColor + halo) * facing;
}

vec3 prismaEndSpaceSky(vec3 worldDir) {
    worldDir = normalize(worldDir);
    vec2 skyUv = prismaSphereUv(worldDir);

    vec3 color = prismaEndFogSky(worldDir);
    color += prismaEndNebula(worldDir);
    color += prismaEndStars(worldDir, skyUv);

    return clamp(color, vec3(0.0), vec3(20.0));
}
