#version 300 es
precision highp float;

uniform float iTime;
uniform vec2 iResolution;
uniform int uLevels;
out vec4 fragColor;

// -------------------------------------------------------------
// Rotations
// -------------------------------------------------------------
mat3 rotX(float a){ float c=cos(a),s=sin(a); return mat3(1.0,0.0,0.0, 0.0,c,-s, 0.0,s,c); }
mat3 rotY(float a){ float c=cos(a),s=sin(a); return mat3(c,0.0,s, 0.0,1.0,0.0, -s,0.0,c); }
mat3 rotZ(float a){ float c=cos(a),s=sin(a); return mat3(c,-s,0.0, s,c,0.0, 0.0,0.0,1.0); }

// -------------------------------------------------------------
// SDFs
// -------------------------------------------------------------
float sdCrispBox(vec3 p, vec3 b){
    vec3 q = abs(p) - b;
    return max(max(q.x,q.y), q.z);
}
float sdRoundBox(vec3 p, vec3 b, float r){
    vec3 q = abs(p) - b;
    return length(max(q, 0.0)) - r;
}

// -------------------------------------------------------------
// Smooth min (melting)
// -------------------------------------------------------------
float smin(float a,float b,float k){
    float h = clamp(0.5 + 0.5*(b-a)/k, 0.0, 1.0);
    return mix(b,a,h) - k*h*(1.0-h);
}
float contactMelt(float da,float db,float baseR,float meltAmp){
    float overlap = da + db;
    float r = baseR * meltAmp;
    float blend = smoothstep(r,0.0,overlap);
    float k = r * blend;
    return smin(da, db, k);
}
vec2 opCM(vec2 a, vec2 b, float baseR, float meltAmp){
    float d = contactMelt(a.x, b.x, baseR, meltAmp);
    return (d == a.x) ? vec2(d,a.y) : vec2(d,b.y);
}

// -------------------------------------------------------------
// Helpers
// -------------------------------------------------------------
float hash(vec3 p){
    p = fract(p*0.3183099 + vec3(0.1,0.2,0.3));
    p += dot(p,p.yzx+19.19);
    return fract(p.x*p.y*p.z);
}
float meltPhase(){
    float w = 0.5 + 0.5 * sin(iTime * 1.1);
    w = pow(w, 3.0);
    return clamp(w, 0.0, 1.0);
}

// level-based selector in [0,1)
float shapeSelector(float f){
    return fract(sin(f * 12.345) * 7654.321);
}

// mild domain wobble that keeps structure readable
vec3 softWobble(vec3 p, float t){
    p += 0.05 * sin(p * 2.3 + t * 0.9);
    p += 0.03 * sin(p.yzx * 3.1 - t * 0.6);
    return p;
}

// kaleidoscopic fold (original-ish)
vec3 kaleido(vec3 p){
    p.xy = abs(p.xy);
    p.xz = abs(p.xz);
    float r = length(p.xy);
    float k = 0.4 + 0.6 * sin(iTime*0.7 + r*2.5);
    p.xy *= mix(1.0, 0.5 + 0.5*sin(r*3.0), 0.3 * k);
    return p;
}

// -------------------------------------------------------------
// Scene with branching appendages + recursive sub-branches
// -------------------------------------------------------------
vec2 scene(vec3 p){
    float t  = iTime;
    float mp = meltPhase();

    // global kaleidoscopic warp + gentle wobble
    p = kaleido(p * 0.95);
    p = softWobble(p, t);

    // gentle global swirl
    float swirl = 0.3 * sin(t*0.25);
    mat3 G = rotY(swirl) * rotZ(swirl*0.6);
    p = G * p;

    const int MAX_L = 9;
    int levels = uLevels;
    if(levels <= 0) levels = MAX_L;
    levels = min(levels, MAX_L);

    vec2 res = vec2(1e9, 0.0);

    float size  = 0.6;
    float scale = 0.45;
    float meltAmpGlobal = 0.8 + 6.0 * pow(mp, 2.3);

    for(int level=0; level<MAX_L; level++){
        if(level >= levels) break;

        float lf = float(level);
        vec3 lp = p;

        // gentle twist
        float ay = sin(lf*0.9 + t*0.4) * (0.30 + 0.22*mp);
        float ax = sin(lf*1.2 + t*0.27) * (0.22 + 0.22*mp);
        float az = sin(lf*0.7 + t*0.5) * (0.18 + 0.18*mp);
        mat3 twist = rotY(ay)*rotX(ax)*rotZ(az);
        lp = twist * lp;

        float grow = 0.25 + 0.75 * abs(sin(t*0.7 + lf*1.37));

        // ------ Parent primitive with controlled variation ------
        float shapeSel = shapeSelector(lf*1.37);
        vec3 halfSize = vec3(size*grow);

        // tiny per-level stretch on Y for visual difference
        halfSize.y *= mix(1.0, 1.3, shapeSel);

        float parentCrisp = sdCrispBox(lp, halfSize);
        float parentRound = sdRoundBox(lp, halfSize,
                                       size*grow*0.25*mp*(0.3 + 0.7*shapeSel));
        // mix between crisp and round depending on level + melt
        float parent = mix(parentCrisp, parentRound,
                           mp * (0.4 + 0.6*shapeSel));

        float baseR = mix(0.06, 0.3, lf/float(MAX_L-1));
        res = opCM(res, vec2(parent, lf), baseR, meltAmpGlobal);

        // ---------- Branch / appendage setup ----------
        float childSize = size * scale * grow;
        float shift = size*grow + childSize;

        // branches get thinner as we go out through levels
        float branchThin = mix(0.9, 0.35, lf / float(MAX_L - 1));
        float branchAngle = 0.35 * (0.3 + 0.7 * mp);

        const int MAX_DIRS = 24;
        vec3 dirs[MAX_DIRS];
        int dCount = 0;

        // Original 6 face directions
        dirs[dCount++] = vec3( shift,0.0,0.0);
        dirs[dCount++] = vec3(-shift,0.0,0.0);
        dirs[dCount++] = vec3(0.0, shift,0.0);
        dirs[dCount++] = vec3(0.0,-shift,0.0);
        dirs[dCount++] = vec3(0.0,0.0, shift);
        dirs[dCount++] = vec3(0.0,0.0,-shift);

        // Additional tilted directions (appendages)
        for(int bi = 0; bi < 3; bi++){
            float a = (lf*1.37 + float(bi)*2.1 + t*0.3);
            float s = sin(a), c = cos(a);

            vec3 baseDir = normalize(vec3(c, 0.0, s));
            vec3 up      = vec3(0.0, 1.0, 0.0);
            vec3 side    = normalize(cross(up, baseDir) + 1e-3);

            vec3 upTilt    = normalize(baseDir + branchAngle * up);
            vec3 downTilt  = normalize(baseDir - branchAngle * up);
            vec3 sideTilt1 = normalize(baseDir + branchAngle * side);
            vec3 sideTilt2 = normalize(baseDir - branchAngle * side);

            if(dCount < MAX_DIRS) dirs[dCount++] = upTilt    * shift;
            if(dCount < MAX_DIRS) dirs[dCount++] = downTilt  * shift;
            if(dCount < MAX_DIRS) dirs[dCount++] = sideTilt1 * shift;
            if(dCount < MAX_DIRS) dirs[dCount++] = sideTilt2 * shift;
        }

        // ---------- First-generation branches / children ----------
        for(int k=0; k<MAX_DIRS; k++){
            if(k >= dCount) break;

            vec3 cq = lp - dirs[k];

            // slightly softer wobble on children
            cq += 0.08 * mp * sin(cq*3.1 + (lf + float(k))*1.13 + t*2.7);

            vec3 chS = vec3(childSize * branchThin);
            // a bit of anisotropy for variety
            float childSel = shapeSelector(lf*3.17 + float(k)*1.11);
            chS.xz *= mix(1.0, 0.85, childSel);

            float roundAmt = mix(0.15, 0.35, childSel);
            float cb = sdRoundBox(cq, chS, chS.x * roundAmt * mp);

            // Slightly stronger melt between branches
            res = opCM(res, vec2(cb, lf), baseR * 1.25, meltAmpGlobal * 1.05);
        }

        // ---------- Second-generation recursive sub-branches ----------
        int subCount = 4;           // moderate amount
        float subLen  = childSize * 0.75;
        float subThin = childSize * 0.45;  // thinner than first-gen

        for(int sk = 0; sk < subCount; sk++){
            float h = hash(vec3(float(sk), lf, 3.17));
            int parentDirIndex = int(floor(h * float(max(dCount-1, 1))));
            parentDirIndex = clamp(parentDirIndex, 0, max(dCount-1, 0));

            vec3 baseDir = normalize(dirs[parentDirIndex]);

            float ang = h * 6.28318 + t * 0.3;
            vec3 randAxis = normalize(vec3(sin(ang), cos(ang*1.3), sin(ang*0.7)));

            vec3 subDir = normalize(
                mix(baseDir, randAxis, 0.35 + 0.35 * mp)
            ) * (shift + subLen);

            vec3 sq = lp - subDir;

            // gentle curving / bending
            sq += 0.06 * sin(sq * 3.1 + t*1.7 + float(sk)*1.2);

            vec3 sbSize = vec3(subThin, subThin * 0.7, subThin);
            float roundR = sbSize.x * (0.25 + 0.35*mp);

            float sd = sdRoundBox(sq, sbSize, roundR);

            // slightly increased melt radius between branch systems
            res = opCM(res, vec2(sd, lf + 0.4), baseR * 1.35, meltAmpGlobal * 1.05);
        }

        // fold to nearest branch / child for next level
        float best = 1e9;
        vec3 bestP = lp;
        for(int k=0; k<MAX_DIRS; k++){
            if(k >= dCount) break;
            vec3 pp = lp - dirs[k];
            float m = dot(pp,pp);
            if(m < best){ best = m; bestP = pp; }
        }

        p = bestP;
        size *= scale;
    }

    return res;
}

// -------------------------------------------------------------
// Raymarch + shading
// -------------------------------------------------------------
vec2 raymarch(vec3 ro, vec3 rd){
    float t = 0.0;
    for(int i=0;i<220;i++){
        vec3 pos = ro + rd*t;
        vec2 d = scene(pos);
        if(d.x < 0.001) return vec2(t,d.y);
        t += d.x;
        if(t > 60.0) break;
    }
    return vec2(-1.0, -1.0);
}
vec3 getNormal(vec3 p){
    float h=0.001;
    const vec3 e=vec3(1.0,-1.0,0.0);
    return normalize(
        e.xyy*scene(p+e.xyy*h).x +
        e.yyx*scene(p+e.yyx*h).x +
        e.yxy*scene(p+e.yxy*h).x +
        e.xxx*scene(p+e.xxx*h).x
    );
}

// rainbow-ish palette (slightly tuned)
vec3 palette(float t){
    vec3 a = vec3(0.52, 0.48, 0.50);
    vec3 b = vec3(0.48, 0.52, 0.50);
    vec3 c = vec3(1.0, 0.7, 0.4);
    vec3 d = vec3(0.0, 0.33, 0.67);
    return a + b * cos(6.28318 * (c * t + d));
}
// thin-film interference (pearlescent)
vec3 thinFilm(float angle, float shift){
    float tt = angle * 1.4 + shift * 2.5;
    return vec3(
        0.5 + 0.5 * sin(tt),
        0.5 + 0.5 * sin(tt + 2.094),
        0.5 + 0.5 * sin(tt + 4.188)
    );
}
// -------------------------------------------------------------
// Skyglobe based on scene's aesthetics
// -------------------------------------------------------------
vec3 skyColor(vec3 rd){
    float mp = meltPhase();
    vec3 p = kaleido(rd * 2.5);
    float ang = atan(rd.z, rd.x);
    float lat = rd.y;
    float bands = 0.5 + 0.5 * sin(ang * 4.0 + iTime*0.3 + p.x*3.0);
    float rings = 0.5 + 0.5 * sin(length(p.xy)*6.0 + iTime*0.5 + p.z*2.0);
    float t = bands*0.4 + rings*0.6 + lat*0.3 + mp*0.25;
    vec3 col = palette(t);
    float h = 0.5 + 0.5*rd.y;
    vec3 skyTop = col * (1.1 + 0.5*h);
    vec3 skyBottom = col * 0.35 + vec3(0.02,0.03,0.05);
    col = mix(skyBottom, skyTop, smoothstep(-0.2, 0.7, rd.y));
    float desat = smoothstep(0.3, 1.0, rd.y);
    float l = dot(col, vec3(0.299, 0.587, 0.114));
    col = mix(col, vec3(l), desat * 0.25);
    return col;
}
// -------------------------------------------------------------
// Main
// -------------------------------------------------------------
void main(){
    vec2 uv = (gl_FragCoord.xy / iResolution - 0.5);
    uv.x *= iResolution.x / iResolution.y;
    // slight camera orbit
    float camT = iTime * 0.25;
    float camR = 4.0;
    vec3 ro = vec3(
        camR * sin(camT),
        1.2 + 0.6*sin(camT*1.7),
        camR * cos(camT)
    );
    vec3 ta = vec3(0.0, 0.0, 0.0);
    vec3 ww = normalize(ta-ro);
    vec3 uu = normalize(cross(vec3(0.0,1.0,0.0), ww));
    vec3 vv = cross(ww, uu);
    vec3 rd = normalize(uv.x*uu + uv.y*vv + ww);
    vec2 hit = raymarch(ro, rd);
    if(hit.x < 0.0){
        vec3 bg = skyColor(rd);
        fragColor = vec4(bg,1.0);
        return;
    }
    vec3 pos = ro + rd*hit.x;
    vec3 nor = getNormal(pos);
    // base lighting
    vec3 lightDir = normalize(vec3(0.7,1.0,0.4));
    float diff = max(dot(nor, lightDir),0.0);
    float mp = meltPhase();
    // iridescent color based on level, position and normal
    float lvl = hit.y;
    float nDot = dot(nor, normalize(vec3(0.3,0.6,0.7)));
    float pHash = hash(pos*1.7);
    float colorT = lvl*0.18 + iTime*0.12 + nDot*0.5 + pHash*0.2;
    vec3 baseCol = palette(colorT);
    // subtle stripe bands along view direction
    float stripes = 0.5 + 0.5*sin(dot(pos, vec3(0.8,1.3,0.9))*1.3 + iTime*1.5);
    baseCol *= mix(0.8, 1.4, stripes);
    // thin-film pearlescent highlight
    vec3 tf = thinFilm(nDot, lvl*0.2 + iTime*0.1);
    baseCol = mix(baseCol, tf, 0.25 + 0.25*mp);
    // melt glow (subsurface-ish)
    vec3 meltGlow = vec3(1.0,0.7,0.9) * pow(mp, 1.8);
    baseCol += meltGlow * 0.4;
    // hemisphere ambient (sky vs ground)
    float hemi = 0.5 + 0.5*nor.y;
    // view-dependent sheen
    float viewSheen = pow(max(dot(reflect(-lightDir, nor), rd),0.0), 6.0);
    baseCol += viewSheen * vec3(1.2, 0.95, 0.85) * 0.3;
    // subtle glitter
    float glitter = pow(max(0.0, sin(dot(pos, nor)*15.0 + iTime*4.0)), 8.0);
    vec3 glitterCol = vec3(1.2, 0.9, 0.8) * glitter * 0.5;
    vec3 col = baseCol * (0.25 + 0.85*diff + 0.25*hemi) + glitterCol;
    // distance fog
    float fog = exp(-0.08 * hit.x);
    vec3 bgFog = vec3(0.02,0.03,0.05);
    col = mix(bgFog, col, fog);
    // gentle chromatic vignette
    float r2 = dot(uv, uv);
    float v = smoothstep(1.1, 0.25, r2);
    col *= vec3(v, v*0.97, v*0.92);
    fragColor = vec4(col,1.0);
}