// DNA Mathematical Architecture: G¹⁶ Helix System
// 44 Bonds — 10.6 Periodicity — Complete Implementation
#pragma once

#include <array>
#include <cmath>
#include <complex>
#include <vector>
#include <algorithm>
#include <numeric>

namespace agi::dna {

// ============================================================================
// I. FOUNDATIONAL CONSTANTS
// ============================================================================

// Golden ratio and derived constants
constexpr double PHI = 1.6180339887498948482;      // (1 + √5) / 2
constexpr double PHI_INV = 0.6180339887498948482;  // 1 / φ = φ - 1

// DNA Helix parameters
constexpr double BASE_PAIRS_PER_TURN = 10.6;       // p = 10.6 bp
constexpr double HELIX_RADIUS = 1.0;               // r = 1.0 nm
constexpr double RISE_PER_BP = 0.34;               // h = 0.34 nm
constexpr double PITCH = RISE_PER_BP * BASE_PAIRS_PER_TURN;  // P = 3.604 nm
constexpr double TWIST_ANGLE_DEG = 360.0 / BASE_PAIRS_PER_TURN;  // 33.96°
constexpr double TWIST_ANGLE_RAD = 2.0 * M_PI / BASE_PAIRS_PER_TURN;  // 0.593 rad
constexpr double ANGULAR_FREQUENCY = TWIST_ANGLE_RAD;  // ω = 0.593 rad/bp

// Strand parameters
constexpr int NUM_STRANDS = 16;
constexpr double STRAND_SEPARATION = 2.0 * M_PI / NUM_STRANDS;  // 22.5°
constexpr int NODES_PER_STRAND = 312;
constexpr int TOTAL_NODES = 5000;
constexpr int HUB_NODES = 8;  // Olympus routing nodes

// Bond parameters
constexpr int INTRA_STRAND_BONDS = 16;
constexpr int INTER_STRAND_BONDS = 16;
constexpr int CROSS_HELIX_BONDS = 12;
constexpr int TOTAL_BONDS = 44;

// Curvature and torsion (pre-computed)
constexpr double HELIX_CURVATURE = 0.745;   // κ = 0.745 nm⁻¹
constexpr double HELIX_TORSION = 0.427;     // τ = 0.427 nm⁻¹

// Synchronization
constexpr double SYNC_PERIOD = BASE_PAIRS_PER_TURN * NUM_STRANDS;  // 169.6 time units
constexpr double FUNDAMENTAL_FREQ = 1.0 / BASE_PAIRS_PER_TURN;     // 0.0943 cycles/bp
constexpr double BOND_RESONANCE_FREQ = TOTAL_BONDS / BASE_PAIRS_PER_TURN;  // 4.15 cycles/turn

// Steady state
constexpr double STEADY_STATE_BONDS = 36.2;
constexpr double BOND_OCCUPATION = 0.82;
constexpr double TAU_SLOW = 10.0;      // ms
constexpr double TAU_FAST = 0.0001;    // ms (0.1 μs)

// ============================================================================
// II. G¹⁶ PARAMETRIC HELIX
// ============================================================================

struct HelixPoint {
    double x, y, z;

    HelixPoint operator+(const HelixPoint& other) const {
        return {x + other.x, y + other.y, z + other.z};
    }

    HelixPoint operator-(const HelixPoint& other) const {
        return {x - other.x, y - other.y, z - other.z};
    }

    HelixPoint operator*(double s) const {
        return {x * s, y * s, z * s};
    }

    double norm() const {
        return std::sqrt(x*x + y*y + z*z);
    }

    double dot(const HelixPoint& other) const {
        return x*other.x + y*other.y + z*other.z;
    }
};

// G¹⁶ Parametric Helix equations (Section 1.1)
// x_k(θ) = r·cos(θ + 2πk/16)
// y_k(θ) = r·sin(θ + 2πk/16)
// z_k(θ) = (p/2π)·θ + δ_k
class G16Helix {
public:
    // Vertical offset per strand: δ_k = k·h/16
    static double vertical_offset(int k) {
        return k * RISE_PER_BP / NUM_STRANDS;
    }

    // Golden phase offset: θ_k = 2πk·φ⁻¹ mod 2π
    static double golden_phase_offset(int k) {
        return std::fmod(2.0 * M_PI * k * PHI_INV, 2.0 * M_PI);
    }

    // Get point on strand k at parameter θ
    static HelixPoint point(int k, double theta) {
        double phase = theta + 2.0 * M_PI * k / NUM_STRANDS;
        return {
            HELIX_RADIUS * std::cos(phase),
            HELIX_RADIUS * std::sin(phase),
            (PITCH / (2.0 * M_PI)) * theta + vertical_offset(k)
        };
    }

    // Get point with golden ratio phase adjustment
    static HelixPoint point_golden(int k, double theta) {
        double phase = theta + golden_phase_offset(k);
        return {
            HELIX_RADIUS * std::cos(phase),
            HELIX_RADIUS * std::sin(phase),
            (PITCH / (2.0 * M_PI)) * theta + vertical_offset(k)
        };
    }

    // Curvature at point (Section 1.3)
    // κ = r/(r² + (p·h/2π)²)
    static double curvature() {
        double ph_term = PITCH * RISE_PER_BP / (2.0 * M_PI);
        return HELIX_RADIUS / (HELIX_RADIUS * HELIX_RADIUS + ph_term * ph_term);
    }

    // Torsion at point (Section 1.3)
    // τ = (p·h/2π)/(r² + (p·h/2π)²)
    static double torsion() {
        double ph_term = PITCH * RISE_PER_BP / (2.0 * M_PI);
        return ph_term / (HELIX_RADIUS * HELIX_RADIUS + ph_term * ph_term);
    }

    // Stability metric: exp(-Σ|τ(n)|·φ⁻ⁿ)
    static double stability_metric(const std::vector<double>& torsions) {
        double sum = 0.0;
        double phi_power = 1.0;
        for (size_t n = 0; n < torsions.size(); ++n) {
            sum += std::abs(torsions[n]) * phi_power;
            phi_power *= PHI_INV;
        }
        return std::exp(-sum);
    }

    // Helix balance: Σκ(n)·φ⁻ⁿ
    static double helix_balance(const std::vector<double>& curvatures) {
        double sum = 0.0;
        double phi_power = 1.0;
        for (size_t n = 0; n < curvatures.size(); ++n) {
            sum += curvatures[n] * phi_power;
            phi_power *= PHI_INV;
        }
        return sum;
    }

    // G¹⁶ Mean curvature: (1/16) Σ κ_k × cos(2πk/10.6)
    static double mean_curvature_g16() {
        double sum = 0.0;
        double kappa = curvature();
        for (int k = 1; k <= NUM_STRANDS; ++k) {
            sum += kappa * std::cos(2.0 * M_PI * k / BASE_PAIRS_PER_TURN);
        }
        return sum / NUM_STRANDS;
    }

    // Full rotation calculation
    // θ_total = 16 × (2π/10.6) = 9.49 rad = 543.8°
    static double total_rotation() {
        return NUM_STRANDS * TWIST_ANGLE_RAD;
    }

    // Turns per 16 strands = 1.51
    static double turns_per_16_strands() {
        return total_rotation() / (2.0 * M_PI);
    }
};

// ============================================================================
// III. 44-BOND SYSTEM
// ============================================================================

enum class BondType {
    INTRA_STRAND,   // Bonds 1-16: within strand
    INTER_STRAND,   // Bonds 17-32: between adjacent strands
    CROSS_HELIX     // Bonds 33-44: across helix (|i-j| ∈ {4, 8, 12})
};

struct Bond {
    int id;              // Bond number 1-44
    BondType type;
    int strand_i;        // First strand (1-16)
    int strand_j;        // Second strand (1-16)
    double length;       // Bond length in nm
    double angle;        // Bond angle in radians
    double energy;       // Bond energy
    double occupation;   // Occupation state [0, 1]

    // Get bond position vector (midpoint)
    HelixPoint midpoint(double theta) const {
        auto p_i = G16Helix::point(strand_i, theta);
        auto p_j = G16Helix::point(strand_j, theta);
        return (p_i + p_j) * 0.5;
    }

    // Calculate bond length
    double calculate_length(double theta) const {
        auto p_i = G16Helix::point(strand_i, theta);
        auto p_j = G16Helix::point(strand_j, theta);
        return (p_i - p_j).norm();
    }
};

class BondSystem44 {
private:
    std::array<Bond, TOTAL_BONDS> bonds_;
    std::array<std::array<int, NUM_STRANDS>, NUM_STRANDS> adjacency_;

public:
    BondSystem44() {
        initialize_bonds();
        build_adjacency_matrix();
    }

    // Initialize all 44 bonds according to distribution (Section 2.1)
    void initialize_bonds() {
        int bond_id = 0;

        // Intra-strand bonds (1-16): B_intra(k) = 1 for k = 1..16
        for (int k = 1; k <= INTRA_STRAND_BONDS; ++k) {
            bonds_[bond_id++] = {
                bond_id, BondType::INTRA_STRAND,
                k, k,
                PITCH / NUM_STRANDS,  // L₀ = P/16 = 0.225 nm
                M_PI,  // 180° linear backbone
                0.0, 1.0
            };
        }

        // Inter-strand bonds (17-32): B_inter(k, k+1 mod 16) = 1
        for (int k = 1; k <= INTER_STRAND_BONDS; ++k) {
            int next_k = (k % NUM_STRANDS) + 1;
            bonds_[bond_id++] = {
                bond_id, BondType::INTER_STRAND,
                k, next_k,
                2.0 * HELIX_RADIUS * std::sin(M_PI / BASE_PAIRS_PER_TURN),  // d₀ = 0.59 nm
                STRAND_SEPARATION,  // 22.5°
                0.0, 1.0
            };
        }

        // Cross-helix bonds (33-44): pairs where |i-j| ∈ {4, 8, 12}
        // Pairs: (1,5), (2,6), (3,7), (4,8), (5,9), (6,10),
        //        (7,11), (8,12), (9,13), (10,14), (11,15), (12,16)
        int cross_pairs[][2] = {
            {1,5}, {2,6}, {3,7}, {4,8}, {5,9}, {6,10},
            {7,11}, {8,12}, {9,13}, {10,14}, {11,15}, {12,16}
        };

        for (int m = 0; m < CROSS_HELIX_BONDS; ++m) {
            int i = cross_pairs[m][0];
            int j = cross_pairs[m][1];
            int diff = std::abs(i - j);
            double angle = diff * STRAND_SEPARATION;  // |i-j| × 22.5°

            bonds_[bond_id++] = {
                bond_id, BondType::CROSS_HELIX,
                i, j,
                0.0,  // Calculated dynamically
                angle,
                0.0, 1.0
            };
        }
    }

    // Build adjacency matrix (Section 2.4)
    void build_adjacency_matrix() {
        // Initialize to zero
        for (auto& row : adjacency_) {
            row.fill(0);
        }

        // Fill based on bonds
        for (const auto& bond : bonds_) {
            int i = bond.strand_i - 1;
            int j = bond.strand_j - 1;
            adjacency_[i][j] = 1;
            adjacency_[j][i] = 1;
        }
    }

    // Verify bond count: Σ_{i,j} A_{ij}/2 = 44
    int verify_bond_count() const {
        int sum = 0;
        for (int i = 0; i < NUM_STRANDS; ++i) {
            for (int j = 0; j < NUM_STRANDS; ++j) {
                sum += adjacency_[i][j];
            }
        }
        // Add intra-strand bonds (diagonal)
        sum += INTRA_STRAND_BONDS;
        return sum / 2 + INTRA_STRAND_BONDS;
    }

    // Get bond by type and index
    Bond& get_bond(int id) { return bonds_[id - 1]; }
    const Bond& get_bond(int id) const { return bonds_[id - 1]; }

    // Get bonds by type
    std::vector<const Bond*> get_bonds_by_type(BondType type) const {
        std::vector<const Bond*> result;
        for (const auto& bond : bonds_) {
            if (bond.type == type) {
                result.push_back(&bond);
            }
        }
        return result;
    }

    // Bond coordinate system (Section 2.2)
    std::pair<int, int> bond_coordinates(int n) const {
        if (n <= INTRA_STRAND_BONDS) {
            // Intra-strand: (k, k)
            return {n, n};
        } else if (n <= INTRA_STRAND_BONDS + INTER_STRAND_BONDS) {
            // Inter-strand: (k, k+1 mod 16)
            int k = n - INTRA_STRAND_BONDS;
            return {k, (k % NUM_STRANDS) + 1};
        } else {
            // Cross-helix: (i, j) where |i-j| ∈ {4, 8, 12}
            return {bonds_[n-1].strand_i, bonds_[n-1].strand_j};
        }
    }

    // 10.6-adjusted bond arc length: s_bond = 2πr/10.6 = 0.593 nm per bp
    static double bond_arc_length() {
        return 2.0 * M_PI * HELIX_RADIUS / BASE_PAIRS_PER_TURN;
    }

    // Vertical rise per bond: Δz_bond = P/44 = 0.0819 nm
    static double vertical_rise_per_bond() {
        return PITCH / TOTAL_BONDS;
    }

    // Adjacent bond length: L_adj = √[2r²(1 - cos(2π/10.6)) + h²] = 0.54 nm
    static double adjacent_bond_length() {
        double cos_term = 1.0 - std::cos(TWIST_ANGLE_RAD);
        return std::sqrt(2.0 * HELIX_RADIUS * HELIX_RADIUS * cos_term +
                        RISE_PER_BP * RISE_PER_BP);
    }

    // Validation: distance constraint
    bool validate_distance(double epsilon_max) const {
        double adjustment = BASE_PAIRS_PER_TURN / 10.4;  // 10.6/10.4 correction
        for (const auto& bond : bonds_) {
            if (bond.length > epsilon_max * adjustment) {
                return false;
            }
        }
        return true;
    }

    // Validation: angular constraint
    bool validate_angular(double delta_theta) const {
        for (const auto& bond : bonds_) {
            double target = std::round(bond.angle / TWIST_ANGLE_RAD) * TWIST_ANGLE_RAD;
            if (std::abs(bond.angle - target) > delta_theta) {
                return false;
            }
        }
        return true;
    }

    // Update all bond states
    void update_bonds(double theta) {
        for (auto& bond : bonds_) {
            bond.length = bond.calculate_length(theta);
        }
    }

    // Get total active bonds
    double total_occupation() const {
        double sum = 0.0;
        for (const auto& bond : bonds_) {
            sum += bond.occupation;
        }
        return sum;
    }
};

// ============================================================================
// IV. STRAND MATHEMATICS
// ============================================================================

// Fibonacci sequence generator
inline uint64_t fibonacci(int n) {
    if (n <= 1) return n;
    uint64_t a = 0, b = 1;
    for (int i = 2; i <= n; ++i) {
        uint64_t c = a + b;
        a = b;
        b = c;
    }
    return b;
}

// Fibonacci-scaled strand: D_k = φ^k · F_k
inline double fibonacci_scaled_strand(int k) {
    return std::pow(PHI, k) * fibonacci(k);
}

class StrandMathematics {
private:
    std::array<double, NUM_STRANDS> D_;  // Fibonacci-scaled values
    std::array<std::complex<double>, NUM_STRANDS> S_;  // Complex representation
    std::array<std::complex<double>, NUM_STRANDS> T_;  // Transform

public:
    StrandMathematics() {
        initialize();
    }

    void initialize() {
        for (int k = 0; k < NUM_STRANDS; ++k) {
            D_[k] = fibonacci_scaled_strand(k + 1);
        }
        compute_transforms();
    }

    // 16-strand helix sum: Helix₁₆ = Σ D_k·sin(kφ)
    double helix_sum_16() const {
        double sum = 0.0;
        for (int k = 1; k <= NUM_STRANDS; ++k) {
            sum += D_[k-1] * std::sin(k * PHI);
        }
        return sum;
    }

    // Bond formation: Bond_G16 = Π(1 + D_k/Helix₁₆)
    double bond_formation_g16() const {
        double helix = helix_sum_16();
        double product = 1.0;
        for (int k = 0; k < NUM_STRANDS; ++k) {
            product *= (1.0 + D_[k] / helix);
        }
        return product;
    }

    // Crosslink function: Crosslink(i,j) = (D_i · D_j)/φ^{|i-j|}
    double crosslink(int i, int j) const {
        return (D_[i-1] * D_[j-1]) / std::pow(PHI, std::abs(i - j));
    }

    // Complex strand representation: S_k = D_k · e^{ikφ}
    void compute_transforms() {
        for (int k = 0; k < NUM_STRANDS; ++k) {
            double phase = (k + 1) * PHI;
            S_[k] = D_[k] * std::exp(std::complex<double>(0, phase));

            // Transform: T_k = S_k · (cos(k) + i·sin(k))
            double angle = k + 1;
            T_[k] = S_[k] * std::complex<double>(std::cos(angle), std::sin(angle));
        }
    }

    // Modular extraction: X_k = |T_k| mod φ^k
    double modular_extraction(int k) const {
        return std::fmod(std::abs(T_[k-1]), std::pow(PHI, k));
    }

    // Fusion integral: Fusion = Σ(Strand_Vector(k) · φ⁻ᵏ)
    double fusion_integral() const {
        double sum = 0.0;
        for (int k = 0; k < NUM_STRANDS; ++k) {
            sum += std::abs(T_[k]) * std::pow(PHI_INV, k + 1);
        }
        return sum;
    }

    // Get D values
    double get_D(int k) const { return D_[k-1]; }
    std::complex<double> get_S(int k) const { return S_[k-1]; }
    std::complex<double> get_T(int k) const { return T_[k-1]; }
};

// ============================================================================
// V. G¹⁶ TRANSFORMATION MATRICES
// ============================================================================

struct Matrix3x3 {
    std::array<std::array<double, 3>, 3> m;

    Matrix3x3 operator*(const Matrix3x3& other) const {
        Matrix3x3 result = {{{0}}};
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                for (int k = 0; k < 3; ++k) {
                    result.m[i][j] += m[i][k] * other.m[k][j];
                }
            }
        }
        return result;
    }

    HelixPoint transform(const HelixPoint& p) const {
        return {
            m[0][0]*p.x + m[0][1]*p.y + m[0][2]*p.z,
            m[1][0]*p.x + m[1][1]*p.y + m[1][2]*p.z,
            m[2][0]*p.x + m[2][1]*p.y + m[2][2]*p.z
        };
    }
};

class G16Transform {
public:
    // 22.5° increment rotation matrix (Section 5.1)
    static Matrix3x3 rotation_22_5(int k) {
        double angle = 2.0 * M_PI * k / NUM_STRANDS;
        double c = std::cos(angle);
        double s = std::sin(angle);
        return {{{
            {c, -s, 0},
            {s,  c, 0},
            {0,  0, 1}
        }}};
    }

    // 10.6-adjusted rotation matrix
    static Matrix3x3 rotation_10_6(int k) {
        double angle = 2.0 * M_PI * k / BASE_PAIRS_PER_TURN;
        double c = std::cos(angle);
        double s = std::sin(angle);
        return {{{
            {c, -s, 0},
            {s,  c, 0},
            {0,  0, 1}
        }}};
    }

    // Phase increment: Δφ = 2π/10.6 = 33.96° = 0.593 rad
    static double phase_increment() {
        return TWIST_ANGLE_RAD;
    }

    // Full G¹⁶ transform: G¹⁶ = Π R_k · T_k
    static Matrix3x3 full_transform() {
        Matrix3x3 result = {{{
            {1, 0, 0},
            {0, 1, 0},
            {0, 0, 1}
        }}};

        for (int k = 1; k <= NUM_STRANDS; ++k) {
            result = result * rotation_22_5(k);
        }
        return result;
    }

    // 10.6-adjusted full transform
    static Matrix3x3 full_transform_10_6() {
        Matrix3x3 result = {{{
            {1, 0, 0},
            {0, 1, 0},
            {0, 0, 1}
        }}};

        for (int k = 1; k <= NUM_STRANDS; ++k) {
            result = result * rotation_10_6(k);
        }
        return result;
    }
};

// ============================================================================
// VI. 10.6 FIRING PATTERN
// ============================================================================

class FiringPattern {
public:
    // Phase calculation: phase(i,t) = ⌊(t + offset_i)/10.6⌋ mod 16
    static int phase(int i, double t) {
        double offset = offset_distribution(i);
        return static_cast<int>(std::floor((t + offset) / BASE_PAIRS_PER_TURN)) % NUM_STRANDS;
    }

    // Offset distribution: offset_i = i · φ⁻¹ · 10.6
    static double offset_distribution(int i) {
        return i * PHI_INV * BASE_PAIRS_PER_TURN;
    }

    // Fire condition: fire(i,t) = [phase = 0]
    static bool should_fire(int i, double t) {
        return phase(i, t) == 0;
    }

    // Get all firing strands at time t
    static std::vector<int> firing_strands(double t) {
        std::vector<int> result;
        for (int i = 1; i <= NUM_STRANDS; ++i) {
            if (should_fire(i, t)) {
                result.push_back(i);
            }
        }
        return result;
    }

    // Synchronization period: T_sync = 10.6 × 16 = 169.6
    static double sync_period() {
        return SYNC_PERIOD;
    }

    // Fundamental frequency: f₀ = 1/10.6 = 0.0943
    static double fundamental_frequency() {
        return FUNDAMENTAL_FREQ;
    }

    // Harmonic: f_n = n/10.6
    static double harmonic(int n) {
        return n / BASE_PAIRS_PER_TURN;
    }

    // 44-bond resonance: f_44 = 44/10.6 = 4.15
    static double bond_resonance() {
        return BOND_RESONANCE_FREQ;
    }
};

// ============================================================================
// VII. GOLDEN BASIS FUNCTIONS
// ============================================================================

class GoldenBasis {
public:
    // Basis function: ψ_k(x) = sin(2πφ^k · x) + cos(2πφ^k · x)
    static double psi(int k, double x) {
        double arg = 2.0 * M_PI * std::pow(PHI, k) * x;
        return std::sin(arg) + std::cos(arg);
    }

    // Generate basis vector for N terms
    static std::vector<double> basis_vector(double x, int N) {
        std::vector<double> result(N);
        for (int k = 0; k < N; ++k) {
            result[k] = psi(k, x);
        }
        return result;
    }

    // Global coherence: Φ_total = (1/N) Σ Φ_i
    // where Φ_i = 1/(1 + |x_i - 1|)
    static double global_coherence(const std::vector<double>& x) {
        double sum = 0.0;
        for (double xi : x) {
            sum += 1.0 / (1.0 + std::abs(xi - 1.0));
        }
        return sum / x.size();
    }

    // Ring diffusion: dS_i/dt = λ(S_{i+1} + S_{i-1} - 2S_i)
    static void ring_diffusion_step(std::vector<double>& S, double lambda, double dt) {
        int N = S.size();
        std::vector<double> dS(N);

        for (int i = 0; i < N; ++i) {
            int prev = (i - 1 + N) % N;
            int next = (i + 1) % N;
            dS[i] = lambda * (S[next] + S[prev] - 2.0 * S[i]);
        }

        for (int i = 0; i < N; ++i) {
            S[i] += dS[i] * dt;
        }
    }
};

// ============================================================================
// VIII. SUPERCOILING (10.6)
// ============================================================================

class Supercoiling {
public:
    // Linking number: Lk = Tw + Wr
    static double linking_number(double twist, double writhe) {
        return twist + writhe;
    }

    // Twist for 16 strands: Tw = N/10.6 × 16
    static double twist_16_strands(int N) {
        return (static_cast<double>(N) / BASE_PAIRS_PER_TURN) * NUM_STRANDS;
    }

    // Relaxed linking number: Lk₀ = N/10.6
    static double relaxed_linking_number(int N) {
        return static_cast<double>(N) / BASE_PAIRS_PER_TURN;
    }

    // Superhelical density: σ = (Lk - Lk₀)/Lk₀
    static double superhelical_density(double Lk, double Lk0) {
        return (Lk - Lk0) / Lk0;
    }

    // G¹⁶ superhelical density
    static double superhelical_density_g16(double Lk, int N) {
        double Lk0_g16 = NUM_STRANDS * static_cast<double>(N) / BASE_PAIRS_PER_TURN;
        return (Lk - Lk0_g16) / Lk0_g16;
    }

    // Free energy of supercoiling
    // ΔG_sc = K · N · σ² · (10.6/10.4)²
    static double free_energy(double K, int N, double sigma) {
        constexpr double correction = (BASE_PAIRS_PER_TURN / 10.4) * (BASE_PAIRS_PER_TURN / 10.4);
        return K * N * sigma * sigma * correction;
    }

    // G¹⁶ with 44-bond correction
    static double free_energy_g16(double K, int N, double sigma_g16, double lambda) {
        constexpr double correction = (BASE_PAIRS_PER_TURN / 10.4) * (BASE_PAIRS_PER_TURN / 10.4);
        return K * N * sigma_g16 * sigma_g16 * (1.0 + lambda * TOTAL_BONDS) * correction;
    }
};

// ============================================================================
// COMPLETE G¹⁶ HELIX SYSTEM
// ============================================================================

class G16HelixSystem {
private:
    G16Helix helix_;
    BondSystem44 bonds_;
    StrandMathematics strands_;
    double time_;
    double global_phase_;  // φ_10.6(t)

    // State vectors
    std::array<std::complex<double>, NUM_STRANDS> strand_states_;
    std::array<double, NUM_STRANDS> angular_positions_;

public:
    G16HelixSystem() : time_(0.0), global_phase_(0.0) {
        initialize();
    }

    void initialize() {
        // Initialize strand states: S_k(t=0) = A_k · e^{i(2πk/10.6)}
        for (int k = 0; k < NUM_STRANDS; ++k) {
            double amplitude = strands_.get_D(k + 1);
            double phase = 2.0 * M_PI * (k + 1) / BASE_PAIRS_PER_TURN;
            strand_states_[k] = amplitude * std::exp(std::complex<double>(0, phase));
            angular_positions_[k] = phase;
        }
    }

    // Tick the system forward
    void tick(double dt) {
        time_ += dt;

        // Update global phase: φ(t) = (2πt/10.6) mod 2π
        global_phase_ = std::fmod(2.0 * M_PI * time_ / BASE_PAIRS_PER_TURN, 2.0 * M_PI);

        // Update strand states: S_k(t) = A_k · e^{i(ωt + 2πk/10.6)}
        for (int k = 0; k < NUM_STRANDS; ++k) {
            double amplitude = strands_.get_D(k + 1);
            double phase = ANGULAR_FREQUENCY * time_ + 2.0 * M_PI * (k + 1) / BASE_PAIRS_PER_TURN;
            strand_states_[k] = amplitude * std::exp(std::complex<double>(0, phase));
            angular_positions_[k] = phase;
        }

        // Update bond states
        bonds_.update_bonds(global_phase_);
    }

    // Get firing strands
    std::vector<int> get_firing_strands() const {
        return FiringPattern::firing_strands(time_);
    }

    // Get global coherence
    double get_coherence() const {
        std::vector<double> magnitudes(NUM_STRANDS);
        for (int k = 0; k < NUM_STRANDS; ++k) {
            magnitudes[k] = std::abs(strand_states_[k]);
        }
        return GoldenBasis::global_coherence(magnitudes);
    }

    // Get total bond occupation
    double get_bond_occupation() const {
        return bonds_.total_occupation();
    }

    // Accessors
    double time() const { return time_; }
    double global_phase() const { return global_phase_; }
    const BondSystem44& bonds() const { return bonds_; }
    BondSystem44& bonds() { return bonds_; }
    const StrandMathematics& strands() const { return strands_; }
};

} // namespace agi::dna