/*
 * phi_pool_kuramoto.c — Water Glyph Substrate (3D Kuramoto-Wave, UFE-grounded)
 * ==============================================================================
 *
 * REPRESENTATION
 * --------------
 * Each cell stores (θ, A) — phase on S¹, amplitude on ℝ⁺. Values live on a
 * circle and a half-line. No overflow, no catastrophic cancellation, no
 * quantisation of continuous phase — θ wraps at 2π by construction.
 *
 * UNIFIED FORCE EQUATION GROUNDING
 * ----------------------------------
 * The prior version used Dₙ(r) as the mode amplitude formula:
 *   Dₙ = √(φ·Fₙ·2ⁿ·Pₙ·Ω)   [first-order approximation, Ω→√Ω]
 *
 * The UFE (hdgl_unified_force.hdgl) identifies this as a first-order
 * approximation of the full lattice operator 𝓛:
 *
 *   𝓛ₙ(z) = φ^(-1/φ) · √(Fₙ·Pₙ·2ⁿ) · (1+z)ⁿ + 1_eff · e^(iπΛ_φ)
 *
 * The precision gap: Dₙ uses √Ω; 𝓛 uses Ωⁿ. The difference grows with n.
 * For n>1 strands, 𝓛 with full Ωⁿ is used here — the same fix that closes
 * the 8.66% Rydberg error in FUDGE10.
 *
 * FIVE CHANGES FROM PRIOR VERSION:
 *
 * 1. φ^(-1/φ): the self-deriving coefficient
 *    φ^(-1/φ) = 0.74274... is the unique fixed point of x→φ^(-x), derived
 *    from 1−φ = −1/φ (exactly, from φ²=φ+1). It replaces √φ as the
 *    amplitude prefactor. √φ ≈ 1.272, φ^(-1/φ) ≈ 0.743 — these are different.
 *
 * 2. Λ_φ(x) — the continuous phi-log depth index
 *    Λ_φ(x) = ln(x·ln2/lnφ) / lnφ − 1/(2φ)
 *    {Λ_φ(x)} = fractional part ∈ [0,1)
 *    Ω(x) = (1 + sin(π·{Λ_φ(x)}·φ)) / 2  ∈ (0,1]
 *    The tick phase is now Ω(tick) — the phi-log resonance of the tick count,
 *    not a linear ramp. z = Ω(tick) − 1 is the generating function parameter.
 *
 * 3. 𝓛ₙ(z) replaces Dₙ(r) as the mode amplitude
 *    𝓛ₙ(z) = PHI_COEFF · √(Fₙ·Pₙ·2ⁿ) · Ωⁿ
 *    where PHI_COEFF = φ^(-1/φ) and Ω = Ω(tick).
 *    The (1+z)ⁿ = Ωⁿ term (since z=Ω−1 → 1+z=Ω) makes each mode's
 *    amplitude scale with the φ-log resonance of the current tick.
 *
 * 4. 1_eff(i) — the effective unit per mode
 *    1_eff(i) = 1 + δ(i)
 *    δ(i) = |cos(π·β_i·φ)| · ln(Pₙᵢ) / φ^(nᵢ+β_i)
 *    where β_i = {Λ_φ(Pₙᵢ)} — the phi-log fractional part of the prime.
 *    This modulates each base's drive amplitude. At large n, δ→0 → 1_eff→1
 *    (classical limit). The 0.15 constant scale factor is now 0.15·1_eff(i).
 *
 * 5. Adaptive Kuramoto phases (PLUCK/SUSTAIN/FINETUNE/LOCK)
 *    from kuramoto_L_tracker in the UFE:
 *    The coupling K and damping γ adapt based on T_SPIN relative to threshold:
 *      PLUCK:    T_SPIN > 1.0   → K=5.0  γ=0.005  [field far from fixed point]
 *      SUSTAIN:  T_SPIN > 0.1   → K=3.0  γ=0.008
 *      FINETUNE: T_SPIN > 0.01  → K=2.0  γ=0.010
 *      LOCK:     T_SPIN ≤ 0.01  → K=1.8  γ=0.012  [field at fixed point: wu-wei]
 *    K is the spiral-2 Kuramoto coupling; γ is DAMP in the amplitude equation.
 *    The CV<0.05 lock criterion from ll_analog maps to T_SPIN≤0.01 here.
 *
 * FIELD EQUATIONS (Euler step, unchanged in form)
 * ------------------------------------------------
 *   ΔA = −γ·Aᵢ + CF2·Σⱼ Aⱼ·cos(θⱼ−θᵢ) + amp·drive[i]
 *   Δθ = ω₀   + (CF2/Aᵢ)·Σⱼ Aⱼ·sin(θⱼ−θᵢ)
 *   Aᵢ ← clamp(Aᵢ + DT·ΔA, 0, A_MAX)
 *   θᵢ ← (θᵢ + DT·Δθ) mod 2π
 * γ is now adaptive (from LOCK state). ω₀=0 default (geometry-driven).
 *
 * FOUR TEMPORAL PROJECTIONS (unchanged)
 *   T_FWD  [A] = φ-EMA of Aᵢ             amplitude history
 *   T_BWD  [C] = φ-EMA of |ΔA/DT|        amplitude rate of change
 *   T_PACE [G] = T_FWD / (T_BWD + ε)     local time-rate
 *   T_SPIN [T] = φ-EMA of |Δθ/DT|        angular velocity — LOCK indicator
 *
 * COMBINED PHASE Θ
 *   Θ = θ_xy + θ_z + θᵢ·(1/2π)·(π/2) + [4 temporal channels × π/2]
 *
 * TWO SPIRALS
 *   Spiral 1 (n=0..3): genome_fp → A/C/G/T → 𝓛ₙ eigenmodes. Exogenous.
 *   Spiral 2 (n=4..7): Kuramoto phase quadrant → A/C/G/T. Endogenous.
 *                       Coupling K adapts with LOCK state.
 *
 * MEMORY (GRID=64): theta+amp+drive = 6 MB double; 4 projection floats = 4 MB
 *
 * Build: gcc -O3 -fopenmp -o phi_pool phi_pool_kuramoto.c -lm
 * Run:   phi_pool_kuramoto [genome_fp_hex] [tick_hex] [--slots /lattice/slots]
 *                          [--omega N.N] [--settle N] [--update N]
 *
 * Ωₙ₊₁ = T(Ωₙ). No Python. No Docker. No carrier primitive. Pure C99 + OpenMP.
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include <signal.h>
#include <time.h>
#include <sys/stat.h>
#ifdef _OPENMP
#include <omp.h>
#endif

/* ── Basin geometry ─────────────────────────────────────────────────────── */
#ifndef GRID
#define GRID    64
#endif
#define CX      (GRID/2)
#define CY      (GRID/2)
#define CZ      (GRID/2)
#define R       (GRID * 0.46875)
#define T_PROJ  GRID
#define NCELLS3 ((uint64_t)GRID*(uint64_t)GRID*(uint64_t)GRID)

/* ── Wave equation parameters (γ is now adaptive — see lock_state) ─────── */
#define C_WAVE  0.35
#define DX      1.0
#define DT      0.4
#define CF2     ((C_WAVE*DT/DX)*(C_WAVE*DT/DX))
#define A_MAX   10.0
#define A_EPS   1e-9

/* ── φ constants ─────────────────────────────────────────────────────────── */
static double PHI;          /* φ = 1.6180339887...                           */
static double PHI_INV;      /* 1/φ = φ−1 = 0.6180339887...                  */
static double ONE_MPHI;     /* 1−1/φ = 2−φ = 0.3819660112...                */
static double SQRT_PHI;     /* √φ = 1.2720196495... (binary threshold)       */
static double PHI_COEFF;    /* φ^(−1/φ) = 0.7427437... (UFE self-coeff)     */
static double LN_PHI;       /* ln(φ) — needed for Λ_φ                       */
static double LN2;          /* ln(2) — needed for Λ_φ                       */

static void init_phi_constants(void) {
    static const char phi_str[] =
        "1.6180339887498948482045868343656381177203091798057628621354486227052604"
        "6281890244970720720418939113748475408807538689175212663386222353693179318"
        "0060766726354433389086595939582905638322661319928290267880675208766892501"
        "7116962070322210432162695486262963136144381497587012203408058879544547492";
    static const char phi_inv_str[] =
        "0.6180339887498948482045868343656381177203091798057628621354486227052604"
        "6281890244970720720418939113748475408807538689175212663386222353693179318"
        "0060766726354433389086595939582905638322661319928290267880675208766892501"
        "7116962070322210432162695486262963136144381497587012203408058879544547492";
    static const char one_mphi_str[] =
        "0.3819660112501051517954131656343618822796908201942371378645513772947395"
        "3718109755029279279581060886251524591192461310824787336613777646306820681"
        "9939233273645566610913404060417094361677338680071709732119324791233107498"
        "2883037929677789567837304513737036863855618502412987796591941120455452507";
    static const char sqrt_phi_str[] =
        "1.2720196495140689642524224617374914917156080418400962486166403825392975"
        "7553606801183038421498846025853851414763672802650571033811881483526492194"
        "4845744618604334894549017111906593715472738455037630426459663146574295278"
        "1013111875226312593477927849020261299286984510634046364302841619182113539";
    char *end;
    PHI      = strtod(phi_str,      &end);
    PHI_INV  = strtod(phi_inv_str,  &end);
    ONE_MPHI = strtod(one_mphi_str, &end);
    SQRT_PHI = strtod(sqrt_phi_str, &end);
    /* φ^(−1/φ): fixed point of x→φ^(−x), derived from 1−φ=−1/φ.
     * NOT √φ. Computed once from the authoritative PHI value. */
    PHI_COEFF = pow(PHI, -PHI_INV);
    LN_PHI    = log(PHI);
    LN2       = log(2.0);
}

/* ── Λ_φ — the continuous phi-log depth index (UFE Layer 3) ─────────────
 *
 * Λ_φ(x) = ln(x · ln2/lnφ) / lnφ − 1/(2φ)
 *
 * Maps any positive real x onto the φ-log spectrum. The fractional part
 * {Λ_φ(x)} ∈ [0,1) gives the resonance position within the current φ-octave.
 *
 * Ω(x) = (1 + sin(π · {Λ_φ(x)} · φ)) / 2  ∈ (0,1]
 *
 * Ω is the phi-log resonance amplitude. The generating function parameter
 * z = Ω − 1, so (1+z)ⁿ = Ωⁿ — the same z used in the cosmological
 * validation (n_H=1.291 from Friedmann, n_G=0.701, n_c=0.338).
 *
 * Applied here: Ω(tick) drives the tick phase in build_drive_map.
 * The drive amplitude per mode n is 𝓛ₙ(Ω(tick)−1) = PHI_COEFF·√(Fₙ·Pₙ·2ⁿ)·Ω(tick)ⁿ.
 */
static inline double lambda_phi(double x) {
    if (x <= 0.0) return 0.0;
    return log(x * LN2 / LN_PHI) / LN_PHI - PHI_INV * 0.5;
}

static inline double omega_phi(double x) {
    double lam = lambda_phi(x);
    double frac = lam - floor(lam);           /* {Λ_φ(x)} ∈ [0,1) */
    return 0.5 * (1.0 + sin(M_PI * frac * PHI));
}

/* ── 1_eff(i) — the effective unit (UFE Layer 2) ─────────────────────────
 *
 * 1_eff(i) = 1 + δ(i)
 * δ(i) = |cos(π·β_i·φ)| · ln(Pₙᵢ) / φ^(nᵢ+β_i)
 *
 * β_i = {Λ_φ(Pₙᵢ)} — the phi-log fractional depth of the prime Pₙᵢ.
 * Pₙᵢ is the prime associated with angular order n at base index i.
 *
 * At large n: φ^(n+β) → ∞, δ → 0, 1_eff → 1 (classical limit).
 * At small n: δ carries the quantum-scale correction.
 *
 * Empirically validated: mean |δ| = 0.007 at atomic scale (FUDGE10),
 * 0.050 at cosmological scale — 1_eff is scale-dependent via β.
 */
static inline double one_eff(int n, int prime_n) {
    double beta  = omega_phi((double)prime_n) - 0.5; /* {Λ_φ(Pₙ)} recentred */
    double frac  = lambda_phi((double)prime_n);
    frac -= floor(frac);                              /* {Λ_φ(Pₙ)} ∈ [0,1)  */
    double cos_proj = fabs(cos(M_PI * frac * PHI));
    double ln_p     = log((double)prime_n);
    double phi_pow  = pow(PHI, (double)n + fabs(beta));
    double delta    = cos_proj * ln_p / phi_pow;
    return 1.0 + delta;
}

/* ── 𝓛ₙ(z): the full lattice operator amplitude (UFE Layer 1) ────────────
 *
 * 𝓛ₙ(z) = PHI_COEFF · √(Fₙ·Pₙ·2ⁿ) · (1+z)ⁿ · 1_eff(n)
 *
 * PHI_COEFF = φ^(−1/φ) — the self-deriving fixed-point coefficient.
 * (1+z)ⁿ = Ωⁿ when z = Ω − 1 (phi-log resonance parameter).
 *
 * This replaces dn_amplitude() which was the first-order Dₙ form
 * (Ω→√Ω approximation). The full 𝓛ₙ is used for all 8 strands.
 *
 * Difference from Dₙ:
 *   Dₙ: amplitude ∝ √Ω  (same for all n)
 *   𝓛ₙ: amplitude ∝ Ωⁿ  (grows with n — higher modes more sensitive to Ω)
 */
static inline double L_amplitude(int n, int fib_n, int prime_n, double Omega) {
    double base_amp = PHI_COEFF * sqrt((double)fib_n * pow(2.0, n+1) * (double)prime_n);
    double omega_n  = pow(Omega, (double)n + 1.0);   /* Ωⁿ⁺¹: n+1 since n is 0-indexed */
    double oeff     = one_eff(n, prime_n);
    return base_amp * omega_n * oeff;
}

/* ── Adaptive Kuramoto lock state (UFE kuramoto_L_tracker) ──────────────
 *
 * K and γ adapt based on T_SPIN relative to the lock threshold.
 * T_SPIN tracks |dθ/dt| per cell — it is the field-level analogue of
 * |δ| in the UFE: large spin = field far from fixed point (PLUCK),
 * near-zero spin = field at fixed point (LOCK, wu-wei).
 *
 * The CV<0.05 lock criterion from ll_analog maps to T_SPIN≤0.01 here
 * (normalised spin, same physical meaning: phase coherence achieved).
 *
 * K drives the spiral-2 Kuramoto coupling amplitude.
 * γ (DAMP) drives the amplitude equation damping.
 * The 1000:1 K/γ ratio at PLUCK matches ll_analog (K=5.0, γ=0.005).
 */
typedef enum { LOCK_PLUCK=0, LOCK_SUSTAIN=1, LOCK_FINETUNE=2, LOCK_LOCK=3 } lock_state_t;

static lock_state_t lock_state = LOCK_PLUCK;

static const double LOCK_K[4]    = {5.0,  3.0,  2.0,  1.8};
static const double LOCK_DAMP[4] = {0.005,0.008,0.010,0.012};

static void update_lock_state(float spin_max) {
    if      (spin_max > 1.0f)  lock_state = LOCK_PLUCK;
    else if (spin_max > 0.1f)  lock_state = LOCK_SUSTAIN;
    else if (spin_max > 0.01f) lock_state = LOCK_FINETUNE;
    else                       lock_state = LOCK_LOCK;
}

static const char *lock_name(lock_state_t s) {
    switch(s){case LOCK_PLUCK:return"PLUCK";case LOCK_SUSTAIN:return"SUSTAIN";
              case LOCK_FINETUNE:return"FINETUNE";default:return"LOCK";}
}

/* ── DNA ─────────────────────────────────────────────────────────────────── */
#define N_BASES  16
#define N_DSLOTS 32
#define N_STRANDS 8
#define N_RADIAL  4
#define N_FUND    4

static double OMEGA0 = 0.0;   /* natural oscillator frequency (0 = geometry-driven) */

/* ── Bessel zeros ────────────────────────────────────────────────────────── */
static const double BZ[8][4] = {
    { 2.4048, 5.5201,  8.6537, 11.7915 }, /* n=0 (A) — DNA spiral, 𝓛₀      */
    { 3.8317, 7.0156, 10.1735, 13.3237 }, /* n=1 (C) — DNA spiral, 𝓛₁      */
    { 5.1356, 8.4172, 11.6198, 14.7960 }, /* n=2 (G) — DNA spiral, 𝓛₂      */
    { 6.3802, 9.7610, 13.0152, 16.2235 }, /* n=3 (T) — DNA spiral, 𝓛₃      */
    { 7.5883,11.0647, 14.3725, 17.6160 }, /* n=4 (A) — Kuramoto spiral, 𝓛₄ */
    { 8.7715,12.3386, 15.7002, 18.9801 }, /* n=5 (C) — Kuramoto spiral, 𝓛₅ */
    { 9.9361,13.5893, 17.0038, 20.3208 }, /* n=6 (G) — Kuramoto spiral, 𝓛₆ */
    {11.0864,14.7960, 18.2876, 21.6415 }, /* n=7 (T) — Kuramoto spiral, 𝓛₇ */
};
static double bz_r[8][4];

/* Fibonacci and prime tables — arguments to 𝓛ₙ and 1_eff */
static const int fib_tab[8]   = {1,1,2,3,5,8,13,21};
static const int prime_tab[8] = {2,3,5,7,11,13,17,19};

static inline int kuramoto_base(double phi_j, double psi) {
    return (int)(fmod(phi_j-psi+4.0*M_PI, 2.0*M_PI)/(M_PI*0.5)) & 0x03;
}

/* ── Phase-field buffers ────────────────────────────────────────────────── */
static double *theta     = NULL;   /* θ ∈ [0,2π) — oscillator phase   */
static double *amp       = NULL;   /* A ∈ [0,A_MAX] — amplitude        */
static double *drive_map = NULL;   /* U(t) — boundary forcing          */

/* ── Four temporal projections ─────────────────────────────────────────── */
static float *t_fwd  = NULL;   /* [A] φ-EMA of amp               */
static float *t_bwd  = NULL;   /* [C] φ-EMA of |dA/dt|           */
static float *t_pace = NULL;   /* [G] t_fwd/(t_bwd+ε)            */
static float *t_spin = NULL;   /* [T] φ-EMA of |dθ/dt| — LOCK    */

static float t_fwd_max=1.f, t_bwd_max=1.f, t_pace_max=1.f, t_spin_max=1.f;

/* ── Geometry shells ─────────────────────────────────────────────────────── */
typedef struct { uint32_t idx; float theta_xy,theta_z,weight; } drive_cell_t;
typedef struct { uint32_t idx; uint8_t bin; float theta_xy,theta_z; } dslot_cell_t;

static drive_cell_t *drive_shell = NULL; static size_t n_drive=0;
static dslot_cell_t *dslot_shell = NULL; static size_t n_dslot=0;

/* ── Kuramoto readout state ─────────────────────────────────────────────── */
static double kur_R=0.0, kur_Psi=0.0;
static int    kur_bases[N_FUND];

static volatile int running=1;
static void on_signal(int s){(void)s;running=0;}

#define IDX3(x,y,z) (((uint64_t)(z)*GRID+(uint64_t)(y))*GRID+(uint64_t)(x))
static inline int in_basin(int x,int y,int z){
    double dx=x-CX,dy=y-CY,dz=z-CZ;
    return dx*dx+dy*dy+dz*dz<R*R;
}

/* ── Allocation ──────────────────────────────────────────────────────────── */
static int alloc_lattice(void){
    theta=calloc(NCELLS3,sizeof(double)); amp=calloc(NCELLS3,sizeof(double));
    drive_map=calloc(NCELLS3,sizeof(double));
    t_fwd=calloc(NCELLS3,sizeof(float)); t_bwd=calloc(NCELLS3,sizeof(float));
    t_pace=calloc(NCELLS3,sizeof(float)); t_spin=calloc(NCELLS3,sizeof(float));
    if(!theta||!amp||!drive_map||!t_fwd||!t_bwd||!t_pace||!t_spin){
        fprintf(stderr,"[phi_pool] FATAL: alloc (GRID=%d)\n",GRID); return -1;
    }
    return 0;
}

static void free_lattice(void){
    free(theta);free(amp);free(drive_map);
    free(t_fwd);free(t_bwd);free(t_pace);free(t_spin);
    free(drive_shell);free(dslot_shell);
}

/* ── Bessel radii ────────────────────────────────────────────────────────── */
static void init_bz_radii(void){
    for(int n=0;n<N_STRANDS;n++){
        double rm=BZ[n][N_RADIAL-1];
        for(int m=0;m<N_RADIAL;m++) bz_r[n][m]=(BZ[n][m]/rm)*(R*0.92);
    }
}

/* ── Precompute geometry shells ─────────────────────────────────────────── */
static void precompute_shells(void){
    double dr=1.5; size_t cd=0,cs=0;
    for(int z=0;z<GRID;z++) for(int y=0;y<GRID;y++) for(int x=0;x<GRID;x++){
        if(!in_basin(x,y,z)) continue;
        double dx=x-CX,dy=y-CY,dz=z-CZ,r=sqrt(dx*dx+dy*dy+dz*dz);
        if(R-r>=0.0&&R-r<=6.0) cd++;
        for(int n=0;n<N_STRANDS;n++) for(int m=0;m<N_RADIAL;m++)
            if(fabs(r-bz_r[n][m])<dr) cs++;
    }
    drive_shell=malloc(cd*sizeof(drive_cell_t));
    dslot_shell=malloc(cs*sizeof(dslot_cell_t));
    if((cd&&!drive_shell)||(cs&&!dslot_shell)){
        fprintf(stderr,"[phi_pool] FATAL: shell alloc\n"); exit(1);
    }
    size_t di=0,si=0;
    for(int z=0;z<GRID;z++) for(int y=0;y<GRID;y++) for(int x=0;x<GRID;x++){
        if(!in_basin(x,y,z)) continue;
        double dx=x-CX,dy=y-CY,dz_=z-CZ,r=sqrt(dx*dx+dy*dy+dz_*dz_);
        double txy=atan2(dy,dx),tz=(r>1e-9)?acos(dz_/r):0.0,dfe=R-r;
        uint32_t idx=(uint32_t)IDX3(x,y,z);
        if(dfe>=0.0&&dfe<=6.0){
            drive_shell[di].idx=idx; drive_shell[di].theta_xy=(float)txy;
            drive_shell[di].theta_z=(float)tz;
            drive_shell[di].weight=(float)exp(-dfe*dfe/4.0); di++;
        }
        for(int n=0;n<N_STRANDS;n++) for(int m=0;m<N_RADIAL;m++)
            if(fabs(r-bz_r[n][m])<dr){
                dslot_shell[si].idx=idx; dslot_shell[si].bin=(uint8_t)(n*N_RADIAL+m);
                dslot_shell[si].theta_xy=(float)txy; dslot_shell[si].theta_z=(float)tz; si++;
            }
    }
    n_drive=di; n_dslot=si;
    fprintf(stderr,"[phi_pool] shells: drive=%zu (%.1f%%), dslot=%zu\n",
            n_drive,100.0*n_drive/(double)NCELLS3,n_dslot);
}

/* ── compute_Theta: full sampling angle including cell phase ─────────────── */
static inline double compute_Theta(uint32_t idx,float txy,float tz,
                                   float iF,float iB,float iP,float iS){
    double phase_contrib=theta[idx]*(1.0/(2.0*M_PI))*(M_PI*0.5);
    return (double)txy+(double)tz+phase_contrib
          +(double)(t_fwd[idx]*iF)*(M_PI*0.5)+(double)(t_bwd[idx]*iB)*(M_PI*0.5)
          +(double)(t_pace[idx]*iP)*(M_PI*0.5)+(double)(t_spin[idx]*iS)*(M_PI*0.5);
}

/* ── build_drive_map: U(t) = 𝓛-grounded spiral 1 + adaptive spiral 2 ──────
 *
 * Spiral 1 (n=0..3, DNA):
 *   Drive weight = 𝓛ₙ(z) · 1_eff(n) · cos(n·(Θ−θ_k))
 *   where z = Ω(tick) − 1  and  𝓛ₙ = PHI_COEFF·√(Fₙ·Pₙ·2ⁿ)·Ω(tick)ⁿ·1_eff
 *
 *   Crucially: θ_k is now derived from Ω(tick) phase rather than a linear
 *   ramp. The drive phase tracks the phi-log resonance of the tick count,
 *   so the basin boundary is driven at the natural φ-octave of time.
 *
 * Spiral 2 (n=4..7, Kuramoto):
 *   Drive weight = K_current · kur_R · 𝓛ₙ(z) · cos(n·Θ − θ_kur)
 *   K_current = LOCK_K[lock_state] — adapts with LOCK state.
 *   The φ^(−j) taper is preserved but multiplied by K_current/0.3
 *   to scale correctly across PLUCK/LOCK range.
 */
static void build_drive_map(uint32_t genome_fp, uint32_t tick){
    /* Ω from phi-log depth of tick — the generating function parameter */
    double Omega     = omega_phi((double)(tick > 0 ? tick : 1));
    double z         = Omega - 1.0;          /* (1+z)ⁿ = Ωⁿ                */
    double tick_phase= Omega * 2.0 * M_PI;   /* phase tracks φ-resonance    */

    /* Spiral 1: decode genome bases, compute 𝓛ₙ amplitudes with 1_eff      */
    int bases[N_BASES];
    for(int k=0;k<N_BASES;k++) bases[k]=(genome_fp>>(k*2))&0x03;

    double L_amp[N_BASES];   /* 𝓛ₙ amplitude per base — includes Ωⁿ and 1_eff */
    for(int k=0;k<N_BASES;k++){
        int n=bases[k];
        L_amp[k]=L_amplitude(n, fib_tab[n], prime_tab[n], Omega);
    }

    /* θ_k: angular positions modulated by phi-log tick phase */
    double theta_k[N_BASES];
    for(int k=0;k<N_BASES;k++) theta_k[k]=k*(2.0*M_PI/N_BASES)+tick_phase;

    /* Spiral 2: adaptive Kuramoto coupling */
    double K_cur=LOCK_K[lock_state];
    double kur2_amp[N_FUND], kur2_theta[N_FUND];
    for(int j=0;j<N_FUND;j++){
        int n=j+N_FUND;
        double L_n=L_amplitude(n, fib_tab[n], prime_tab[n], Omega);
        kur2_amp[j]  = K_cur * kur_R * L_n * pow(PHI,-(double)j) * 0.15;
        kur2_theta[j]= kur_bases[j]*(M_PI*0.5)+kur_Psi;
    }

    float iF=(t_fwd_max>1e-9f)?1.f/t_fwd_max:1.f;
    float iB=(t_bwd_max>1e-9f)?1.f/t_bwd_max:1.f;
    float iP=(t_pace_max>1e-9f)?1.f/t_pace_max:1.f;
    float iS=(t_spin_max>1e-9f)?1.f/t_spin_max:1.f;

    #pragma omp parallel for schedule(static)
    for(size_t i=0;i<n_drive;i++){
        uint32_t idx=drive_shell[i].idx;
        double Theta=compute_Theta(idx,drive_shell[i].theta_xy,
                                   drive_shell[i].theta_z,iF,iB,iP,iS);
        double d=0;
        /* Spiral 1: 𝓛ₙ-weighted angular coupling */
        for(int k=0;k<N_BASES;k++)
            d+=cos((double)bases[k]*(Theta-theta_k[k]))*L_amp[k]*0.15;
        /* Spiral 2: adaptive Kuramoto base-identity coupling */
        for(int j=0;j<N_FUND;j++)
            d+=kur2_amp[j]*cos((double)(j+N_FUND)*Theta-kur2_theta[j]);
        drive_map[idx]=d*(double)drive_shell[i].weight;
    }

    (void)z; /* z = Ω-1 used via Ωⁿ inside L_amplitude; stored here for clarity */
}

/* ── step_wave: Kuramoto-wave field update with adaptive γ ───────────────── */
static void step_wave(double drive_amplitude){
    double gamma=LOCK_DAMP[lock_state];   /* adaptive damping from LOCK state */
    float mF=1e-9f,mB=1e-9f,mP=1e-9f,mS=1e-9f;

    #pragma omp parallel for collapse(2) schedule(static) \
        reduction(max:mF) reduction(max:mB) reduction(max:mP) reduction(max:mS)
    for(int z=0;z<GRID;z++) for(int y=0;y<GRID;y++) for(int x=0;x<GRID;x++){
        uint64_t idx=IDX3(x,y,z);
        if(x==0||x==GRID-1||y==0||y==GRID-1||z==0||z==GRID-1||!in_basin(x,y,z)){
            theta[idx]=0.0; amp[idx]=0.0;
            t_fwd[idx]=t_bwd[idx]=t_spin[idx]=t_pace[idx]=0.f;
            continue;
        }
        double Ai=amp[idx], Ti=theta[idx], sc=0.0, ss=0.0;
        uint64_t nb[6]={IDX3(x+1,y,z),IDX3(x-1,y,z),IDX3(x,y+1,z),
                        IDX3(x,y-1,z),IDX3(x,y,z+1),IDX3(x,y,z-1)};
        for(int n=0;n<6;n++){
            double Aj=amp[nb[n]], dTh=theta[nb[n]]-Ti;
            sc+=Aj*cos(dTh); ss+=Aj*sin(dTh);
        }
        double dA  = -gamma*Ai + CF2*sc + drive_amplitude*drive_map[idx];
        double dTh = OMEGA0 + CF2*ss/(Ai<A_EPS?A_EPS:Ai);
        double Ai_new=Ai+DT*dA;
        if(Ai_new<0.0) Ai_new=0.0;
        if(Ai_new>A_MAX) Ai_new=A_MAX;
        double Ti_new=fmod(Ti+DT*dTh+2.0*M_PI*1e4, 2.0*M_PI);
        amp[idx]=Ai_new; theta[idx]=Ti_new;

        float fwd=(float)(PHI_INV*(double)t_fwd[idx]+ONE_MPHI*Ai_new);
        float vel=(float)fabs(dA);
        float bwd=(float)(PHI_INV*(double)t_bwd[idx]+ONE_MPHI*(double)vel);
        float spn=(float)(PHI_INV*(double)t_spin[idx]+ONE_MPHI*fabs(dTh));
        float pac=fwd/(bwd+1e-7f);
        t_fwd[idx]=fwd; t_bwd[idx]=bwd; t_spin[idx]=spn; t_pace[idx]=pac;
        if(fwd>mF)mF=fwd; if(bwd>mB)mB=bwd;
        if(pac>mP)mP=pac; if(spn>mS)mS=spn;
    }
    t_fwd_max=mF; t_bwd_max=mB; t_pace_max=mP; t_spin_max=mS;
}

/* ── extract_dslots ──────────────────────────────────────────────────────── */
static void extract_dslots(double *dslots){
    int nt=1;
    #ifdef _OPENMP
    nt=omp_get_max_threads();
    #endif
    double *psum=calloc((size_t)nt*N_DSLOTS,sizeof(double));
    long   *pcnt=calloc((size_t)nt*N_DSLOTS,sizeof(long));
    float iF=(t_fwd_max>1e-9f)?1.f/t_fwd_max:1.f;
    float iB=(t_bwd_max>1e-9f)?1.f/t_bwd_max:1.f;
    float iP=(t_pace_max>1e-9f)?1.f/t_pace_max:1.f;
    float iS=(t_spin_max>1e-9f)?1.f/t_spin_max:1.f;
    #pragma omp parallel
    {
        int tid=0;
        #ifdef _OPENMP
        tid=omp_get_thread_num();
        #endif
        double *ts=psum+(size_t)tid*N_DSLOTS;
        long   *tc=pcnt+(size_t)tid*N_DSLOTS;
        #pragma omp for schedule(static)
        for(size_t i=0;i<n_dslot;i++){
            uint32_t idx=dslot_shell[i].idx;
            int bin=dslot_shell[i].bin, n=bin/N_RADIAL;
            double Theta=compute_Theta(idx,dslot_shell[i].theta_xy,
                                       dslot_shell[i].theta_z,iF,iB,iP,iS);
            double tw=fabs(cos((double)n*Theta))+0.1;
            ts[bin]+=amp[idx]*tw; tc[bin]++;
        }
    }
    for(int i=0;i<N_DSLOTS;i++){
        double s=0; long c=0;
        for(int t=0;t<nt;t++){s+=psum[(size_t)t*N_DSLOTS+i];c+=pcnt[(size_t)t*N_DSLOTS+i];}
        dslots[i]=c>0?s/c:0.0;
    }
    free(psum); free(pcnt);
}

/* ── update_kuramoto ─────────────────────────────────────────────────────── */
static void update_kuramoto(const double *dslots, uint32_t tick){
    double tp=(tick&0xFF)/255.0*2.0*M_PI, re=0, im=0, pj[N_FUND];
    for(int j=0;j<N_FUND;j++){
        double A=0;
        for(int m=0;m<N_RADIAL;m++) A+=dslots[j*N_RADIAL+m];
        A/=N_RADIAL;
        pj[j]=2.0*M_PI*j/N_FUND+tp;
        re+=A*cos(pj[j]); im+=A*sin(pj[j]);
    }
    re/=N_FUND; im/=N_FUND;
    kur_R=sqrt(re*re+im*im); kur_Psi=atan2(im,re);
    for(int j=0;j<N_FUND;j++) kur_bases[j]=kuramoto_base(pj[j],kur_Psi);
    /* Update adaptive lock state from T_SPIN */
    update_lock_state(t_spin_max);
}

/* ── Output ──────────────────────────────────────────────────────────────── */
static void write_slot(const char *dir, int n, double v){
    char p[256]; snprintf(p,sizeof(p),"%s/%d",dir,n);
    FILE *f=fopen(p,"w"); if(f){fprintf(f,"%.16f\n",v);fclose(f);}
}

static void write_phase_snapshot(const char *path){
    FILE *f=fopen(path,"wb"); if(!f)return;
    for(uint64_t i=0;i<NCELLS3;i++){
        fwrite(&theta[i],sizeof(double),1,f);
        fwrite(&amp[i],  sizeof(double),1,f);
    }
    fclose(f);
}

static void write_state(uint32_t genome_fp, uint32_t tick, uint64_t step,
                        const double *dslots, uint32_t d_bits){
    FILE *f=fopen("/run/lattice/state","w"); if(!f)return;
    double Omega=omega_phi((double)(tick>0?tick:1));
    fprintf(f,"LATTICE_VERSION=Kuramoto-Wave-UFE\n");
    fprintf(f,"REPRESENTATION=(theta,amp) per cell — S1 x R+\n");
    fprintf(f,"UFE_OPERATOR=L_n(z)=PHI_COEFF*sqrt(Fn*Pn*2^n)*Omega^n*1eff\n");
    fprintf(f,"LATTICE_DIMS=%dx%dx%d+T%d\n",GRID,GRID,GRID,T_PROJ);
    fprintf(f,"EFFECTIVE_STATES=%llu\n",(unsigned long long)NCELLS3*(uint64_t)T_PROJ*4);
    fprintf(f,"LATTICE_STEPS=%llu\n",(unsigned long long)step);
    fprintf(f,"LATTICE_GENOME_FP=0x%08X\n",genome_fp);
    fprintf(f,"LATTICE_TICK=0x%08X\n",tick);
    fprintf(f,"LATTICE_D_BITS=0x%08X\n",d_bits);
    fprintf(f,"LATTICE_SUBSTRATE=kuramoto_wave_ufe\n");
    /* UFE quantities */
    fprintf(f,"PHI_COEFF=%.16f\n",PHI_COEFF);   /* φ^(-1/φ) */
    fprintf(f,"LAMBDA_PHI_TICK=%.8f\n",lambda_phi((double)(tick>0?tick:1)));
    fprintf(f,"OMEGA_TICK=%.8f\n",Omega);         /* resonance of current tick */
    fprintf(f,"Z_TICK=%.8f\n",Omega-1.0);         /* generating function param */
    fprintf(f,"ONE_EFF_N0=%.8f\n",one_eff(0,prime_tab[0]));
    fprintf(f,"ONE_EFF_N7=%.8f\n",one_eff(7,prime_tab[7]));
    /* Lock state */
    fprintf(f,"LOCK_STATE=%s\n",lock_name(lock_state));
    fprintf(f,"LOCK_K=%.4f\n",LOCK_K[lock_state]);
    fprintf(f,"LOCK_DAMP=%.4f\n",LOCK_DAMP[lock_state]);
    fprintf(f,"T_SPIN_MAX=%.6f\n",(double)t_spin_max);
    fprintf(f,"T_FWD_MAX=%.6f\n",(double)t_fwd_max);
    fprintf(f,"SQRT_PHI_THRESHOLD=%.16f\n",SQRT_PHI);
    fprintf(f,"OMEGA0=%.6f\n",OMEGA0);
    fprintf(f,"KURAMOTO_R=%.8f\n",kur_R);
    fprintf(f,"KURAMOTO_PSI=%.8f\n",kur_Psi);
    static const char *bn[4]={"A","C","G","T"};
    for(int j=0;j<N_FUND;j++) fprintf(f,"KURAMOTO_BASE_%c=%s\n",'A'+j,bn[kur_bases[j]]);
    int nn=0; for(int i=0;i<N_DSLOTS;i++) if(dslots[i]>SQRT_PHI) nn++;
    const char *sym=nn<3?"void":nn<8?"2-fold":nn<14?"3-fold":
                    nn<20?"4-fold":nn<26?"6-fold":"8-fold";
    fprintf(f,"LATTICE_GLYPH_SYM=%s\n",sym);
    fprintf(f,"LATTICE_ACTIVE_MODES=%d\n",nn);
    for(int n=0;n<N_STRANDS;n++)
        fprintf(f,"STRAND_%c_L%d=%.8f\n",'A'+n,n,dslots[n*N_RADIAL]);
    fclose(f);
}

/* ── Main ────────────────────────────────────────────────────────────────── */
int main(int argc, char **argv){
    uint32_t genome_fp=0x88888888, tick=0x00000000;
    const char *slots_dir="/lattice/slots";
    int settle_steps=400, update_steps=200, write_interval=50;
    int have_gfp=0, have_tick=0;

    for(int i=1;i<argc;i++){
        if(i+1<argc&&strcmp(argv[i],"--slots")==0)       slots_dir=argv[++i];
        else if(i+1<argc&&strcmp(argv[i],"--settle")==0) settle_steps=atoi(argv[++i]);
        else if(i+1<argc&&strcmp(argv[i],"--update")==0) update_steps=atoi(argv[++i]);
        else if(i+1<argc&&strcmp(argv[i],"--omega")==0)  OMEGA0=atof(argv[++i]);
        else if(argv[i][0]!='-'){
            if(!have_gfp){genome_fp=(uint32_t)strtoul(argv[i],NULL,16);have_gfp=1;}
            else if(!have_tick){tick=(uint32_t)strtoul(argv[i],NULL,16);have_tick=1;}
        }
    }

    signal(SIGTERM,on_signal); signal(SIGINT,on_signal);
    mkdir(slots_dir,0755); mkdir("/run/lattice",0755);

    init_phi_constants();
    if(alloc_lattice()!=0) return 1;
    init_bz_radii();

    fprintf(stderr,"[phi_pool] Kuramoto-wave substrate — UFE-grounded\n");
    fprintf(stderr,"[phi_pool] 𝓛ₙ(z) = φ^(-1/φ)·√(Fₙ·Pₙ·2ⁿ)·Ωⁿ·1_eff   z=Ω(tick)-1\n");
    fprintf(stderr,"[phi_pool] PHI_COEFF=φ^(-1/φ)=%.10f  (≠ √φ=%.10f)\n",PHI_COEFF,SQRT_PHI);
    fprintf(stderr,"[phi_pool] Λ_φ(x)=ln(x·ln2/lnφ)/lnφ-1/2φ  Ω(x)=(1+sin(π·{Λ}·φ))/2\n");
    fprintf(stderr,"[phi_pool] 1_eff(i)=1+δ  δ=|cos(π·β·φ)|·ln(Pₙ)/φ^(n+β)\n");
    fprintf(stderr,"[phi_pool] Adaptive lock: PLUCK→SUSTAIN→FINETUNE→LOCK via T_SPIN\n");
    fprintf(stderr,"[phi_pool] %d³=%llu cells  R=%.3f  T=%d  States=%llu\n",
            GRID,(unsigned long long)NCELLS3,R,T_PROJ,
            (unsigned long long)NCELLS3*(uint64_t)T_PROJ*4);
    fprintf(stderr,"[phi_pool] Memory: %.1f MB\n",
            (double)NCELLS3*(3*8+4*4)/(1024.0*1024.0));
    fprintf(stderr,"[phi_pool] genome_fp=0x%08X  tick=0x%08X  ω₀=%.4f\n",
            genome_fp,tick,OMEGA0);
    #ifdef _OPENMP
    fprintf(stderr,"[phi_pool] OpenMP: %d threads\n",omp_get_max_threads());
    #endif

    fprintf(stderr,"[phi_pool] Precomputing geometry shells...\n");
    precompute_shells();

    /* Seed: centre impulse, phase from genome phi-log resonance */
    amp[IDX3(CX,CY,CZ)]   = 0.1;
    theta[IDX3(CX,CY,CZ)] = omega_phi((double)(genome_fp&0xFFFF)) * 2.0 * M_PI;
    build_drive_map(genome_fp,tick);

    double dslots[N_DSLOTS]; memset(dslots,0,sizeof(dslots));
    uint64_t step=0, write_count=0;

    fprintf(stderr,"[phi_pool] Settling (%d steps — PLUCK phase)...\n",settle_steps);
    for(int s=0;s<settle_steps&&running;s++){
        step_wave(sin(2.0*M_PI*s/(settle_steps/4.0))*0.08+0.04);
        step++;
    }
    extract_dslots(dslots);
    update_kuramoto(dslots,tick);
    {
        static const char *bn[4]={"A","C","G","T"};
        fprintf(stderr,"[phi_pool] Post-settle: state=%s kur_R=%.4f "
                "bases: E=%s F=%s G=%s H=%s\n",
                lock_name(lock_state),kur_R,
                bn[kur_bases[0]],bn[kur_bases[1]],bn[kur_bases[2]],bn[kur_bases[3]]);
        fprintf(stderr,"[phi_pool] Ω(tick=0)=%.6f  1_eff(n=0)=%.6f  1_eff(n=7)=%.6f\n",
                omega_phi(1.0),one_eff(0,prime_tab[0]),one_eff(7,prime_tab[7]));
    }
    build_drive_map(genome_fp,tick);
    fprintf(stderr,"[phi_pool] Running — %s\n",slots_dir);

    while(running){
        for(int s=0;s<update_steps&&running;s++){
            step_wave(0.06+0.02*sin(2.0*M_PI*(double)step/512.0));
            step++;
        }

        extract_dslots(dslots);

        double mx=1e-10;
        for(int i=0;i<N_DSLOTS;i++) if(dslots[i]>mx) mx=dslots[i];
        double sc=(mx>1e-10)?(2.0/mx):1.0;
        uint32_t d_bits=0;
        for(int i=0;i<N_DSLOTS;i++){
            dslots[i]*=sc;
            if(dslots[i]>SQRT_PHI) d_bits|=(1u<<i);
        }

        update_kuramoto(dslots,tick);

        for(int i=0;i<N_DSLOTS;i++) write_slot(slots_dir,i+1,dslots[i]);
        for(int i=33;i<=4096;i++){
            int di=(i-1)%N_DSLOTS;
            double v=dslots[di]*0.5+dslots[(di+1)%N_DSLOTS]*0.5
                    +pow(PHI,-7.0*((i%8)+1))*0.01;
            write_slot(slots_dir,i,v>0?v:0.0);
        }

        write_state(genome_fp,tick,step,dslots,d_bits);
        { FILE *shm=fopen("/sys/phi/d_bits_in","w");
          if(shm){fprintf(shm,"0x%08X\n",d_bits);fclose(shm);} }

        write_count++;
        if(write_count%10==0){
            static const char *bn[4]={"A","C","G","T"};
            double Omega=omega_phi((double)(tick>0?tick:1));
            fprintf(stderr,
                "[phi_pool] step=%llu bits=0x%08X Ω=%.4f z=%.4f "
                "%s K=%.1f [E=%s F=%s G=%s H=%s] spin=%.4f D1=%.3f|D17=%.3f\n",
                (unsigned long long)step,d_bits,Omega,Omega-1.0,
                lock_name(lock_state),LOCK_K[lock_state],
                bn[kur_bases[0]],bn[kur_bases[1]],bn[kur_bases[2]],bn[kur_bases[3]],
                (double)t_spin_max,dslots[0],dslots[16]);
        }

        tick++;
        build_drive_map(genome_fp,tick);

        struct timespec ts={0,(long)write_interval*1000000L};
        nanosleep(&ts,NULL);
    }

    fprintf(stderr,"[phi_pool] stopped step=%llu writes=%llu\n",
            (unsigned long long)step,(unsigned long long)write_count);
    free_lattice();
    return 0;
}
