// Package dna - Bond Dynamics, Rate Equations, and Thermodynamics
// Sections IX, X, XII of DNA Mathematical Architecture
package dna

import (
	"math"
)

// =============================================================================
// IX. 44-BOND ENERGY FUNCTIONS
// =============================================================================

// BondEnergy calculates energy functions for the 44-bond system
type BondEnergy struct {
	// Spring constants
	KappaIntra float64 // Intra-strand stiffness
	KappaInter float64 // Inter-strand stiffness
	Epsilon    float64 // Cosine modulation amplitude

	// Lennard-Jones parameters
	EpsilonCross float64
	SigmaCross   float64

	// 10.6 periodicity coupling
	KappaPeriod float64
}

// NewBondEnergy creates a new BondEnergy with default parameters
func NewBondEnergy() *BondEnergy {
	return &BondEnergy{
		KappaIntra:   100.0,
		KappaInter:   50.0,
		Epsilon:      1.0,
		EpsilonCross: 2.0,
		SigmaCross:   0.5,
		KappaPeriod:  10.0,
	}
}

// Equilibrium lengths
const (
	L0Intra = Pitch / NumStrands                                   // 0.225 nm
	D0Inter = 2.0 * HelixRadius * 0.296198 // sin(π/10.6) ≈ 0.59 nm
)

// IntraStrandEnergy: E_intra(k) = ½κ_intra(L_k - L₀)²
func (be *BondEnergy) IntraStrandEnergy(Lk float64) float64 {
	dL := Lk - L0Intra
	return 0.5 * be.KappaIntra * dL * dL
}

// InterStrandEnergy: E_inter(k, k+1) = ½κ_inter(d_{k,k+1} - d₀)² + ε·cos(2π/10.6 · k)
func (be *BondEnergy) InterStrandEnergy(dk float64, k int) float64 {
	dd := dk - D0Inter
	harmonic := 0.5 * be.KappaInter * dd * dd
	periodic := be.Epsilon * math.Cos(TwistAngleRad*float64(k))
	return harmonic + periodic
}

// CrossHelixEnergy: E_cross(i,j) = -ε_cross/d_{ij}⁶ + σ_cross/d_{ij}¹²
func (be *BondEnergy) CrossHelixEnergy(dij float64) float64 {
	if dij < 0.01 {
		return 0.0
	}
	d6 := math.Pow(dij, 6)
	d12 := d6 * d6
	return -be.EpsilonCross/d6 + be.SigmaCross/d12
}

// PeriodicityEnergy: E_10.6 = -κ_period · cos(2πN/10.6)
func (be *BondEnergy) PeriodicityEnergy(N int) float64 {
	return -be.KappaPeriod * math.Cos(2.0*math.Pi*float64(N)/BasePairsPerTurn)
}

// TotalEnergy calculates total system energy
func (be *BondEnergy) TotalEnergy(bonds *BondSystem44, N int) float64 {
	ETotal := 0.0

	// Intra-strand (bonds 1-16)
	for n := 1; n <= IntraStrandBonds; n++ {
		bond := bonds.GetBond(n)
		ETotal += be.IntraStrandEnergy(bond.Length)
	}

	// Inter-strand (bonds 17-32)
	for n := IntraStrandBonds + 1; n <= IntraStrandBonds+InterStrandBonds; n++ {
		bond := bonds.GetBond(n)
		k := n - IntraStrandBonds
		ETotal += be.InterStrandEnergy(bond.Length, k)
	}

	// Cross-helix (bonds 33-44)
	for n := IntraStrandBonds + InterStrandBonds + 1; n <= TotalBonds; n++ {
		bond := bonds.GetBond(n)
		if bond.Length > 0.01 {
			ETotal += be.CrossHelixEnergy(bond.Length)
		}
	}

	// Add 10.6 periodicity
	ETotal += be.PeriodicityEnergy(N)

	return ETotal
}

// EnergyGradient calculates energy gradient for bond n
func (be *BondEnergy) EnergyGradient(bond *Bond, n int) float64 {
	if n <= IntraStrandBonds {
		return be.KappaIntra * (bond.Length - L0Intra)
	} else if n <= IntraStrandBonds+InterStrandBonds {
		return be.KappaInter * (bond.Length - D0Inter)
	} else {
		// LJ gradient: dE/dr = 6ε/r⁷ - 12σ/r¹³
		r := bond.Length
		if r < 0.01 {
			return 0.0
		}
		r7 := math.Pow(r, 7)
		r13 := math.Pow(r, 13)
		return 6.0*be.EpsilonCross/r7 - 12.0*be.SigmaCross/r13
	}
}

// =============================================================================
// X. THERMODYNAMICS
// =============================================================================

// Thermodynamics handles thermodynamic calculations
type Thermodynamics struct {
	Temperature float64 // Kelvin
}

const (
	kB = 1.380649e-23 // Boltzmann constant J/K
	R  = 8.314        // Gas constant J/(mol·K)
)

// NewThermodynamics creates a new Thermodynamics with default temperature
func NewThermodynamics(T float64) *Thermodynamics {
	return &Thermodynamics{Temperature: T}
}

// Theta10_6 is the 10.6-periodic weighting function
func Theta10_6(configParam float64) float64 {
	return math.Exp(-math.Pow(math.Mod(configParam, BasePairsPerTurn), 2) / 2.0)
}

// PartitionFunction: Z = Σ exp(-E/k_B T) · Θ(10.6)
func (t *Thermodynamics) PartitionFunction(energies, configParams []float64) float64 {
	Z := 0.0
	for i, E := range energies {
		boltzmann := math.Exp(-E / (kB * t.Temperature))
		weight := Theta10_6(configParams[i])
		Z += boltzmann * weight
	}
	return Z
}

// PerBondFreeEnergy: ΔG_bond = -k_B T · ln(Z)/44
func (t *Thermodynamics) PerBondFreeEnergy(Z float64) float64 {
	return -kB * t.Temperature * math.Log(Z) / TotalBonds
}

// GibbsFreeEnergy: ΔG°_total = Σ ΔG°_bond(n) + ΔG°_coop + ΔG°_10.6
func (t *Thermodynamics) GibbsFreeEnergy(bondEnergies []float64, GCoop float64) float64 {
	GTotal := 0.0
	for _, E := range bondEnergies {
		GTotal += E
	}
	GTotal += GCoop
	// ΔG°_10.6 = -RT · ln(10.6/10.4) = -0.019 RT
	G10_6 := -R * t.Temperature * math.Log(BasePairsPerTurn/10.4)
	GTotal += G10_6
	return GTotal
}

// MeltingTemperature: T_m = ΔH°_44/(ΔS°_44 + R·ln(C_t)) × (10.6/10.4)
func (t *Thermodynamics) MeltingTemperature(deltaH, deltaS, Ct float64) float64 {
	correction := BasePairsPerTurn / 10.4
	return (deltaH / (deltaS + R*math.Log(Ct))) * correction
}

// BoltzmannProbability for state with energy E
func (t *Thermodynamics) BoltzmannProbability(E, Z float64) float64 {
	return math.Exp(-E/(kB*t.Temperature)) / Z
}

// =============================================================================
// XII. KINETIC RATE EQUATIONS
// =============================================================================

// KineticRates handles bond formation and breaking rates
type KineticRates struct {
	// Rate constants
	Af   float64 // Formation pre-factor
	Ab   float64 // Breaking pre-factor
	E0f  float64 // Base formation activation energy
	E0b  float64 // Base breaking activation energy
	DEf  float64 // Formation energy variation
	DEb  float64 // Breaking energy variation
	E106 float64 // 10.6 modulation energy

	// Modulation amplitudes
	MuF float64 // Formation modulation
	MuB float64 // Breaking modulation

	// Cooperativity
	SigmaF float64 // Formation cooperativity
	SigmaB float64 // Breaking cooperativity

	// Coupling
	Alpha float64 // Diffusion coefficient
	Beta  float64 // Inter-type coupling
	Gamma float64 // Cross-helix coupling
	Delta float64 // Cross-helix secondary

	Temperature float64
}

const kBRate = 0.001987 // kcal/(mol·K)

// NewKineticRates creates a new KineticRates with default parameters
func NewKineticRates() *KineticRates {
	return &KineticRates{
		Af:          1e6,
		Ab:          1e4,
		E0f:         10.0,
		E0b:         20.0,
		DEf:         2.0,
		DEb:         2.0,
		E106:        1.0,
		MuF:         0.2,
		MuB:         0.2,
		SigmaF:      0.5,
		SigmaB:      0.3,
		Alpha:       0.1,
		Beta:        0.05,
		Gamma:       0.02,
		Delta:       0.01,
		Temperature: 300.0,
	}
}

// ActivationEnergyFormation: E_a^f(n) = E_0^f + ΔE_f·sin²(πn/44) + E_10.6·cos(2πn/10.6)
func (kr *KineticRates) ActivationEnergyFormation(n int) float64 {
	sinTerm := math.Sin(math.Pi * float64(n) / TotalBonds)
	cosTerm := math.Cos(2.0 * math.Pi * float64(n) / BasePairsPerTurn)
	return kr.E0f + kr.DEf*sinTerm*sinTerm + kr.E106*cosTerm
}

// ActivationEnergyBreaking: E_a^b(n) = E_0^b - ΔE_b·sin²(πn/44) - E_10.6·cos(2πn/10.6)
func (kr *KineticRates) ActivationEnergyBreaking(n int) float64 {
	sinTerm := math.Sin(math.Pi * float64(n) / TotalBonds)
	cosTerm := math.Cos(2.0 * math.Pi * float64(n) / BasePairsPerTurn)
	return kr.E0b - kr.DEb*sinTerm*sinTerm - kr.E106*cosTerm
}

// BaseFormationRate: k_f^(0)(n) = A_f · exp(-E_a^f(n)/k_B T)
func (kr *KineticRates) BaseFormationRate(n int) float64 {
	Ea := kr.ActivationEnergyFormation(n)
	return kr.Af * math.Exp(-Ea/(kBRate*kr.Temperature))
}

// BaseBreakingRate: k_b^(0)(n) = A_b · exp(-E_a^b(n)/k_B T)
func (kr *KineticRates) BaseBreakingRate(n int) float64 {
	Ea := kr.ActivationEnergyBreaking(n)
	return kr.Ab * math.Exp(-Ea/(kBRate*kr.Temperature))
}

// FormationRate: k_f(n,t) = k_f^(0)(n) · [1 + μ_f·cos(2πn/10.6 + φ(t))]
func (kr *KineticRates) FormationRate(n int, phiT float64) float64 {
	k0 := kr.BaseFormationRate(n)
	modulation := 1.0 + kr.MuF*math.Cos(2.0*math.Pi*float64(n)/BasePairsPerTurn+phiT)
	return k0 * modulation
}

// BreakingRate: k_b(n,t) = k_b^(0)(n) · [1 - μ_b·cos(2πn/10.6 + φ(t))]
func (kr *KineticRates) BreakingRate(n int, phiT float64) float64 {
	k0 := kr.BaseBreakingRate(n)
	modulation := 1.0 - kr.MuB*math.Cos(2.0*math.Pi*float64(n)/BasePairsPerTurn+phiT)
	return k0 * math.Max(0.0, modulation)
}

// GetNeighbors returns neighbor bond indices for bond n
func (kr *KineticRates) GetNeighbors(n int) []int {
	var neighbors []int

	if n <= IntraStrandBonds {
		// N_intra(k) = {k-1, k+1} mod 16
		k := n
		neighbors = append(neighbors, ((k-2+NumStrands)%NumStrands)+1)
		neighbors = append(neighbors, (k%NumStrands)+1)
	} else if n <= IntraStrandBonds+InterStrandBonds {
		// N_inter(k) = {k, k+16}
		k := n - IntraStrandBonds
		neighbors = append(neighbors, k)
		neighbors = append(neighbors, k+IntraStrandBonds)
	} else {
		// N_cross(m) = {m-1, m+1, m+4}
		m := n - IntraStrandBonds - InterStrandBonds
		if m > 1 {
			neighbors = append(neighbors, n-1)
		}
		if m < CrossHelixBonds {
			neighbors = append(neighbors, n+1)
		}
		if m+4 <= CrossHelixBonds {
			neighbors = append(neighbors, n+4)
		}
	}

	return neighbors
}

// CooperativeFormationRate: k_f^coop(n,t) = k_f(n,t) · Π[1 + σ_f·b_m(t)]
func (kr *KineticRates) CooperativeFormationRate(n int, phiT float64, occupations []float64) float64 {
	k := kr.FormationRate(n, phiT)
	neighbors := kr.GetNeighbors(n)
	for _, m := range neighbors {
		if m > 0 && m <= TotalBonds {
			k *= (1.0 + kr.SigmaF*occupations[m-1])
		}
	}
	return k
}

// CooperativeBreakingRate: k_b^coop(n,t) = k_b(n,t) · Π[1 - σ_b·b_m(t)]
func (kr *KineticRates) CooperativeBreakingRate(n int, phiT float64, occupations []float64) float64 {
	k := kr.BreakingRate(n, phiT)
	neighbors := kr.GetNeighbors(n)
	for _, m := range neighbors {
		if m > 0 && m <= TotalBonds {
			k *= (1.0 - kr.SigmaB*occupations[m-1])
		}
	}
	return math.Max(0.0, k)
}

// EquilibriumOccupation: b_n^eq = k_f(n)/(k_f(n) + k_b(n))
func (kr *KineticRates) EquilibriumOccupation(n int, phiT float64) float64 {
	kf := kr.FormationRate(n, phiT)
	kb := kr.BreakingRate(n, phiT)
	return kf / (kf + kb)
}

// SlowRelaxationTime: τ_slow = 1/(k_f^min + k_b^min) · [1 + 43·σ_eff]
func (kr *KineticRates) SlowRelaxationTime() float64 {
	kMin := math.Min(kr.BaseFormationRate(22), kr.BaseBreakingRate(22))
	sigmaEff := (kr.SigmaF + kr.SigmaB) / 2.0
	return (1.0 / kMin) * (1.0 + 43.0*sigmaEff)
}

// FastRelaxationTime: τ_fast = 1/(k_f^max + k_b^max)
func (kr *KineticRates) FastRelaxationTime() float64 {
	kMax := kr.BaseFormationRate(1) + kr.BaseBreakingRate(1)
	return 1.0 / kMax
}

// =============================================================================
// BOND DYNAMICS
// =============================================================================

// BondDynamics manages the dynamics of all 44 bonds
type BondDynamics struct {
	B      [TotalBonds]float64 // Bond occupations
	Rates  *KineticRates
	Phase  float64 // φ_10.6(t)
	Time   float64
	Alpha  float64 // Diffusion
	Beta   float64 // Inter-coupling
	Gamma  float64 // Cross-coupling
	Delta  float64 // Cross secondary
}

// NewBondDynamics creates a new BondDynamics
func NewBondDynamics() *BondDynamics {
	bd := &BondDynamics{
		Rates: NewKineticRates(),
		Alpha: 0.1,
		Beta:  0.05,
		Gamma: 0.02,
		Delta: 0.01,
	}
	// Initialize at equilibrium
	for n := 0; n < TotalBonds; n++ {
		bd.B[n] = bd.Rates.EquilibriumOccupation(n+1, 0.0)
	}
	return bd
}

// SingleBondDerivative: db_n/dt = k_f·(1-b_n) - k_b·b_n
func (bd *BondDynamics) SingleBondDerivative(n int) float64 {
	occupations := bd.B[:]
	kf := bd.Rates.CooperativeFormationRate(n, bd.Phase, occupations)
	kb := bd.Rates.CooperativeBreakingRate(n, bd.Phase, occupations)
	return kf*(1.0-bd.B[n-1]) - kb*bd.B[n-1]
}

// ComputeDerivatives computes all bond derivatives
func (bd *BondDynamics) ComputeDerivatives() [TotalBonds]float64 {
	var db [TotalBonds]float64

	// Intra-strand (n = 1..16)
	for n := 1; n <= IntraStrandBonds; n++ {
		idx := n - 1
		prev := ((n - 2 + NumStrands) % NumStrands)
		next := (n % NumStrands)

		db[idx] = bd.SingleBondDerivative(n)
		// Diffusion: α·(b_{n-1} + b_{n+1} - 2b_n)
		db[idx] += bd.Alpha * (bd.B[prev] + bd.B[next] - 2.0*bd.B[idx])
		// Inter-coupling
		if n+IntraStrandBonds <= TotalBonds {
			db[idx] += bd.Beta * math.Sin(2.0*math.Pi*float64(n)/BasePairsPerTurn) *
				(bd.B[n+IntraStrandBonds-1] - bd.B[idx])
		}
	}

	// Inter-strand (n = 17..32)
	for n := IntraStrandBonds + 1; n <= IntraStrandBonds+InterStrandBonds; n++ {
		idx := n - 1
		prev := idx - 1
		next := idx + 1
		if next >= IntraStrandBonds+InterStrandBonds {
			next = IntraStrandBonds
		}

		db[idx] = bd.SingleBondDerivative(n)
		// Diffusion
		db[idx] += bd.Alpha * (bd.B[prev] + bd.B[next] - 2.0*bd.B[idx])
		// Cross-coupling
		intraIdx := n - IntraStrandBonds - 1
		db[idx] += bd.Gamma * math.Cos(2.0*math.Pi*float64(n)/BasePairsPerTurn) *
			(bd.B[intraIdx] - bd.B[idx])
	}

	// Cross-helix (n = 33..44)
	for n := IntraStrandBonds + InterStrandBonds + 1; n <= TotalBonds; n++ {
		idx := n - 1
		m := n - IntraStrandBonds - InterStrandBonds

		db[idx] = bd.SingleBondDerivative(n)

		// Diffusion with boundary handling
		if idx > IntraStrandBonds+InterStrandBonds {
			prev := idx - 1
			next := idx + 1
			if next >= TotalBonds {
				next = IntraStrandBonds + InterStrandBonds
			}
			db[idx] += bd.Alpha * (bd.B[prev] + bd.B[next] - 2.0*bd.B[idx])
		}

		// Cross-helix coupling
		db[idx] += bd.Delta * math.Cos(2.0*math.Pi*float64(m)*4.0/BasePairsPerTurn)
	}

	return db
}

// Step integrates one step using RK4
func (bd *BondDynamics) Step(dt float64) {
	// k1
	k1 := bd.ComputeDerivatives()

	// k2
	var bTemp [TotalBonds]float64
	for i := 0; i < TotalBonds; i++ {
		bTemp[i] = bd.B[i] + 0.5*dt*k1[i]
	}
	bBackup := bd.B
	bd.B = bTemp
	k2 := bd.ComputeDerivatives()

	// k3
	for i := 0; i < TotalBonds; i++ {
		bTemp[i] = bBackup[i] + 0.5*dt*k2[i]
	}
	bd.B = bTemp
	k3 := bd.ComputeDerivatives()

	// k4
	for i := 0; i < TotalBonds; i++ {
		bTemp[i] = bBackup[i] + dt*k3[i]
	}
	bd.B = bTemp
	k4 := bd.ComputeDerivatives()

	// Combine
	for i := 0; i < TotalBonds; i++ {
		bd.B[i] = bBackup[i] + (dt/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i])
		if bd.B[i] < 0.0 {
			bd.B[i] = 0.0
		} else if bd.B[i] > 1.0 {
			bd.B[i] = 1.0
		}
	}

	// Update phase
	bd.Time += dt
	bd.Phase = math.Mod(2.0*math.Pi*bd.Time/BasePairsPerTurn, 2.0*math.Pi)
}

// SteadyStateBonds returns time-averaged steady state
func (bd *BondDynamics) SteadyStateBonds() float64 {
	sum := 0.0
	for n := 1; n <= TotalBonds; n++ {
		sum += bd.Rates.EquilibriumOccupation(n, 0.0)
	}
	return sum
}

// SteadyStateOccupation returns θ* = B*/44
func (bd *BondDynamics) SteadyStateOccupation() float64 {
	return bd.SteadyStateBonds() / TotalBonds
}

// TotalOccupation returns current total occupation
func (bd *BondDynamics) TotalOccupation() float64 {
	sum := 0.0
	for _, b := range bd.B {
		sum += b
	}
	return sum
}