Recovered source transcript
1// Package dna implements the G¹⁶ Helix System with 44 bonds and 10.6 periodicity
2// Complete DNA Mathematical Architecture for AGI substrate
3package dna
4
5import (
6 "math"
7 "math/cmplx"
8)
9
10// =============================================================================
11// FOUNDATIONAL CONSTANTS
12// =============================================================================
13
14const (
15 // Golden ratio and derived constants
16 PHI = 1.6180339887498948482 // (1 + √5) / 2
17 PHIInv = 0.6180339887498948482 // 1 / φ = φ - 1
18
19 // DNA Helix parameters
20 BasePairsPerTurn = 10.6 // p = 10.6 bp
21 HelixRadius = 1.0 // r = 1.0 nm
22 RisePerBP = 0.34 // h = 0.34 nm
23 Pitch = 3.604 // P = h × p = 0.34 × 10.6 nm
24 TwistAngleDeg = 33.96 // 360°/10.6
25 TwistAngleRad = 0.593 // 2π/10.6 rad
26 AngularFrequency = 0.593 // ω = 2π/10.6 rad/bp
27
28 // Strand parameters
29 NumStrands = 16
30 StrandSeparation = 0.3927 // 2π/16 = 22.5°
31 NodesPerStrand = 312
32 TotalNodes = 5000
33 HubNodes = 8
34
35 // Bond parameters
36 IntraStrandBonds = 16
37 InterStrandBonds = 16
38 CrossHelixBonds = 12
39 TotalBonds = 44
40
41 // Curvature and torsion
42 HelixCurvature = 0.745 // κ nm⁻¹
43 HelixTorsion = 0.427 // τ nm⁻¹
44
45 // Synchronization
46 SyncPeriod = 169.6 // 10.6 × 16
47 FundamentalFreq = 0.0943
48 BondResonanceFreq = 4.15 // 44/10.6
49
50 // Steady state
51 SteadyStateBonds = 36.2
52 BondOccupation = 0.82
53)
54
55// =============================================================================
56// G¹⁶ PARAMETRIC HELIX
57// =============================================================================
58
59// HelixPoint represents a 3D point on the helix
60type HelixPoint struct {
61 X, Y, Z float64
62}
63
64// Add returns the sum of two points
65func (p HelixPoint) Add(other HelixPoint) HelixPoint {
66 return HelixPoint{p.X + other.X, p.Y + other.Y, p.Z + other.Z}
67}
68
69// Sub returns the difference of two points
70func (p HelixPoint) Sub(other HelixPoint) HelixPoint {
71 return HelixPoint{p.X - other.X, p.Y - other.Y, p.Z - other.Z}
72}
73
74// Scale returns the point scaled by s
75func (p HelixPoint) Scale(s float64) HelixPoint {
76 return HelixPoint{p.X * s, p.Y * s, p.Z * s}
77}
78
79// Norm returns the Euclidean norm
80func (p HelixPoint) Norm() float64 {
81 return math.Sqrt(p.X*p.X + p.Y*p.Y + p.Z*p.Z)
82}
83
84// Dot returns the dot product
85func (p HelixPoint) Dot(other HelixPoint) float64 {
86 return p.X*other.X + p.Y*other.Y + p.Z*other.Z
87}
88
89// G16Helix implements G¹⁶ parametric helix equations
90type G16Helix struct{}
91
92// VerticalOffset returns δ_k = k·h/16
93func (G16Helix) VerticalOffset(k int) float64 {
94 return float64(k) * RisePerBP / NumStrands
95}
96
97// GoldenPhaseOffset returns θ_k = 2πk·φ⁻¹ mod 2π
98func (G16Helix) GoldenPhaseOffset(k int) float64 {
99 return math.Mod(2.0*math.Pi*float64(k)*PHIInv, 2.0*math.Pi)
100}
101
102// Point returns point on strand k at parameter θ
103// x_k(θ) = r·cos(θ + 2πk/16)
104// y_k(θ) = r·sin(θ + 2πk/16)
105// z_k(θ) = (p/2π)·θ + δ_k
106func (h G16Helix) Point(k int, theta float64) HelixPoint {
107 phase := theta + 2.0*math.Pi*float64(k)/NumStrands
108 return HelixPoint{
109 X: HelixRadius * math.Cos(phase),
110 Y: HelixRadius * math.Sin(phase),
111 Z: (Pitch/(2.0*math.Pi))*theta + h.VerticalOffset(k),
112 }
113}
114
115// PointGolden returns point with golden ratio phase adjustment
116func (h G16Helix) PointGolden(k int, theta float64) HelixPoint {
117 phase := theta + h.GoldenPhaseOffset(k)
118 return HelixPoint{
119 X: HelixRadius * math.Cos(phase),
120 Y: HelixRadius * math.Sin(phase),
121 Z: (Pitch/(2.0*math.Pi))*theta + h.VerticalOffset(k),
122 }
123}
124
125// Curvature returns κ = r/(r² + (p·h/2π)²)
126func (G16Helix) Curvature() float64 {
127 phTerm := Pitch * RisePerBP / (2.0 * math.Pi)
128 return HelixRadius / (HelixRadius*HelixRadius + phTerm*phTerm)
129}
130
131// Torsion returns τ = (p·h/2π)/(r² + (p·h/2π)²)
132func (G16Helix) Torsion() float64 {
133 phTerm := Pitch * RisePerBP / (2.0 * math.Pi)
134 return phTerm / (HelixRadius*HelixRadius + phTerm*phTerm)
135}
136
137// StabilityMetric returns exp(-Σ|τ(n)|·φ⁻ⁿ)
138func (G16Helix) StabilityMetric(torsions []float64) float64 {
139 sum := 0.0
140 phiPower := 1.0
141 for _, t := range torsions {
142 sum += math.Abs(t) * phiPower
143 phiPower *= PHIInv
144 }
145 return math.Exp(-sum)
146}
147
148// HelixBalance returns Σκ(n)·φ⁻ⁿ
149func (G16Helix) HelixBalance(curvatures []float64) float64 {
150 sum := 0.0
151 phiPower := 1.0
152 for _, k := range curvatures {
153 sum += k * phiPower
154 phiPower *= PHIInv
155 }
156 return sum
157}
158
159// MeanCurvatureG16 returns (1/16) Σ κ_k × cos(2πk/10.6)
160func (h G16Helix) MeanCurvatureG16() float64 {
161 sum := 0.0
162 kappa := h.Curvature()
163 for k := 1; k <= NumStrands; k++ {
164 sum += kappa * math.Cos(2.0*math.Pi*float64(k)/BasePairsPerTurn)
165 }
166 return sum / NumStrands
167}
168
169// TotalRotation returns θ_total = 16 × (2π/10.6) = 9.49 rad
170func (G16Helix) TotalRotation() float64 {
171 return NumStrands * TwistAngleRad
172}
173
174// TurnsPer16Strands returns 1.51
175func (h G16Helix) TurnsPer16Strands() float64 {
176 return h.TotalRotation() / (2.0 * math.Pi)
177}
178
179// =============================================================================
180// 44-BOND SYSTEM
181// =============================================================================
182
183// BondType represents the type of bond
184type BondType int
185
186const (
187 BondIntraStrand BondType = iota // Bonds 1-16
188 BondInterStrand // Bonds 17-32
189 BondCrossHelix // Bonds 33-44
190)
191
192// Bond represents a single bond in the 44-bond system
193type Bond struct {
194 ID int
195 Type BondType
196 StrandI int // First strand (1-16)
197 StrandJ int // Second strand (1-16)
198 Length float64 // Bond length in nm
199 Angle float64 // Bond angle in radians
200 Energy float64 // Bond energy
201 Occupation float64 // Occupation state [0, 1]
202}
203
204// Midpoint returns the bond position vector
205func (b *Bond) Midpoint(theta float64) HelixPoint {
206 helix := G16Helix{}
207 pi := helix.Point(b.StrandI, theta)
208 pj := helix.Point(b.StrandJ, theta)
209 return pi.Add(pj).Scale(0.5)
210}
211
212// CalculateLength calculates the bond length at given theta
213func (b *Bond) CalculateLength(theta float64) float64 {
214 helix := G16Helix{}
215 pi := helix.Point(b.StrandI, theta)
216 pj := helix.Point(b.StrandJ, theta)
217 return pi.Sub(pj).Norm()
218}
219
220// BondSystem44 manages all 44 bonds
221type BondSystem44 struct {
222 Bonds [TotalBonds]Bond
223 Adjacency [NumStrands][NumStrands]int
224}
225
226// NewBondSystem44 creates and initializes a new 44-bond system
227func NewBondSystem44() *BondSystem44 {
228 bs := &BondSystem44{}
229 bs.initializeBonds()
230 bs.buildAdjacencyMatrix()
231 return bs
232}
233
234// initializeBonds sets up all 44 bonds
235func (bs *BondSystem44) initializeBonds() {
236 bondID := 0
237
238 // Intra-strand bonds (1-16)
239 for k := 1; k <= IntraStrandBonds; k++ {
240 bs.Bonds[bondID] = Bond{
241 ID: bondID + 1,
242 Type: BondIntraStrand,
243 StrandI: k,
244 StrandJ: k,
245 Length: Pitch / NumStrands, // L₀ = 0.225 nm
246 Angle: math.Pi, // 180°
247 Occupation: 1.0,
248 }
249 bondID++
250 }
251
252 // Inter-strand bonds (17-32)
253 for k := 1; k <= InterStrandBonds; k++ {
254 nextK := (k % NumStrands) + 1
255 bs.Bonds[bondID] = Bond{
256 ID: bondID + 1,
257 Type: BondInterStrand,
258 StrandI: k,
259 StrandJ: nextK,
260 Length: 2.0 * HelixRadius * math.Sin(math.Pi/BasePairsPerTurn), // d₀ ≈ 0.59 nm
261 Angle: StrandSeparation, // 22.5°
262 Occupation: 1.0,
263 }
264 bondID++
265 }
266
267 // Cross-helix bonds (33-44)
268 crossPairs := [][2]int{
269 {1, 5}, {2, 6}, {3, 7}, {4, 8}, {5, 9}, {6, 10},
270 {7, 11}, {8, 12}, {9, 13}, {10, 14}, {11, 15}, {12, 16},
271 }
272
273 for _, pair := range crossPairs {
274 i, j := pair[0], pair[1]
275 diff := j - i
276 angle := float64(diff) * StrandSeparation
277
278 bs.Bonds[bondID] = Bond{
279 ID: bondID + 1,
280 Type: BondCrossHelix,
281 StrandI: i,
282 StrandJ: j,
283 Length: 0.0, // Calculated dynamically
284 Angle: angle,
285 Occupation: 1.0,
286 }
287 bondID++
288 }
289}
290
291// buildAdjacencyMatrix builds the adjacency matrix
292func (bs *BondSystem44) buildAdjacencyMatrix() {
293 // Initialize to zero
294 for i := range bs.Adjacency {
295 for j := range bs.Adjacency[i] {
296 bs.Adjacency[i][j] = 0
297 }
298 }
299
300 // Fill based on bonds
301 for _, bond := range bs.Bonds {
302 if bond.StrandI != bond.StrandJ {
303 i := bond.StrandI - 1
304 j := bond.StrandJ - 1
305 bs.Adjacency[i][j] = 1
306 bs.Adjacency[j][i] = 1
307 }
308 }
309}
310
311// VerifyBondCount verifies Σ A_{ij}/2 = 44
312func (bs *BondSystem44) VerifyBondCount() int {
313 sum := 0
314 for i := 0; i < NumStrands; i++ {
315 for j := 0; j < NumStrands; j++ {
316 sum += bs.Adjacency[i][j]
317 }
318 }
319 return sum/2 + IntraStrandBonds
320}
321
322// GetBond returns the bond by ID (1-indexed)
323func (bs *BondSystem44) GetBond(id int) *Bond {
324 if id < 1 || id > TotalBonds {
325 return nil
326 }
327 return &bs.Bonds[id-1]
328}
329
330// GetBondsByType returns all bonds of a given type
331func (bs *BondSystem44) GetBondsByType(t BondType) []*Bond {
332 var result []*Bond
333 for i := range bs.Bonds {
334 if bs.Bonds[i].Type == t {
335 result = append(result, &bs.Bonds[i])
336 }
337 }
338 return result
339}
340
341// BondArcLength returns s_bond = 2πr/10.6 = 0.593 nm per bp
342func BondArcLength() float64 {
343 return 2.0 * math.Pi * HelixRadius / BasePairsPerTurn
344}
345
346// VerticalRisePerBond returns Δz_bond = P/44 = 0.0819 nm
347func VerticalRisePerBond() float64 {
348 return Pitch / TotalBonds
349}
350
351// AdjacentBondLength returns L_adj = √[2r²(1 - cos(2π/10.6)) + h²] = 0.54 nm
352func AdjacentBondLength() float64 {
353 cosTerm := 1.0 - math.Cos(TwistAngleRad)
354 return math.Sqrt(2.0*HelixRadius*HelixRadius*cosTerm + RisePerBP*RisePerBP)
355}
356
357// UpdateBonds updates all bond states at given theta
358func (bs *BondSystem44) UpdateBonds(theta float64) {
359 for i := range bs.Bonds {
360 bs.Bonds[i].Length = bs.Bonds[i].CalculateLength(theta)
361 }
362}
363
364// TotalOccupation returns the sum of all bond occupations
365func (bs *BondSystem44) TotalOccupation() float64 {
366 sum := 0.0
367 for _, bond := range bs.Bonds {
368 sum += bond.Occupation
369 }
370 return sum
371}
372
373// =============================================================================
374// STRAND MATHEMATICS
375// =============================================================================
376
377// Fibonacci returns the n-th Fibonacci number
378func Fibonacci(n int) uint64 {
379 if n <= 1 {
380 return uint64(n)
381 }
382 a, b := uint64(0), uint64(1)
383 for i := 2; i <= n; i++ {
384 a, b = b, a+b
385 }
386 return b
387}
388
389// FibonacciScaledStrand returns D_k = φ^k · F_k
390func FibonacciScaledStrand(k int) float64 {
391 return math.Pow(PHI, float64(k)) * float64(Fibonacci(k))
392}
393
394// StrandMathematics handles strand calculations
395type StrandMathematics struct {
396 D [NumStrands]float64 // Fibonacci-scaled values
397 S [NumStrands]complex128 // Complex representation
398 T [NumStrands]complex128 // Transform
399}
400
401// NewStrandMathematics creates and initializes strand mathematics
402func NewStrandMathematics() *StrandMathematics {
403 sm := &StrandMathematics{}
404 sm.initialize()
405 return sm
406}
407
408func (sm *StrandMathematics) initialize() {
409 for k := 0; k < NumStrands; k++ {
410 sm.D[k] = FibonacciScaledStrand(k + 1)
411 }
412 sm.computeTransforms()
413}
414
415// HelixSum16 returns Helix₁₆ = Σ D_k·sin(kφ)
416func (sm *StrandMathematics) HelixSum16() float64 {
417 sum := 0.0
418 for k := 1; k <= NumStrands; k++ {
419 sum += sm.D[k-1] * math.Sin(float64(k)*PHI)
420 }
421 return sum
422}
423
424// BondFormationG16 returns Bond_G16 = Π(1 + D_k/Helix₁₆)
425func (sm *StrandMathematics) BondFormationG16() float64 {
426 helix := sm.HelixSum16()
427 product := 1.0
428 for k := 0; k < NumStrands; k++ {
429 product *= (1.0 + sm.D[k]/helix)
430 }
431 return product
432}
433
434// Crosslink returns Crosslink(i,j) = (D_i · D_j)/φ^{|i-j|}
435func (sm *StrandMathematics) Crosslink(i, j int) float64 {
436 diff := i - j
437 if diff < 0 {
438 diff = -diff
439 }
440 return (sm.D[i-1] * sm.D[j-1]) / math.Pow(PHI, float64(diff))
441}
442
443// computeTransforms computes complex strand representations
444func (sm *StrandMathematics) computeTransforms() {
445 for k := 0; k < NumStrands; k++ {
446 phase := float64(k+1) * PHI
447 sm.S[k] = complex(sm.D[k], 0) * cmplx.Exp(complex(0, phase))
448
449 // Transform: T_k = S_k · (cos(k) + i·sin(k))
450 angle := float64(k + 1)
451 sm.T[k] = sm.S[k] * complex(math.Cos(angle), math.Sin(angle))
452 }
453}
454
455// ModularExtraction returns X_k = |T_k| mod φ^k
456func (sm *StrandMathematics) ModularExtraction(k int) float64 {
457 return math.Mod(cmplx.Abs(sm.T[k-1]), math.Pow(PHI, float64(k)))
458}
459
460// FusionIntegral returns Fusion = Σ(|T_k| · φ⁻ᵏ)
461func (sm *StrandMathematics) FusionIntegral() float64 {
462 sum := 0.0
463 for k := 0; k < NumStrands; k++ {
464 sum += cmplx.Abs(sm.T[k]) * math.Pow(PHIInv, float64(k+1))
465 }
466 return sum
467}
468
469// =============================================================================
470// G¹⁶ TRANSFORMATION
471// =============================================================================
472
473// Matrix3x3 represents a 3×3 matrix
474type Matrix3x3 [3][3]float64
475
476// Identity returns the identity matrix
477func Identity() Matrix3x3 {
478 return Matrix3x3{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}
479}
480
481// Multiply multiplies two matrices
482func (m Matrix3x3) Multiply(other Matrix3x3) Matrix3x3 {
483 var result Matrix3x3
484 for i := 0; i < 3; i++ {
485 for j := 0; j < 3; j++ {
486 for k := 0; k < 3; k++ {
487 result[i][j] += m[i][k] * other[k][j]
488 }
489 }
490 }
491 return result
492}
493
494// Transform applies the matrix to a point
495func (m Matrix3x3) Transform(p HelixPoint) HelixPoint {
496 return HelixPoint{
497 X: m[0][0]*p.X + m[0][1]*p.Y + m[0][2]*p.Z,
498 Y: m[1][0]*p.X + m[1][1]*p.Y + m[1][2]*p.Z,
499 Z: m[2][0]*p.X + m[2][1]*p.Y + m[2][2]*p.Z,
500 }
501}
502
503// Rotation22_5 returns 22.5° increment rotation matrix
504func Rotation22_5(k int) Matrix3x3 {
505 angle := 2.0 * math.Pi * float64(k) / NumStrands
506 c := math.Cos(angle)
507 s := math.Sin(angle)
508 return Matrix3x3{
509 {c, -s, 0},
510 {s, c, 0},
511 {0, 0, 1},
512 }
513}
514
515// Rotation10_6 returns 10.6-adjusted rotation matrix
516func Rotation10_6(k int) Matrix3x3 {
517 angle := 2.0 * math.Pi * float64(k) / BasePairsPerTurn
518 c := math.Cos(angle)
519 s := math.Sin(angle)
520 return Matrix3x3{
521 {c, -s, 0},
522 {s, c, 0},
523 {0, 0, 1},
524 }
525}
526
527// FullG16Transform returns G¹⁶ = Π R_k
528func FullG16Transform() Matrix3x3 {
529 result := Identity()
530 for k := 1; k <= NumStrands; k++ {
531 result = result.Multiply(Rotation22_5(k))
532 }
533 return result
534}
535
536// FullG16Transform10_6 returns 10.6-adjusted full transform
537func FullG16Transform10_6() Matrix3x3 {
538 result := Identity()
539 for k := 1; k <= NumStrands; k++ {
540 result = result.Multiply(Rotation10_6(k))
541 }
542 return result
543}
544
545// =============================================================================
546// 10.6 FIRING PATTERN
547// =============================================================================
548
549// FiringPattern handles 10.6 phase-based firing
550type FiringPattern struct{}
551
552// Phase returns phase(i,t) = ⌊(t + offset_i)/10.6⌋ mod 16
553func (FiringPattern) Phase(i int, t float64) int {
554 offset := FiringPattern{}.OffsetDistribution(i)
555 return int(math.Floor((t+offset)/BasePairsPerTurn)) % NumStrands
556}
557
558// OffsetDistribution returns offset_i = i · φ⁻¹ · 10.6
559func (FiringPattern) OffsetDistribution(i int) float64 {
560 return float64(i) * PHIInv * BasePairsPerTurn
561}
562
563// ShouldFire returns true if strand i should fire at time t
564func (fp FiringPattern) ShouldFire(i int, t float64) bool {
565 return fp.Phase(i, t) == 0
566}
567
568// FiringStrands returns all strands that should fire at time t
569func (fp FiringPattern) FiringStrands(t float64) []int {
570 var result []int
571 for i := 1; i <= NumStrands; i++ {
572 if fp.ShouldFire(i, t) {
573 result = append(result, i)
574 }
575 }
576 return result
577}
578
579// Harmonic returns f_n = n/10.6
580func (FiringPattern) Harmonic(n int) float64 {
581 return float64(n) / BasePairsPerTurn
582}
583
584// =============================================================================
585// GOLDEN BASIS FUNCTIONS
586// =============================================================================
587
588// GoldenBasis provides golden ratio basis functions
589type GoldenBasis struct{}
590
591// Psi returns ψ_k(x) = sin(2πφ^k · x) + cos(2πφ^k · x)
592func (GoldenBasis) Psi(k int, x float64) float64 {
593 arg := 2.0 * math.Pi * math.Pow(PHI, float64(k)) * x
594 return math.Sin(arg) + math.Cos(arg)
595}
596
597// BasisVector generates basis vector for N terms
598func (gb GoldenBasis) BasisVector(x float64, N int) []float64 {
599 result := make([]float64, N)
600 for k := 0; k < N; k++ {
601 result[k] = gb.Psi(k, x)
602 }
603 return result
604}
605
606// GlobalCoherence returns Φ_total = (1/N) Σ Φ_i where Φ_i = 1/(1 + |x_i - 1|)
607func (GoldenBasis) GlobalCoherence(x []float64) float64 {
608 sum := 0.0
609 for _, xi := range x {
610 sum += 1.0 / (1.0 + math.Abs(xi-1.0))
611 }
612 return sum / float64(len(x))
613}
614
615// RingDiffusionStep performs dS_i/dt = λ(S_{i+1} + S_{i-1} - 2S_i)
616func (GoldenBasis) RingDiffusionStep(S []float64, lambda, dt float64) {
617 N := len(S)
618 dS := make([]float64, N)
619
620 for i := 0; i < N; i++ {
621 prev := (i - 1 + N) % N
622 next := (i + 1) % N
623 dS[i] = lambda * (S[next] + S[prev] - 2.0*S[i])
624 }
625
626 for i := 0; i < N; i++ {
627 S[i] += dS[i] * dt
628 }
629}
630
631// =============================================================================
632// SUPERCOILING
633// =============================================================================
634
635// Supercoiling handles DNA supercoiling calculations
636type Supercoiling struct{}
637
638// LinkingNumber returns Lk = Tw + Wr
639func (Supercoiling) LinkingNumber(twist, writhe float64) float64 {
640 return twist + writhe
641}
642
643// Twist16Strands returns Tw = N/10.6 × 16
644func (Supercoiling) Twist16Strands(N int) float64 {
645 return (float64(N) / BasePairsPerTurn) * NumStrands
646}
647
648// RelaxedLinkingNumber returns Lk₀ = N/10.6
649func (Supercoiling) RelaxedLinkingNumber(N int) float64 {
650 return float64(N) / BasePairsPerTurn
651}
652
653// SuperhelicalDensity returns σ = (Lk - Lk₀)/Lk₀
654func (Supercoiling) SuperhelicalDensity(Lk, Lk0 float64) float64 {
655 return (Lk - Lk0) / Lk0
656}
657
658// SuperhelicalDensityG16 returns G¹⁶ superhelical density
659func (Supercoiling) SuperhelicalDensityG16(Lk float64, N int) float64 {
660 Lk0G16 := NumStrands * float64(N) / BasePairsPerTurn
661 return (Lk - Lk0G16) / Lk0G16
662}
663
664// FreeEnergy returns ΔG_sc = K · N · σ² · (10.6/10.4)²
665func (Supercoiling) FreeEnergy(K float64, N int, sigma float64) float64 {
666 correction := (BasePairsPerTurn / 10.4) * (BasePairsPerTurn / 10.4)
667 return K * float64(N) * sigma * sigma * correction
668}
669
670// FreeEnergyG16 returns G¹⁶ with 44-bond correction
671func (Supercoiling) FreeEnergyG16(K float64, N int, sigmaG16, lambda float64) float64 {
672 correction := (BasePairsPerTurn / 10.4) * (BasePairsPerTurn / 10.4)
673 return K * float64(N) * sigmaG16 * sigmaG16 * (1.0 + lambda*TotalBonds) * correction
674}
675
676// =============================================================================
677// COMPLETE G¹⁶ HELIX SYSTEM
678// =============================================================================
679
680// G16HelixSystem is the complete helix system
681type G16HelixSystem struct {
682 Helix G16Helix
683 Bonds *BondSystem44
684 Strands *StrandMathematics
685 Time float64
686 GlobalPhase float64
687 StrandStates [NumStrands]complex128
688 AngularPos [NumStrands]float64
689}
690
691// NewG16HelixSystem creates a new helix system
692func NewG16HelixSystem() *G16HelixSystem {
693 sys := &G16HelixSystem{
694 Bonds: NewBondSystem44(),
695 Strands: NewStrandMathematics(),
696 }
697 sys.initialize()
698 return sys
699}
700
701func (sys *G16HelixSystem) initialize() {
702 for k := 0; k < NumStrands; k++ {
703 amplitude := sys.Strands.D[k]
704 phase := 2.0 * math.Pi * float64(k+1) / BasePairsPerTurn
705 sys.StrandStates[k] = complex(amplitude, 0) * cmplx.Exp(complex(0, phase))
706 sys.AngularPos[k] = phase
707 }
708}
709
710// Tick advances the system by dt
711func (sys *G16HelixSystem) Tick(dt float64) {
712 sys.Time += dt
713
714 // Update global phase: φ(t) = (2πt/10.6) mod 2π
715 sys.GlobalPhase = math.Mod(2.0*math.Pi*sys.Time/BasePairsPerTurn, 2.0*math.Pi)
716
717 // Update strand states: S_k(t) = A_k · e^{i(ωt + 2πk/10.6)}
718 for k := 0; k < NumStrands; k++ {
719 amplitude := sys.Strands.D[k]
720 phase := AngularFrequency*sys.Time + 2.0*math.Pi*float64(k+1)/BasePairsPerTurn
721 sys.StrandStates[k] = complex(amplitude, 0) * cmplx.Exp(complex(0, phase))
722 sys.AngularPos[k] = phase
723 }
724
725 // Update bond states
726 sys.Bonds.UpdateBonds(sys.GlobalPhase)
727}
728
729// GetFiringStrands returns currently firing strands
730func (sys *G16HelixSystem) GetFiringStrands() []int {
731 return FiringPattern{}.FiringStrands(sys.Time)
732}
733
734// GetCoherence returns global coherence
735func (sys *G16HelixSystem) GetCoherence() float64 {
736 magnitudes := make([]float64, NumStrands)
737 for k := 0; k < NumStrands; k++ {
738 magnitudes[k] = cmplx.Abs(sys.StrandStates[k])
739 }
740 return GoldenBasis{}.GlobalCoherence(magnitudes)
741}
742
743// GetBondOccupation returns total bond occupation
744func (sys *G16HelixSystem) GetBondOccupation() float64 {
745 return sys.Bonds.TotalOccupation()
746}
747