/*
Event Horizon relativistic renderer for Iris / GLSL 1.30.

GENERATED from blackhole_engine.py (SHA-256 7ba6cfa881a363cbe4a1407ef70f82bdb6025d1fa2d58783220a2fcaa7885c83).
Do not hand-edit the simulation core; run tools/port_blackhole.py instead.

The noise, black-body emission, accretion density, plunging region, Doppler /
gravitational redshift, photon ring and RK4 geodesic loop below come directly
from SCENE_FS_TMPL in the supplied Python simulation. Only the camera adapter,
Minecraft sky compositing and compile-time profile selection are Iris-specific.
*/

#if BLACK_HOLE_QUALITY == 1
    #define OCT_MAIN 3
    #define OCT_SUB 2
    #define OCT_PLUNGE 3
    #define OCT_NEB 2
    #define STAR_LAYERS 1
    #define USE_SHADOW 0
    #define USE_DISPERSION 0
    #define USE_VEIL 0
    #define DT_SCALE 0.225
    #define BH_STEPS 260
#elif BLACK_HOLE_QUALITY == 2
    #define OCT_MAIN 4
    #define OCT_SUB 3
    #define OCT_PLUNGE 4
    #define OCT_NEB 3
    #define STAR_LAYERS 2
    #define USE_SHADOW 1
    #define USE_DISPERSION 0
    #define USE_VEIL 1
    #define DT_SCALE 0.165
    #define BH_STEPS 440
#else
    #define OCT_MAIN 5
    #define OCT_SUB 4
    #define OCT_PLUNGE 5
    #define OCT_NEB 4
    #define STAR_LAYERS 3
    #define USE_SHADOW 1
    #define USE_DISPERSION 1
    #define USE_VEIL 1
    #define DT_SCALE 0.120
    #define BH_STEPS 720
#endif

// SceneSettings defaults from blackhole_engine.py.
#define uBhScale 1.10
#define uDiskIn 2.55
#define uDiskOut 20.0
#define uDiskH 0.055
#define uDiskBright 3.50
#define uTemp 12500.0
#define uTempGrad 0.72
#define uDopp 0.24
#define uSpin 0.40
#define uAbsorb 0.92
#define uFil 1.48
#define uPatch 1.12
#define uConc 1.65
#define uPlunge 0.76
#define uDust 0.92
#define uRing 0.72
#define uCore 0.24
#define uHalo 0.072
#define uStars 0.0015
#define uSat 0.36
#define uTime ehAnimationTime
#define uSteps BH_STEPS
#define uCamPos BH_CAM_POS

#define RS (uBhScale)
#define M (0.5 * uBhScale)
#define RPLUNGE (1.045 * uBhScale)
#define RPHOTON (1.5 * uBhScale)
const float ESCAPE = 90.0;
const float TNORM = 0.488;

const vec3 BH_CAM_POS = vec3(16.299244, 1.814519, 36.481541);
float ehAnimationTime = 0.0;

float hash13(vec3 p3){
    p3 = fract(p3 * 0.1031);
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.x + p3.y) * p3.z);
}
vec2 hash23(vec3 p3){
    p3 = fract(p3 * vec3(0.1031, 0.1030, 0.0973));
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.xx + p3.yz) * p3.zy);
}
float vnoise(vec3 x){
    vec3 i = floor(x), f = fract(x);
    f = f * f * (3.0 - 2.0 * f);
    return mix(mix(mix(hash13(i+vec3(0,0,0)), hash13(i+vec3(1,0,0)), f.x),
                   mix(hash13(i+vec3(0,1,0)), hash13(i+vec3(1,1,0)), f.x), f.y),
               mix(mix(hash13(i+vec3(0,0,1)), hash13(i+vec3(1,0,1)), f.x),
                   mix(hash13(i+vec3(0,1,1)), hash13(i+vec3(1,1,1)), f.x), f.y), f.z);
}
float fbm(vec3 p, int oct){
    float a = 0.5, s = 0.0;
    for(int i = 0; i < oct; i++){
        s += a * vnoise(p);
        p = p * 2.03 + vec3(11.3, 7.7, 3.1);
        a *= 0.5;
    }
    return s;
}
float ridge(vec3 p){
    float n = 1.0 - abs(vnoise(p) * 2.0 - 1.0);
    return n * n;
}
float rfbm(vec3 p, int oct){
    float a = 0.5, s = 0.0, w = 1.0, norm = 0.0;
    for(int i = 0; i < oct; i++){
        float r = ridge(p) * w;
        s += a * r;
        norm += a;
        w = clamp(r * 2.0, 0.0, 1.0);
        p = p * 2.14 + vec3(9.1, 4.7, 6.3);
        a *= 0.55;
    }
    return s / max(norm, 1e-4);
}

// ── Couleur corps noir : normalisée en LUMINANCE (chroma préservée) ───────
vec3 blackbody(float T){
    T = clamp(T, 900.0, 40000.0);
    vec3 c;
    c.x = 56100000.0 * pow(T, -1.5) + 148.0;
    c.y = (T > 6500.0) ? 35200000.0 * pow(T, -1.5) + 184.0
                       : 100.04 * log(T) - 623.6;
    c.z = 194.18 * log(T) - 1448.6;
    c = clamp(c, 0.0, 255.0) / 255.0;

    float l = dot(c, vec3(0.2126, 0.7152, 0.0722));
    c /= max(l, 1e-3);
    c = mix(vec3(1.0), c, clamp(uSat, 0.0, 3.0));
    return max(c, vec3(0.0));
}

// ── Fond ──────────────────────────────────────────────────────────────────
vec3 starfield(vec3 d){
    vec3 col = vec3(0.0);
    for(int k = 0; k < STAR_LAYERS; k++){
        float sc = 105.0 * exp2(float(k));
        vec3 p  = d * sc;
        vec3 id = floor(p);
        vec3 f  = fract(p) - 0.5;
        float h = hash13(id + vec3(37.0 * float(k) + 1.0));
        float thr = 0.979 - 0.005 * float(k);
        if (h > thr){
            vec2  r2  = hash23(id + 11.0);
            vec3  off = (vec3(r2, hash13(id + 5.0)) - 0.5) * 0.5;
            float dd  = length(f - off);
            float br  = pow((h - thr) / (1.0 - thr), 3.0);
            float core = exp(-dd * dd * 145.0);
            float glow = exp(-dd * dd * 11.0) * 0.075;
            float tint = hash13(id + 3.0);
            vec3  sc2  = mix(vec3(0.66,0.78,1.0), vec3(1.0,0.82,0.62), tint);
            col += sc2 * br * (core + glow) * (1.75 - 0.42 * float(k));
        }
    }
    float n = fbm(d * 2.2 + 4.0, OCT_NEB);
    col += mix(vec3(0.0032,0.0042,0.0105), vec3(0.0105,0.0055,0.0145),
               smoothstep(0.35, 0.75, n)) * smoothstep(0.38, 0.90, n);
    return col * uStars;
}

// ── Densité : disque képlérien + plunging region ──────────────────────────
float diskDensity(vec3 p, out float rr, out float surfaceDetail){
    rr = length(p.xz);
    surfaceDetail = 1.0;
    // Let a small amount of optically thin material escape beyond the nominal
    // disk radius.  The smooth noisy falloff below decides its actual edge.
    if (rr < RPLUNGE || rr > uDiskOut * 1.08) return 0.0;

    float lr = log(max(rr, 0.05));
    float diskT = clamp((rr - uDiskIn) / max(uDiskOut - uDiskIn, 1e-3), 0.0, 1.0);
    float omega = uSpin / (rr * sqrt(rr));
    // Differential Keplerian shear winds the procedural density field into
    // coherent spiral ribbons instead of screen-space-looking random blobs.
    float a  = -uTime * omega - mix(6.5, 2.8, diskT) * lr;
    float ca = cos(a), sa = sin(a);
    vec3  q  = vec3(ca * p.x - sa * p.z, p.y, sa * p.x + ca * p.z);
    vec3  dir = normalize(vec3(q.x, 0.0, q.z));
    float phi = atan(q.z, q.x);

    // ══ 1. PLUNGING REGION ════════════════════════════════════════════════
    if (rr < uDiskIn){
        float pl = clamp(uPlunge, 0.0, 2.0);
        if (pl <= 0.0) return 0.0;

        float u = clamp((uDiskIn - rr) / max(uDiskIn - RPLUNGE, 1e-3), 0.0, 1.0);

        float hp = uDiskH * (0.95 - 0.48 * u);
        float y  = p.y / max(hp, 1e-4);
        float vert = exp(-y * y * 1.5);
        if (vert < 0.0025) return 0.0;

        // spirales très cisaillées + turbulence fine
        vec3 nc = vec3(dir.x * 3.6, dir.z * 3.6, lr * 27.0 + p.y * 3.0);
        float f = rfbm(nc, OCT_PLUNGE);
        float turb = fbm(nc * vec3(2.1, 2.1, 1.6) + 19.0, 2);
        f = pow(clamp(f * 1.22 + turb * 0.18, 0.0, 1.0), 1.42);
        float dens = (0.52 + 0.48 * f);

        float fall = exp(-2.25 * u) * (1.0 - smoothstep(0.82, 1.0, u));
        float join = smoothstep(0.0, 0.16, u) * 0.85 + 0.15;

        surfaceDetail = clamp(0.52 + 0.92 * f, 0.35, 1.42);
        return vert * dens * fall * join * 1.62 * pl;
    }

    // ══ 2. DISQUE ═════════════════════════════════════════════════════════
    float t = diskT;

    // Keep the inner photosphere razor thin (the reference has a narrow white
    // spine) while strongly flaring only the dusty outer disk.  A single
    // thickness multiplier made the centre into a broad white rectangle as
    // soon as the arms became volumetric enough.
    float h = uDiskH * (1.0 + 5.20 * pow(t, 1.34));

    // ── instabilite magneto-rotationnelle (MRI) + flottabilite de Parker ──
    // Dans les regions externes le rapport pression magnetique / pression
    // gazeuse grimpe : les tubes de flux deviennent flottants et soulevent
    // aleatoirement la matiere hors du plan equatorial, vers le haut ET le
    // bas. On module donc l'altitude locale du plan (warp) et l'echelle de
    // hauteur (puff), sans toucher a l'interieur ou la gravite domine.
    float wEdge = smoothstep(0.26, 1.0, t);
    float mri1  = fbm(vec3(dir.x * 2.1, dir.z * 2.1, lr * 5.0) + 91.3, 3) - 0.5;
    float mri2  = fbm(vec3(dir.x * 5.4, dir.z * 5.4, lr * 9.6) + 47.9, 2) - 0.5;
    float warp  = (mri1 * 1.95 + mri2 * 0.95) * wEdge;
    float puff  = 1.0 + 1.35 * wEdge
                * smoothstep(0.32, 0.95,
                    fbm(vec3(dir.x * 3.2, dir.z * 3.2, lr * 6.4) + 13.9, 3));
    h *= puff;

    float y = (p.y - warp * uDiskH * 2.6) / h;
    float vert = exp(-y * y * 1.28);

#if USE_VEIL
    // voiles ténus soulevés au-dessus du plan : donne du volume aux bords
    float veil = exp(-abs(y) * 0.62) * 0.075
               * smoothstep(0.10, 0.75, t)
               * smoothstep(0.30, 0.85, fbm(vec3(dir.x*1.9, dir.z*1.9, lr*2.4)+7.1, 3));
    vert = max(vert, veil);
#endif
    // panaches magnetiques : la ou la pression magnetique l'emporte, des
    // langues de plasma sont ejectees verticalement (haut ET bas) puis
    // retombent -> bord du disque instable, jamais deux fois pareil.
    float pl0    = fbm(vec3(dir.x * 4.4, dir.z * 4.4, lr * 3.4) + 123.4, 3);
    float pl1    = fbm(vec3(dir.x * 9.1, dir.z * 9.1, lr * 7.7) + 205.1, 2);
    float tongue = smoothstep(0.52, 0.94, pl0 * 0.70 + pl1 * 0.30) * wEdge;
    float hplume = uDiskH * (2.5 + 10.0 * t);
    float plume  = tongue * exp(-abs(p.y) / max(hplume, 1e-4))
                 * (0.16 + 0.14 * pl1);
    vert = max(vert, plume);
    if (vert < 0.0025) return 0.0;

    float fr = mix(16.5, 6.0, t);
    float fa = mix(2.05, 3.45, t);
    vec3 nc = vec3(dir.x * fa, dir.z * fa, lr * fr) + vec3(0.0, 0.0, p.y * 2.0);
    // Domain-warped log-polar turbulence.  MRI stretches eddies along the
    // Keplerian flow, so the gas forms nested hair-thin ribbons rather than
    // isotropic procedural clouds.
    float domainWarp = fbm(vec3(dir.x * 1.8, dir.z * 1.8, lr * 4.1) + 83.2, 3) - 0.5;
    nc.z += domainWarp * mix(5.2, 2.4, t);
    float f1 = rfbm(nc, OCT_MAIN);
    float f2 = rfbm(nc * vec3(2.3, 2.3, 2.9) + 13.0, OCT_SUB);
    float fil = mix(f1, f1 * 0.60 + f2 * 0.40, 0.62);
    fil = pow(clamp(fil * 1.30, 0.0, 1.0), mix(1.0, 3.2, clamp(uFil, 0.0, 2.0)));

    float ribbonPhase = phi * mix(11.0, 6.0, t)
                      + lr * mix(23.0, 10.0, t) + domainWarp * 5.5;
    float ribbons = pow(0.5 + 0.5 * cos(ribbonPhase), mix(5.2, 2.0, t));
    float hairline = pow(0.5 + 0.5 * cos(ribbonPhase * 2.07 + f2 * 4.0), 8.0);
    float flowStructure = 0.48 + 0.72 * ribbons + 0.28 * hairline;
    fil *= mix(1.0, flowStructure, 0.64 * clamp(uFil, 0.0, 1.5));

    float base = mix(0.44, 0.035, smoothstep(0.05, 0.85, t));

    float cn = clamp(uConc, 0.0, 2.0);
    float radial = pow(uDiskIn / rr, mix(0.55, 1.95, cn * 0.5))
                 * exp(-mix(0.70, 2.30, cn * 0.5) * t);
    radial = mix(1.0, radial, clamp(cn, 0.0, 1.0)) * 1.72;

    // bourrelet interne
    float d0 = (rr - uDiskIn * 1.13) / (uDiskIn * 0.42);
    radial *= 1.0 + 1.40 * exp(-d0 * d0);

    // fragmentation externe
    float clump = fbm(vec3(dir.x * 2.9, dir.z * 2.9, lr * 3.1) + 31.7, OCT_SUB + 1);
    float mott  = fbm(vec3(dir.x * 7.5, dir.z * 7.5, lr * 8.4) + 57.3, OCT_SUB);
    float torn  = smoothstep(0.34, 0.66, clump * 0.72 + mott * 0.28);
    float amt   = smoothstep(0.22, 0.95, t) * clamp(uPatch, 0.0, 1.5);
    float mask  = mix(1.0, mix(0.06, 1.35, torn), clamp(amt, 0.0, 1.0));

    // rubans de poussière : creusent des sillons sombres, donnent la profondeur
    float lane = rfbm(vec3(dir.x * 1.35, dir.z * 1.35, lr * 4.6) + 71.0, 3);
    mask *= 1.0 - clamp(uDust, 0.0, 1.0) * 0.55
                * smoothstep(0.42, 0.86, lane) * smoothstep(0.04, 0.35, t);

    float inner = smoothstep(uDiskIn * 0.995, uDiskIn * 1.06, rr) * 0.8 + 0.2;
    // A noisy photospheric edge avoids a perfectly cut-out ellipse and leaves
    // sparse, magnetically lifted streamers at large radii.
    float edgeNoise = fbm(vec3(dir.x * 5.7, dir.z * 5.7, lr * 12.0) + 141.0, 3);
    float edgeCoord = rr / uDiskOut - (edgeNoise - 0.5) * 0.075;
    float outer = 1.0 - smoothstep(0.84, 1.035, edgeCoord);
    outer += 0.12 * tongue * smoothstep(0.72, 0.96, t)
           * (1.0 - smoothstep(0.96, 1.0, t));
    // Density alone stops carrying contrast once a ray reaches optical
    // thickness. Preserve the MRI temperature/emissivity texture separately
    // so bright strands and cool dust lanes survive radiative saturation.
    surfaceDetail = clamp(0.30 + 1.12 * fil, 0.34, 1.48);
    surfaceDetail *= mix(1.0, 0.68 + 0.52 * torn, clamp(amt, 0.0, 1.0));
    surfaceDetail *= 1.0 - clamp(uDust, 0.0, 1.0) * 0.30
                           * smoothstep(0.46, 0.84, lane);
    return vert * inner * outer * radial * mask
         * (base + (1.0 - base) * fil);
}

vec3 accel(vec3 p, float h2){
    float r2 = dot(p, p);
    float r5 = r2 * r2 * sqrt(r2);
    return -1.5 * RS * h2 * p / r5;
}

#if USE_SHADOW
// self-shadowing : quelques pas verticaux vers le « haut » du disque pour
// savoir combien de matière surplombe le point -> le disque cesse d'être
// une nappe plate uniformément lumineuse.
float shadowUp(vec3 p){
    float acc = 0.0;
    float dummy, detail;
    float sgn = (p.y >= 0.0) ? 1.0 : -1.0;
    for(int i = 1; i <= 4; i++){
        vec3 sp = p + vec3(0.0, sgn * uDiskH * 1.15 * float(i), 0.0);
        acc += diskDensity(sp, dummy, detail);
    }
    return exp(-acc * uDiskH * 1.15 * uAbsorb * 1.5);
}
#endif

vec3 EHTraceBlackHole(vec3 rd, vec2 uv, vec3 sceneColor, float animationTime) {
    vec3 p = uCamPos;
    vec3 v = normalize(rd);
    float dith = hash13(vec3(gl_FragCoord.xy, animationTime));
    float uSeed = fract(animationTime * 0.61803398875 + dith);

#if USE_DISPERSION
    // dispersion chromatique du lensing : chaque échantillon dévie très
    // légèrement la constante de déviation selon une longueur d'onde tirée
    // au sort ; l'accumulation temporelle fabrique un liseré spectral
    // sur l'anneau d'Einstein.
    float lam = fract(dith + uSeed * 0.618);
    // Gravity is achromatic; this tiny residual models only the camera optics.
    float bend = 1.0 + (lam - 0.5) * 0.003;
    vec3  lamTint = normalize(mix(vec3(1.25,0.85,0.60), vec3(0.70,0.90,1.35), lam)) * 1.732;
#else
    float bend = 1.0;
    vec3  lamTint = vec3(1.0);
#endif

    vec3  L  = cross(p, v);
    float h2 = dot(L, L) * bend;

    vec3  col   = vec3(0.0);
    vec3  cor   = vec3(0.0);
    float trans = 1.0;
    bool  eaten = false;
    bool  escaped = false;
    float ringAcc = 0.0;
    float minRadius = 1e9;
    float orbitSweep = 0.0;

    for(int i = 0; i < uSteps; i++){
        float r = length(p);
        if (r < RS * 1.003) { eaten = true; break; }
        if (r > ESCAPE && dot(p, v) > 0.0) { escaped = true; break; }
        if (trans < 0.003) break;

        float dt = DT_SCALE * max(r - RS, 0.02) * (0.35 + 0.65 * smoothstep(1.0, 6.0, r));
        dt = min(dt, 0.5);
        if (r < uDiskOut * 1.3){
            float dy = abs(p.y);
            dt = min(dt, max(uDiskH * 0.58, dy * 0.40));
        }
        dt = max(dt, 0.004);

        // RK4
        vec3 k1p = v,               k1v = accel(p, h2);
        vec3 k2p = v + 0.5*dt*k1v,  k2v = accel(p + 0.5*dt*k1p, h2);
        vec3 k3p = v + 0.5*dt*k2v,  k3v = accel(p + 0.5*dt*k2p, h2);
        vec3 k4p = v + dt*k3v,      k4v = accel(p + dt*k3p, h2);
        vec3 np = p + dt / 6.0 * (k1p + 2.0*k2p + 2.0*k3p + k4p);
        vec3 nv = v + dt / 6.0 * (k1v + 2.0*k2v + 2.0*k3v + k4v);

        vec3  mid = 0.5 * (p + np);
        float seg = length(np - p);
        float rm  = length(mid);
        minRadius = min(minRadius, rm);
        orbitSweep += seg * sqrt(max(h2, 0.0)) / max(rm * rm, 1e-4);

        float rr, surfaceDetail;
        // Do not miss the razor-thin equatorial disk when an adaptive step
        // crosses y=0. This makes primary and higher-order lensed images both
        // sharper and more stable at low sample counts.
        vec3 samplePos = mid;
        float sampleSeg = seg;
        if (p.y * np.y <= 0.0 && abs(p.y - np.y) > 1e-7){
            float crossT = clamp(p.y / (p.y - np.y), 0.0, 1.0);
            samplePos = mix(p, np, crossT);
            sampleSeg = max(seg, uDiskH * 0.65);
        }
        float dens = diskDensity(samplePos, rr, surfaceDetail);
        if (dens > 0.0){
            bool plunging = (rr < uDiskIn);

            float x    = uDiskIn / max(rr, RPLUNGE);
            float prof = pow(min(x, 1.0), 0.75)
                       * pow(max(1.0 - sqrt(min(x, 1.0)), 0.0), 0.25) / TNORM;
            prof = max(prof, 0.10);
            if (plunging) prof = mix(prof, 1.05, 0.65);

            // Doppler + redshift gravitationnel
            vec3  rot  = normalize(vec3(-samplePos.z, 0.0, samplePos.x));
            float beta = min(0.88, sqrt(M / max(rr, 1.15)));
            float gam  = inversesqrt(max(1.0 - beta*beta, 1e-4));
            vec3  n    = normalize(-v);
            float g    = sqrt(max(1.0 - RS / max(rr, RS*1.01), 0.0))
                       / (gam * (1.0 - beta * dot(rot, n)));
            g = mix(1.0, g, uDopp);
            g = clamp(g, 0.22, 2.85);

            // GRADIENT DE TEMPÉRATURE RADIAL : T ~ r^(-3/4)
            //  -> bord interne blanc-bleuté, extérieur ambre/cuivre.
            float tr   = pow(clamp(uDiskIn / max(rr, RPLUNGE), 0.0, 1.0), 0.75);
            float grad = mix(1.0, mix(0.55, 1.85, tr), clamp(uTempGrad, 0.0, 2.0) * 0.72);
            float T = uTemp * grad * pow(prof, 0.55) * mix(1.0, g, 0.90);
            if (plunging) T *= 0.78;

            float hotThread = smoothstep(0.34, 1.42, surfaceDetail);
            T *= mix(0.91, 1.09, hotThread);

            vec3  emi = blackbody(T) * pow(prof, 1.70) * pow(g, 3.0);
            emi *= mix(0.54, 1.34, hotThread);
            // bords externes ambre/cuivre (poussiere chaude) comme sur la ref
            float warmT = 1.0 - clamp(tr, 0.0, 1.0);
            emi *= mix(vec3(1.0), vec3(1.22, 0.78, 0.70), pow(warmT, 1.05));
            // Limb brightening of the optically thin atmosphere.  At the very
            // shallow reference angle it produces the brilliant horizontal
            // spine while retaining turbulent detail above and below it.
            float grazing = 1.0 - abs(dot(normalize(v), vec3(0.0, 1.0, 0.0)));
            emi *= 0.88 + 0.34 * pow(grazing, 3.0);
            if (plunging) emi *= 0.86;
#if USE_SHADOW
            emi *= mix(1.0, shadowUp(samplePos), 0.55);
#endif

            float k  = dens * uAbsorb;
            float tk = exp(-k * sampleSeg);
            float w  = (k > 1e-5) ? (1.0 - tk) / k : sampleSeg;
            col   += trans * emi * dens * w * uDiskBright * lamTint;
            trans *= tk;
        }

        // anneau de photons : rehaut fin sur la sphère de photons
        float dph = (rm - RPHOTON) / 0.115;
        float dph2 = (rm - RPHOTON * 1.06) / 0.34;
        ringAcc += (exp(-dph * dph) + 0.45 * exp(-dph2 * dph2)) * seg * trans;

        // couronne chaude, uniquement HORS de l'ombre
        if (rm > RS * 2.35){
            float d1 = (rm - 3.05) / 0.85;
            cor += trans * exp(-d1 * d1) * seg * uHalo * 0.10
                 * vec3(1.00, 0.72, 0.50);
        }

        p = np; v = nv;
    }

    // parametre d'impact : sert au voile lumineux qui baigne l'ombre
    float bImp = length(cross(uCamPos, rd));
    float bc   = 2.598;                    // 3*sqrt(3)*M  (rayon critique)
    float u    = clamp(bImp / bc, 0.0, 1.0);
    float nearCritical = exp(-pow((bImp / bc - 1.0) / 0.24, 2.0));
    float orbitGlow = smoothstep(1.7, 4.8, orbitSweep)
                    * exp(-pow((minRadius - RPHOTON) / 0.42, 2.0));
    float screenAngle = atan(uv.y, uv.x);
    float ringGrain = 0.92 + 0.08 * fbm(vec3(cos(screenAngle) * 3.1,
                                              sin(screenAngle) * 3.1,
                                              uTime * 0.035 + 17.0), 3);

    if (!eaten){
        col += cor;
        if (escaped) col += trans * (sceneColor + starfield(normalize(v)));
        else col += trans * sceneColor;
        col += vec3(1.00, 0.93, 0.95) * ringAcc * uRing * 0.038 * lamTint;
        // glare de lentille juste au-dessus du parametre critique : supprime
        // l'anneau sombre entre l'ombre et l'anneau de photons
        float critical = bImp / bc - 1.0;
        float narrowRing = exp(-pow(critical / 0.043, 2.0));
        float outerShoulder = exp(-pow((critical - 0.082) / 0.155, 2.0));
        col += vec3(1.00, 0.955, 0.975) * uRing
             * (0.115 * narrowRing * ringGrain
                + 0.025 * outerShoulder + 0.018 * orbitGlow);
        // Wide photographic scattering around the critical curve. The narrow
        // photon ring stays resolved while this layer produces the soft,
        // lavender Einstein halo visible in the reference image.
        col += vec3(1.00, 0.91, 0.96) * nearCritical * uHalo * 1.28;
    } else {
        // voile de diffusion / lumiere diffractee autour de l'horizon :
        // l'interieur de l'anneau reste lumineux (aspect photographique).
        // l'horizon ne rayonne pas : l'ombre reste quasi noire. Seule
        // subsiste la lumiere diffractee rasante juste sous le parametre
        // d'impact critique -> fin liseré, pas de disque blanc.
        float fill = 0.010 + 0.07 * pow(u, 5.0) + 1.55 * pow(u, 18.0);
        col += vec3(1.00, 0.955, 0.975) * uCore * fill * 0.72;
        col += cor * 0.02;
        col += vec3(1.00, 0.92, 0.88) * ringAcc * uRing * 0.012 * lamTint;
    }

    return col * BLACK_HOLE_BRIGHTNESS;
}

vec3 DrawEventHorizon(vec3 sceneColor, vec3 viewDirection,
                      mat4 viewToWorld, float animationTime,
                      float visibility) {
    if (visibility < 0.001) return sceneColor;

    vec3 worldRay = normalize(mat3(viewToWorld) * normalize(viewDirection));
    vec3 centerDirection = normalize(vec3(-0.56, 0.34, 0.76));
    float facing = dot(worldRay, centerDirection);
    if (facing < 0.72) return sceneColor;

    vec3 celestialRight = normalize(cross(vec3(0.0, 1.0, 0.0), centerDirection));
    vec3 celestialUp = normalize(cross(centerDirection, celestialRight));
    vec2 tangent = vec2(dot(worldRay, celestialRight), dot(worldRay, celestialUp))
                 / max(facing, 0.001);

    // Angular placement is an Iris-only camera adapter. The Overworld uses a
    // much smaller apparent scale so the same physical Python simulation reads
    // as a genuinely distant celestial object; The End keeps the hero scale.
    float angularScale = BLACK_HOLE_SIZE;
    #ifdef OVERWORLD
        angularScale *= BLACK_HOLE_OVERWORLD_SCALE;
    #endif
    tangent /= max(angularScale, 0.05);
    float bound = length(tangent);
    float objectMask = 1.0 - smoothstep(0.52, 0.64, bound);
    if (objectMask <= 0.0) return sceneColor;

    float rotation = BLACK_HOLE_ROTATION * 0.01745329251994;
    float cr = cos(rotation);
    float sr = sin(rotation);
    tangent = mat2(cr, -sr, sr, cr) * tangent;

    vec3 cameraForward = normalize(-BH_CAM_POS);
    vec3 cameraRight = normalize(cross(cameraForward, vec3(0.0, 1.0, 0.0)));
    vec3 cameraUp = normalize(cross(cameraRight, cameraForward));
    vec3 simulationRay = normalize(cameraForward
                         + cameraRight * tangent.x
                         + cameraUp * tangent.y);

    ehAnimationTime = animationTime * BLACK_HOLE_ANIMATION_SPEED;
    vec3 traced = EHTraceBlackHole(simulationRay, tangent, sceneColor,
                                   ehAnimationTime);
    return mix(sceneColor, traced,
               clamp(visibility, 0.0, 1.0) * objectMask);
}

#undef OCT_MAIN
#undef OCT_SUB
#undef OCT_PLUNGE
#undef OCT_NEB
#undef STAR_LAYERS
#undef USE_SHADOW
#undef USE_DISPERSION
#undef USE_VEIL
#undef DT_SCALE
#undef BH_STEPS
#undef uBhScale
#undef uDiskIn
#undef uDiskOut
#undef uDiskH
#undef uDiskBright
#undef uTemp
#undef uTempGrad
#undef uDopp
#undef uSpin
#undef uAbsorb
#undef uFil
#undef uPatch
#undef uConc
#undef uPlunge
#undef uDust
#undef uRing
#undef uCore
#undef uHalo
#undef uStars
#undef uSat
#undef uTime
#undef uSteps
#undef uCamPos
#undef RS
#undef M
#undef RPLUNGE
#undef RPHOTON
