// Package dna implements the G¹⁶ Helix System with 44 bonds and 10.6 periodicity
// Complete DNA Mathematical Architecture for AGI substrate
package dna

import (
	"math"
	"math/cmplx"
)

// =============================================================================
// FOUNDATIONAL CONSTANTS
// =============================================================================

const (
	// Golden ratio and derived constants
	PHI    = 1.6180339887498948482 // (1 + √5) / 2
	PHIInv = 0.6180339887498948482 // 1 / φ = φ - 1

	// DNA Helix parameters
	BasePairsPerTurn = 10.6       // p = 10.6 bp
	HelixRadius      = 1.0        // r = 1.0 nm
	RisePerBP        = 0.34       // h = 0.34 nm
	Pitch            = 3.604      // P = h × p = 0.34 × 10.6 nm
	TwistAngleDeg    = 33.96      // 360°/10.6
	TwistAngleRad    = 0.593      // 2π/10.6 rad
	AngularFrequency = 0.593      // ω = 2π/10.6 rad/bp

	// Strand parameters
	NumStrands       = 16
	StrandSeparation = 0.3927 // 2π/16 = 22.5°
	NodesPerStrand   = 312
	TotalNodes       = 5000
	HubNodes         = 8

	// Bond parameters
	IntraStrandBonds = 16
	InterStrandBonds = 16
	CrossHelixBonds  = 12
	TotalBonds       = 44

	// Curvature and torsion
	HelixCurvature = 0.745 // κ nm⁻¹
	HelixTorsion   = 0.427 // τ nm⁻¹

	// Synchronization
	SyncPeriod        = 169.6 // 10.6 × 16
	FundamentalFreq   = 0.0943
	BondResonanceFreq = 4.15 // 44/10.6

	// Steady state
	SteadyStateBonds = 36.2
	BondOccupation   = 0.82
)

// =============================================================================
// G¹⁶ PARAMETRIC HELIX
// =============================================================================

// HelixPoint represents a 3D point on the helix
type HelixPoint struct {
	X, Y, Z float64
}

// Add returns the sum of two points
func (p HelixPoint) Add(other HelixPoint) HelixPoint {
	return HelixPoint{p.X + other.X, p.Y + other.Y, p.Z + other.Z}
}

// Sub returns the difference of two points
func (p HelixPoint) Sub(other HelixPoint) HelixPoint {
	return HelixPoint{p.X - other.X, p.Y - other.Y, p.Z - other.Z}
}

// Scale returns the point scaled by s
func (p HelixPoint) Scale(s float64) HelixPoint {
	return HelixPoint{p.X * s, p.Y * s, p.Z * s}
}

// Norm returns the Euclidean norm
func (p HelixPoint) Norm() float64 {
	return math.Sqrt(p.X*p.X + p.Y*p.Y + p.Z*p.Z)
}

// Dot returns the dot product
func (p HelixPoint) Dot(other HelixPoint) float64 {
	return p.X*other.X + p.Y*other.Y + p.Z*other.Z
}

// G16Helix implements G¹⁶ parametric helix equations
type G16Helix struct{}

// VerticalOffset returns δ_k = k·h/16
func (G16Helix) VerticalOffset(k int) float64 {
	return float64(k) * RisePerBP / NumStrands
}

// GoldenPhaseOffset returns θ_k = 2πk·φ⁻¹ mod 2π
func (G16Helix) GoldenPhaseOffset(k int) float64 {
	return math.Mod(2.0*math.Pi*float64(k)*PHIInv, 2.0*math.Pi)
}

// Point returns point on strand k at parameter θ
// x_k(θ) = r·cos(θ + 2πk/16)
// y_k(θ) = r·sin(θ + 2πk/16)
// z_k(θ) = (p/2π)·θ + δ_k
func (h G16Helix) Point(k int, theta float64) HelixPoint {
	phase := theta + 2.0*math.Pi*float64(k)/NumStrands
	return HelixPoint{
		X: HelixRadius * math.Cos(phase),
		Y: HelixRadius * math.Sin(phase),
		Z: (Pitch/(2.0*math.Pi))*theta + h.VerticalOffset(k),
	}
}

// PointGolden returns point with golden ratio phase adjustment
func (h G16Helix) PointGolden(k int, theta float64) HelixPoint {
	phase := theta + h.GoldenPhaseOffset(k)
	return HelixPoint{
		X: HelixRadius * math.Cos(phase),
		Y: HelixRadius * math.Sin(phase),
		Z: (Pitch/(2.0*math.Pi))*theta + h.VerticalOffset(k),
	}
}

// Curvature returns κ = r/(r² + (p·h/2π)²)
func (G16Helix) Curvature() float64 {
	phTerm := Pitch * RisePerBP / (2.0 * math.Pi)
	return HelixRadius / (HelixRadius*HelixRadius + phTerm*phTerm)
}

// Torsion returns τ = (p·h/2π)/(r² + (p·h/2π)²)
func (G16Helix) Torsion() float64 {
	phTerm := Pitch * RisePerBP / (2.0 * math.Pi)
	return phTerm / (HelixRadius*HelixRadius + phTerm*phTerm)
}

// StabilityMetric returns exp(-Σ|τ(n)|·φ⁻ⁿ)
func (G16Helix) StabilityMetric(torsions []float64) float64 {
	sum := 0.0
	phiPower := 1.0
	for _, t := range torsions {
		sum += math.Abs(t) * phiPower
		phiPower *= PHIInv
	}
	return math.Exp(-sum)
}

// HelixBalance returns Σκ(n)·φ⁻ⁿ
func (G16Helix) HelixBalance(curvatures []float64) float64 {
	sum := 0.0
	phiPower := 1.0
	for _, k := range curvatures {
		sum += k * phiPower
		phiPower *= PHIInv
	}
	return sum
}

// MeanCurvatureG16 returns (1/16) Σ κ_k × cos(2πk/10.6)
func (h G16Helix) MeanCurvatureG16() float64 {
	sum := 0.0
	kappa := h.Curvature()
	for k := 1; k <= NumStrands; k++ {
		sum += kappa * math.Cos(2.0*math.Pi*float64(k)/BasePairsPerTurn)
	}
	return sum / NumStrands
}

// TotalRotation returns θ_total = 16 × (2π/10.6) = 9.49 rad
func (G16Helix) TotalRotation() float64 {
	return NumStrands * TwistAngleRad
}

// TurnsPer16Strands returns 1.51
func (h G16Helix) TurnsPer16Strands() float64 {
	return h.TotalRotation() / (2.0 * math.Pi)
}

// =============================================================================
// 44-BOND SYSTEM
// =============================================================================

// BondType represents the type of bond
type BondType int

const (
	BondIntraStrand BondType = iota // Bonds 1-16
	BondInterStrand                 // Bonds 17-32
	BondCrossHelix                  // Bonds 33-44
)

// Bond represents a single bond in the 44-bond system
type Bond struct {
	ID         int
	Type       BondType
	StrandI    int     // First strand (1-16)
	StrandJ    int     // Second strand (1-16)
	Length     float64 // Bond length in nm
	Angle      float64 // Bond angle in radians
	Energy     float64 // Bond energy
	Occupation float64 // Occupation state [0, 1]
}

// Midpoint returns the bond position vector
func (b *Bond) Midpoint(theta float64) HelixPoint {
	helix := G16Helix{}
	pi := helix.Point(b.StrandI, theta)
	pj := helix.Point(b.StrandJ, theta)
	return pi.Add(pj).Scale(0.5)
}

// CalculateLength calculates the bond length at given theta
func (b *Bond) CalculateLength(theta float64) float64 {
	helix := G16Helix{}
	pi := helix.Point(b.StrandI, theta)
	pj := helix.Point(b.StrandJ, theta)
	return pi.Sub(pj).Norm()
}

// BondSystem44 manages all 44 bonds
type BondSystem44 struct {
	Bonds     [TotalBonds]Bond
	Adjacency [NumStrands][NumStrands]int
}

// NewBondSystem44 creates and initializes a new 44-bond system
func NewBondSystem44() *BondSystem44 {
	bs := &BondSystem44{}
	bs.initializeBonds()
	bs.buildAdjacencyMatrix()
	return bs
}

// initializeBonds sets up all 44 bonds
func (bs *BondSystem44) initializeBonds() {
	bondID := 0

	// Intra-strand bonds (1-16)
	for k := 1; k <= IntraStrandBonds; k++ {
		bs.Bonds[bondID] = Bond{
			ID:         bondID + 1,
			Type:       BondIntraStrand,
			StrandI:    k,
			StrandJ:    k,
			Length:     Pitch / NumStrands, // L₀ = 0.225 nm
			Angle:      math.Pi,            // 180°
			Occupation: 1.0,
		}
		bondID++
	}

	// Inter-strand bonds (17-32)
	for k := 1; k <= InterStrandBonds; k++ {
		nextK := (k % NumStrands) + 1
		bs.Bonds[bondID] = Bond{
			ID:         bondID + 1,
			Type:       BondInterStrand,
			StrandI:    k,
			StrandJ:    nextK,
			Length:     2.0 * HelixRadius * math.Sin(math.Pi/BasePairsPerTurn), // d₀ ≈ 0.59 nm
			Angle:      StrandSeparation,                                       // 22.5°
			Occupation: 1.0,
		}
		bondID++
	}

	// Cross-helix bonds (33-44)
	crossPairs := [][2]int{
		{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 _, pair := range crossPairs {
		i, j := pair[0], pair[1]
		diff := j - i
		angle := float64(diff) * StrandSeparation

		bs.Bonds[bondID] = Bond{
			ID:         bondID + 1,
			Type:       BondCrossHelix,
			StrandI:    i,
			StrandJ:    j,
			Length:     0.0, // Calculated dynamically
			Angle:      angle,
			Occupation: 1.0,
		}
		bondID++
	}
}

// buildAdjacencyMatrix builds the adjacency matrix
func (bs *BondSystem44) buildAdjacencyMatrix() {
	// Initialize to zero
	for i := range bs.Adjacency {
		for j := range bs.Adjacency[i] {
			bs.Adjacency[i][j] = 0
		}
	}

	// Fill based on bonds
	for _, bond := range bs.Bonds {
		if bond.StrandI != bond.StrandJ {
			i := bond.StrandI - 1
			j := bond.StrandJ - 1
			bs.Adjacency[i][j] = 1
			bs.Adjacency[j][i] = 1
		}
	}
}

// VerifyBondCount verifies Σ A_{ij}/2 = 44
func (bs *BondSystem44) VerifyBondCount() int {
	sum := 0
	for i := 0; i < NumStrands; i++ {
		for j := 0; j < NumStrands; j++ {
			sum += bs.Adjacency[i][j]
		}
	}
	return sum/2 + IntraStrandBonds
}

// GetBond returns the bond by ID (1-indexed)
func (bs *BondSystem44) GetBond(id int) *Bond {
	if id < 1 || id > TotalBonds {
		return nil
	}
	return &bs.Bonds[id-1]
}

// GetBondsByType returns all bonds of a given type
func (bs *BondSystem44) GetBondsByType(t BondType) []*Bond {
	var result []*Bond
	for i := range bs.Bonds {
		if bs.Bonds[i].Type == t {
			result = append(result, &bs.Bonds[i])
		}
	}
	return result
}

// BondArcLength returns s_bond = 2πr/10.6 = 0.593 nm per bp
func BondArcLength() float64 {
	return 2.0 * math.Pi * HelixRadius / BasePairsPerTurn
}

// VerticalRisePerBond returns Δz_bond = P/44 = 0.0819 nm
func VerticalRisePerBond() float64 {
	return Pitch / TotalBonds
}

// AdjacentBondLength returns L_adj = √[2r²(1 - cos(2π/10.6)) + h²] = 0.54 nm
func AdjacentBondLength() float64 {
	cosTerm := 1.0 - math.Cos(TwistAngleRad)
	return math.Sqrt(2.0*HelixRadius*HelixRadius*cosTerm + RisePerBP*RisePerBP)
}

// UpdateBonds updates all bond states at given theta
func (bs *BondSystem44) UpdateBonds(theta float64) {
	for i := range bs.Bonds {
		bs.Bonds[i].Length = bs.Bonds[i].CalculateLength(theta)
	}
}

// TotalOccupation returns the sum of all bond occupations
func (bs *BondSystem44) TotalOccupation() float64 {
	sum := 0.0
	for _, bond := range bs.Bonds {
		sum += bond.Occupation
	}
	return sum
}

// =============================================================================
// STRAND MATHEMATICS
// =============================================================================

// Fibonacci returns the n-th Fibonacci number
func Fibonacci(n int) uint64 {
	if n <= 1 {
		return uint64(n)
	}
	a, b := uint64(0), uint64(1)
	for i := 2; i <= n; i++ {
		a, b = b, a+b
	}
	return b
}

// FibonacciScaledStrand returns D_k = φ^k · F_k
func FibonacciScaledStrand(k int) float64 {
	return math.Pow(PHI, float64(k)) * float64(Fibonacci(k))
}

// StrandMathematics handles strand calculations
type StrandMathematics struct {
	D [NumStrands]float64         // Fibonacci-scaled values
	S [NumStrands]complex128      // Complex representation
	T [NumStrands]complex128      // Transform
}

// NewStrandMathematics creates and initializes strand mathematics
func NewStrandMathematics() *StrandMathematics {
	sm := &StrandMathematics{}
	sm.initialize()
	return sm
}

func (sm *StrandMathematics) initialize() {
	for k := 0; k < NumStrands; k++ {
		sm.D[k] = FibonacciScaledStrand(k + 1)
	}
	sm.computeTransforms()
}

// HelixSum16 returns Helix₁₆ = Σ D_k·sin(kφ)
func (sm *StrandMathematics) HelixSum16() float64 {
	sum := 0.0
	for k := 1; k <= NumStrands; k++ {
		sum += sm.D[k-1] * math.Sin(float64(k)*PHI)
	}
	return sum
}

// BondFormationG16 returns Bond_G16 = Π(1 + D_k/Helix₁₆)
func (sm *StrandMathematics) BondFormationG16() float64 {
	helix := sm.HelixSum16()
	product := 1.0
	for k := 0; k < NumStrands; k++ {
		product *= (1.0 + sm.D[k]/helix)
	}
	return product
}

// Crosslink returns Crosslink(i,j) = (D_i · D_j)/φ^{|i-j|}
func (sm *StrandMathematics) Crosslink(i, j int) float64 {
	diff := i - j
	if diff < 0 {
		diff = -diff
	}
	return (sm.D[i-1] * sm.D[j-1]) / math.Pow(PHI, float64(diff))
}

// computeTransforms computes complex strand representations
func (sm *StrandMathematics) computeTransforms() {
	for k := 0; k < NumStrands; k++ {
		phase := float64(k+1) * PHI
		sm.S[k] = complex(sm.D[k], 0) * cmplx.Exp(complex(0, phase))

		// Transform: T_k = S_k · (cos(k) + i·sin(k))
		angle := float64(k + 1)
		sm.T[k] = sm.S[k] * complex(math.Cos(angle), math.Sin(angle))
	}
}

// ModularExtraction returns X_k = |T_k| mod φ^k
func (sm *StrandMathematics) ModularExtraction(k int) float64 {
	return math.Mod(cmplx.Abs(sm.T[k-1]), math.Pow(PHI, float64(k)))
}

// FusionIntegral returns Fusion = Σ(|T_k| · φ⁻ᵏ)
func (sm *StrandMathematics) FusionIntegral() float64 {
	sum := 0.0
	for k := 0; k < NumStrands; k++ {
		sum += cmplx.Abs(sm.T[k]) * math.Pow(PHIInv, float64(k+1))
	}
	return sum
}

// =============================================================================
// G¹⁶ TRANSFORMATION
// =============================================================================

// Matrix3x3 represents a 3×3 matrix
type Matrix3x3 [3][3]float64

// Identity returns the identity matrix
func Identity() Matrix3x3 {
	return Matrix3x3{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}
}

// Multiply multiplies two matrices
func (m Matrix3x3) Multiply(other Matrix3x3) Matrix3x3 {
	var result Matrix3x3
	for i := 0; i < 3; i++ {
		for j := 0; j < 3; j++ {
			for k := 0; k < 3; k++ {
				result[i][j] += m[i][k] * other[k][j]
			}
		}
	}
	return result
}

// Transform applies the matrix to a point
func (m Matrix3x3) Transform(p HelixPoint) HelixPoint {
	return HelixPoint{
		X: m[0][0]*p.X + m[0][1]*p.Y + m[0][2]*p.Z,
		Y: m[1][0]*p.X + m[1][1]*p.Y + m[1][2]*p.Z,
		Z: m[2][0]*p.X + m[2][1]*p.Y + m[2][2]*p.Z,
	}
}

// Rotation22_5 returns 22.5° increment rotation matrix
func Rotation22_5(k int) Matrix3x3 {
	angle := 2.0 * math.Pi * float64(k) / NumStrands
	c := math.Cos(angle)
	s := math.Sin(angle)
	return Matrix3x3{
		{c, -s, 0},
		{s, c, 0},
		{0, 0, 1},
	}
}

// Rotation10_6 returns 10.6-adjusted rotation matrix
func Rotation10_6(k int) Matrix3x3 {
	angle := 2.0 * math.Pi * float64(k) / BasePairsPerTurn
	c := math.Cos(angle)
	s := math.Sin(angle)
	return Matrix3x3{
		{c, -s, 0},
		{s, c, 0},
		{0, 0, 1},
	}
}

// FullG16Transform returns G¹⁶ = Π R_k
func FullG16Transform() Matrix3x3 {
	result := Identity()
	for k := 1; k <= NumStrands; k++ {
		result = result.Multiply(Rotation22_5(k))
	}
	return result
}

// FullG16Transform10_6 returns 10.6-adjusted full transform
func FullG16Transform10_6() Matrix3x3 {
	result := Identity()
	for k := 1; k <= NumStrands; k++ {
		result = result.Multiply(Rotation10_6(k))
	}
	return result
}

// =============================================================================
// 10.6 FIRING PATTERN
// =============================================================================

// FiringPattern handles 10.6 phase-based firing
type FiringPattern struct{}

// Phase returns phase(i,t) = ⌊(t + offset_i)/10.6⌋ mod 16
func (FiringPattern) Phase(i int, t float64) int {
	offset := FiringPattern{}.OffsetDistribution(i)
	return int(math.Floor((t+offset)/BasePairsPerTurn)) % NumStrands
}

// OffsetDistribution returns offset_i = i · φ⁻¹ · 10.6
func (FiringPattern) OffsetDistribution(i int) float64 {
	return float64(i) * PHIInv * BasePairsPerTurn
}

// ShouldFire returns true if strand i should fire at time t
func (fp FiringPattern) ShouldFire(i int, t float64) bool {
	return fp.Phase(i, t) == 0
}

// FiringStrands returns all strands that should fire at time t
func (fp FiringPattern) FiringStrands(t float64) []int {
	var result []int
	for i := 1; i <= NumStrands; i++ {
		if fp.ShouldFire(i, t) {
			result = append(result, i)
		}
	}
	return result
}

// Harmonic returns f_n = n/10.6
func (FiringPattern) Harmonic(n int) float64 {
	return float64(n) / BasePairsPerTurn
}

// =============================================================================
// GOLDEN BASIS FUNCTIONS
// =============================================================================

// GoldenBasis provides golden ratio basis functions
type GoldenBasis struct{}

// Psi returns ψ_k(x) = sin(2πφ^k · x) + cos(2πφ^k · x)
func (GoldenBasis) Psi(k int, x float64) float64 {
	arg := 2.0 * math.Pi * math.Pow(PHI, float64(k)) * x
	return math.Sin(arg) + math.Cos(arg)
}

// BasisVector generates basis vector for N terms
func (gb GoldenBasis) BasisVector(x float64, N int) []float64 {
	result := make([]float64, N)
	for k := 0; k < N; k++ {
		result[k] = gb.Psi(k, x)
	}
	return result
}

// GlobalCoherence returns Φ_total = (1/N) Σ Φ_i where Φ_i = 1/(1 + |x_i - 1|)
func (GoldenBasis) GlobalCoherence(x []float64) float64 {
	sum := 0.0
	for _, xi := range x {
		sum += 1.0 / (1.0 + math.Abs(xi-1.0))
	}
	return sum / float64(len(x))
}

// RingDiffusionStep performs dS_i/dt = λ(S_{i+1} + S_{i-1} - 2S_i)
func (GoldenBasis) RingDiffusionStep(S []float64, lambda, dt float64) {
	N := len(S)
	dS := make([]float64, N)

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

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

// =============================================================================
// SUPERCOILING
// =============================================================================

// Supercoiling handles DNA supercoiling calculations
type Supercoiling struct{}

// LinkingNumber returns Lk = Tw + Wr
func (Supercoiling) LinkingNumber(twist, writhe float64) float64 {
	return twist + writhe
}

// Twist16Strands returns Tw = N/10.6 × 16
func (Supercoiling) Twist16Strands(N int) float64 {
	return (float64(N) / BasePairsPerTurn) * NumStrands
}

// RelaxedLinkingNumber returns Lk₀ = N/10.6
func (Supercoiling) RelaxedLinkingNumber(N int) float64 {
	return float64(N) / BasePairsPerTurn
}

// SuperhelicalDensity returns σ = (Lk - Lk₀)/Lk₀
func (Supercoiling) SuperhelicalDensity(Lk, Lk0 float64) float64 {
	return (Lk - Lk0) / Lk0
}

// SuperhelicalDensityG16 returns G¹⁶ superhelical density
func (Supercoiling) SuperhelicalDensityG16(Lk float64, N int) float64 {
	Lk0G16 := NumStrands * float64(N) / BasePairsPerTurn
	return (Lk - Lk0G16) / Lk0G16
}

// FreeEnergy returns ΔG_sc = K · N · σ² · (10.6/10.4)²
func (Supercoiling) FreeEnergy(K float64, N int, sigma float64) float64 {
	correction := (BasePairsPerTurn / 10.4) * (BasePairsPerTurn / 10.4)
	return K * float64(N) * sigma * sigma * correction
}

// FreeEnergyG16 returns G¹⁶ with 44-bond correction
func (Supercoiling) FreeEnergyG16(K float64, N int, sigmaG16, lambda float64) float64 {
	correction := (BasePairsPerTurn / 10.4) * (BasePairsPerTurn / 10.4)
	return K * float64(N) * sigmaG16 * sigmaG16 * (1.0 + lambda*TotalBonds) * correction
}

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

// G16HelixSystem is the complete helix system
type G16HelixSystem struct {
	Helix         G16Helix
	Bonds         *BondSystem44
	Strands       *StrandMathematics
	Time          float64
	GlobalPhase   float64
	StrandStates  [NumStrands]complex128
	AngularPos    [NumStrands]float64
}

// NewG16HelixSystem creates a new helix system
func NewG16HelixSystem() *G16HelixSystem {
	sys := &G16HelixSystem{
		Bonds:   NewBondSystem44(),
		Strands: NewStrandMathematics(),
	}
	sys.initialize()
	return sys
}

func (sys *G16HelixSystem) initialize() {
	for k := 0; k < NumStrands; k++ {
		amplitude := sys.Strands.D[k]
		phase := 2.0 * math.Pi * float64(k+1) / BasePairsPerTurn
		sys.StrandStates[k] = complex(amplitude, 0) * cmplx.Exp(complex(0, phase))
		sys.AngularPos[k] = phase
	}
}

// Tick advances the system by dt
func (sys *G16HelixSystem) Tick(dt float64) {
	sys.Time += dt

	// Update global phase: φ(t) = (2πt/10.6) mod 2π
	sys.GlobalPhase = math.Mod(2.0*math.Pi*sys.Time/BasePairsPerTurn, 2.0*math.Pi)

	// Update strand states: S_k(t) = A_k · e^{i(ωt + 2πk/10.6)}
	for k := 0; k < NumStrands; k++ {
		amplitude := sys.Strands.D[k]
		phase := AngularFrequency*sys.Time + 2.0*math.Pi*float64(k+1)/BasePairsPerTurn
		sys.StrandStates[k] = complex(amplitude, 0) * cmplx.Exp(complex(0, phase))
		sys.AngularPos[k] = phase
	}

	// Update bond states
	sys.Bonds.UpdateBonds(sys.GlobalPhase)
}

// GetFiringStrands returns currently firing strands
func (sys *G16HelixSystem) GetFiringStrands() []int {
	return FiringPattern{}.FiringStrands(sys.Time)
}

// GetCoherence returns global coherence
func (sys *G16HelixSystem) GetCoherence() float64 {
	magnitudes := make([]float64, NumStrands)
	for k := 0; k < NumStrands; k++ {
		magnitudes[k] = cmplx.Abs(sys.StrandStates[k])
	}
	return GoldenBasis{}.GlobalCoherence(magnitudes)
}

// GetBondOccupation returns total bond occupation
func (sys *G16HelixSystem) GetBondOccupation() float64 {
	return sys.Bonds.TotalOccupation()
}