630 lines
28 KiB
Python
630 lines
28 KiB
Python
|
#!/usr/bin/env python3
|
|||
|
|
"""
|
|||
|
|
Nuclear Magic Number Analyzer
|
|||
|
|
Analyzes EM sweep data for signatures of nuclear shell structure
|
|||
|
|
in the coherence mode spectrum of the Resonance Engine lattice.
|
|||
|
|
|
|||
|
|
Six analyses:
|
|||
|
|
1. Coherence peak clustering — find & cluster coherence maxima
|
|||
|
|
2. Mode counting vs shell degeneracy — compare distinct mode counts to magic numbers
|
|||
|
|
3. Gap structure — coherence gap ratios vs nuclear shell gaps
|
|||
|
|
4. Omega-resolved shell occupancy — occupied-state count per omega slice
|
|||
|
|
5. 2D torus mode comparison — lattice mode degeneracies vs observed peaks
|
|||
|
|
6. GUE pair correlation — nearest-neighbor spacing vs Wigner surmise
|
|||
|
|
|
|||
|
|
Usage:
|
|||
|
|
python3 nuclear_magic_analyzer.py <sweep_csv>
|
|||
|
|
|
|||
|
|
Output saved to: ../results/nuclear_magic_analysis_<timestamp>.txt
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
import sys
|
|||
|
|
import os
|
|||
|
|
import numpy as np
|
|||
|
|
import pandas as pd
|
|||
|
|
from datetime import datetime
|
|||
|
|
from collections import Counter
|
|||
|
|
|
|||
|
|
# ── Nuclear physics constants ──────────────────────────────────────
|
|||
|
|
MAGIC_NUMBERS = [2, 8, 20, 28, 50, 82, 126]
|
|||
|
|
# Shell degeneracies (2j+1 for each filled subshell up to each magic closure)
|
|||
|
|
SHELL_DEGENERACIES = {
|
|||
|
|
2: [2], # 1s1/2
|
|||
|
|
8: [2, 4, 2], # 1s, 1p3/2, 1p1/2
|
|||
|
|
20: [2, 4, 2, 6, 4, 2], # sd shell
|
|||
|
|
28: [2, 4, 2, 6, 4, 2, 8], # f7/2
|
|||
|
|
50: [2, 4, 2, 6, 4, 2, 8, 6, 10, 4, 2],
|
|||
|
|
82: [2, 4, 2, 6, 4, 2, 8, 6, 10, 4, 2, 12, 8, 6, 4, 2],
|
|||
|
|
}
|
|||
|
|
# Gap ratios between successive magic numbers
|
|||
|
|
MAGIC_GAPS = np.diff(MAGIC_NUMBERS[:6]).astype(float)
|
|||
|
|
MAGIC_GAP_RATIOS = MAGIC_GAPS / MAGIC_GAPS[0] # normalized to first gap
|
|||
|
|
|
|||
|
|
|
|||
|
|
def load_sweep(csv_path):
|
|||
|
|
"""Load and validate sweep CSV."""
|
|||
|
|
df = pd.read_csv(csv_path)
|
|||
|
|
required = ['omega', 'khra_amp', 'gixx_amp', 'coherence', 'asymmetry', 'vorticity_mean']
|
|||
|
|
missing = [c for c in required if c not in df.columns]
|
|||
|
|
if missing:
|
|||
|
|
print(f"ERROR: Missing columns: {missing}")
|
|||
|
|
sys.exit(1)
|
|||
|
|
return df
|
|||
|
|
|
|||
|
|
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
# Analysis 1: Coherence Peak Clustering
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
def analysis_coherence_peaks(df, out):
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
out.append("ANALYSIS 1: Coherence Peak Clustering")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
|
|||
|
|
# Group by omega, find max coherence per omega slice
|
|||
|
|
omega_groups = df.groupby('omega')
|
|||
|
|
omega_vals = sorted(df['omega'].unique())
|
|||
|
|
|
|||
|
|
peak_data = []
|
|||
|
|
for omega in omega_vals:
|
|||
|
|
group = omega_groups.get_group(omega)
|
|||
|
|
idx_max = group['coherence'].idxmax()
|
|||
|
|
row = group.loc[idx_max]
|
|||
|
|
peak_data.append({
|
|||
|
|
'omega': omega,
|
|||
|
|
'coherence': row['coherence'],
|
|||
|
|
'khra_amp': row['khra_amp'],
|
|||
|
|
'gixx_amp': row['gixx_amp'],
|
|||
|
|
'asymmetry': row['asymmetry'],
|
|||
|
|
'vorticity': row['vorticity_mean'],
|
|||
|
|
})
|
|||
|
|
|
|||
|
|
peaks_df = pd.DataFrame(peak_data)
|
|||
|
|
coh_values = peaks_df['coherence'].values
|
|||
|
|
global_mean = coh_values.mean()
|
|||
|
|
global_std = coh_values.std()
|
|||
|
|
|
|||
|
|
out.append(f"\nPeak coherence per omega slice:")
|
|||
|
|
out.append(f" Mean: {global_mean:.6f} Std: {global_std:.6f}")
|
|||
|
|
out.append(f" Range: [{coh_values.min():.6f}, {coh_values.max():.6f}]")
|
|||
|
|
out.append("")
|
|||
|
|
|
|||
|
|
# Identify significant peaks (> mean + 1 sigma)
|
|||
|
|
threshold = global_mean + global_std
|
|||
|
|
strong_peaks = peaks_df[peaks_df['coherence'] > threshold]
|
|||
|
|
out.append(f"Strong peaks (>{threshold:.6f}):")
|
|||
|
|
if len(strong_peaks) == 0:
|
|||
|
|
out.append(" None above threshold — trying mean + 0.5*sigma...")
|
|||
|
|
threshold = global_mean + 0.5 * global_std
|
|||
|
|
strong_peaks = peaks_df[peaks_df['coherence'] > threshold]
|
|||
|
|
|
|||
|
|
for _, row in strong_peaks.iterrows():
|
|||
|
|
out.append(f" Ω={row['omega']:.1f} Coh={row['coherence']:.6f} "
|
|||
|
|
f"K={row['khra_amp']:.3f} G={row['gixx_amp']:.4f}")
|
|||
|
|
|
|||
|
|
# Cluster adjacent peaks
|
|||
|
|
if len(strong_peaks) > 0:
|
|||
|
|
clusters = []
|
|||
|
|
current_cluster = [strong_peaks.iloc[0]['omega']]
|
|||
|
|
for i in range(1, len(strong_peaks)):
|
|||
|
|
if strong_peaks.iloc[i]['omega'] - strong_peaks.iloc[i-1]['omega'] <= 0.15:
|
|||
|
|
current_cluster.append(strong_peaks.iloc[i]['omega'])
|
|||
|
|
else:
|
|||
|
|
clusters.append(current_cluster)
|
|||
|
|
current_cluster = [strong_peaks.iloc[i]['omega']]
|
|||
|
|
clusters.append(current_cluster)
|
|||
|
|
|
|||
|
|
out.append(f"\n {len(clusters)} cluster(s) of strong peaks:")
|
|||
|
|
for i, cl in enumerate(clusters):
|
|||
|
|
center = np.mean(cl)
|
|||
|
|
out.append(f" Cluster {i+1}: Ω ∈ [{min(cl):.1f}, {max(cl):.1f}], center={center:.2f}, width={len(cl)}")
|
|||
|
|
|
|||
|
|
out.append(f"\nFull peak table:")
|
|||
|
|
out.append(f" {'Omega':>6} {'Coherence':>10} {'Khra':>6} {'Gixx':>7} {'Asym':>8} {'Vort':>10}")
|
|||
|
|
for _, row in peaks_df.iterrows():
|
|||
|
|
marker = " *" if row['coherence'] > threshold else " "
|
|||
|
|
out.append(f" {row['omega']:6.1f} {row['coherence']:10.6f} {row['khra_amp']:6.3f} "
|
|||
|
|
f"{row['gixx_amp']:7.4f} {row['asymmetry']:8.4f} {row['vorticity']:10.6f}{marker}")
|
|||
|
|
|
|||
|
|
return peaks_df
|
|||
|
|
|
|||
|
|
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
# Analysis 2: Mode Counting vs Nuclear Shell Degeneracies
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
def analysis_mode_counting(df, out):
|
|||
|
|
out.append("")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
out.append("ANALYSIS 2: Mode Counting vs Nuclear Shell Degeneracies")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
|
|||
|
|
omega_vals = sorted(df['omega'].unique())
|
|||
|
|
|
|||
|
|
# For each omega slice, count distinct coherence levels
|
|||
|
|
# "Distinct" = separated by more than a tolerance
|
|||
|
|
all_coherences = df['coherence'].values
|
|||
|
|
resolution = np.std(all_coherences) * 0.1 # adaptive resolution
|
|||
|
|
if resolution < 1e-6:
|
|||
|
|
resolution = 1e-4
|
|||
|
|
out.append(f"\nCoherence resolution (tolerance): {resolution:.6f}")
|
|||
|
|
|
|||
|
|
mode_counts = {}
|
|||
|
|
for omega in omega_vals:
|
|||
|
|
group = df[df['omega'] == omega]
|
|||
|
|
coh_sorted = np.sort(group['coherence'].values)
|
|||
|
|
# Count distinct levels: merge values within resolution
|
|||
|
|
modes = [coh_sorted[0]]
|
|||
|
|
for c in coh_sorted[1:]:
|
|||
|
|
if c - modes[-1] > resolution:
|
|||
|
|
modes.append(c)
|
|||
|
|
mode_counts[omega] = len(modes)
|
|||
|
|
|
|||
|
|
out.append(f"\nDistinct coherence modes per omega slice:")
|
|||
|
|
out.append(f" {'Omega':>6} {'Modes':>6} {'Nearest Magic':>14} {'Δ':>4}")
|
|||
|
|
total_modes = []
|
|||
|
|
for omega in omega_vals:
|
|||
|
|
n = mode_counts[omega]
|
|||
|
|
total_modes.append(n)
|
|||
|
|
nearest_magic = min(MAGIC_NUMBERS, key=lambda m: abs(m - n))
|
|||
|
|
delta = n - nearest_magic
|
|||
|
|
marker = " <<<" if delta == 0 else ""
|
|||
|
|
out.append(f" {omega:6.1f} {n:6d} {nearest_magic:14d} {delta:+4d}{marker}")
|
|||
|
|
|
|||
|
|
# Overall statistics
|
|||
|
|
mode_arr = np.array(total_modes)
|
|||
|
|
out.append(f"\n Mode count range: [{mode_arr.min()}, {mode_arr.max()}]")
|
|||
|
|
out.append(f" Mean modes: {mode_arr.mean():.1f}")
|
|||
|
|
|
|||
|
|
# Cumulative mode count across all omega
|
|||
|
|
all_coh = np.sort(df['coherence'].unique())
|
|||
|
|
distinct_global = [all_coh[0]]
|
|||
|
|
for c in all_coh[1:]:
|
|||
|
|
if c - distinct_global[-1] > resolution:
|
|||
|
|
distinct_global.append(c)
|
|||
|
|
out.append(f" Total distinct global modes: {len(distinct_global)}")
|
|||
|
|
|
|||
|
|
# Compare to magic numbers
|
|||
|
|
out.append(f"\n Magic number proximity:")
|
|||
|
|
for mn in MAGIC_NUMBERS[:6]:
|
|||
|
|
hits = [omega for omega, n in mode_counts.items() if n == mn]
|
|||
|
|
if hits:
|
|||
|
|
out.append(f" N={mn}: matched at Ω = {', '.join(f'{h:.1f}' for h in hits)}")
|
|||
|
|
else:
|
|||
|
|
closest = min(mode_counts.items(), key=lambda x: abs(x[1] - mn))
|
|||
|
|
out.append(f" N={mn}: no exact match (closest: Ω={closest[0]:.1f} with {closest[1]} modes)")
|
|||
|
|
|
|||
|
|
return mode_counts
|
|||
|
|
|
|||
|
|
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
# Analysis 3: Gap Structure
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
def analysis_gap_structure(df, out):
|
|||
|
|
out.append("")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
out.append("ANALYSIS 3: Gap Structure (Coherence Gaps vs Nuclear Shell Gaps)")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
|
|||
|
|
# Global coherence spectrum: sort all unique values, compute gaps
|
|||
|
|
coh_all = np.sort(df['coherence'].unique())
|
|||
|
|
gaps = np.diff(coh_all)
|
|||
|
|
|
|||
|
|
out.append(f"\nGlobal coherence spectrum: {len(coh_all)} unique values")
|
|||
|
|
out.append(f" Value range: [{coh_all[0]:.6f}, {coh_all[-1]:.6f}]")
|
|||
|
|
out.append(f" Total span: {coh_all[-1] - coh_all[0]:.6f}")
|
|||
|
|
|
|||
|
|
if len(gaps) > 0:
|
|||
|
|
out.append(f"\nGap statistics:")
|
|||
|
|
out.append(f" Mean gap: {gaps.mean():.6f}")
|
|||
|
|
out.append(f" Std gap: {gaps.std():.6f}")
|
|||
|
|
out.append(f" Min gap: {gaps.min():.6f}")
|
|||
|
|
out.append(f" Max gap: {gaps.max():.6f}")
|
|||
|
|
|
|||
|
|
# Find the largest gaps — these correspond to "shell closures"
|
|||
|
|
n_top = min(10, len(gaps))
|
|||
|
|
top_idx = np.argsort(gaps)[-n_top:][::-1]
|
|||
|
|
out.append(f"\n Top {n_top} largest gaps (shell boundaries):")
|
|||
|
|
out.append(f" {'Rank':>4} {'Gap':>10} {'Below':>10} {'Above':>10} {'Ratio':>8}")
|
|||
|
|
gap_ratios = []
|
|||
|
|
for rank, idx in enumerate(top_idx):
|
|||
|
|
ratio = gaps[idx] / gaps.mean() if gaps.mean() > 0 else 0
|
|||
|
|
gap_ratios.append(gaps[idx])
|
|||
|
|
out.append(f" {rank+1:4d} {gaps[idx]:10.6f} {coh_all[idx]:10.6f} "
|
|||
|
|
f"{coh_all[idx+1]:10.6f} {ratio:8.2f}x")
|
|||
|
|
|
|||
|
|
# Compare gap ratios to nuclear shell gap ratios
|
|||
|
|
if len(gap_ratios) >= 5:
|
|||
|
|
observed_ratios = np.array(gap_ratios[:5]) / gap_ratios[0]
|
|||
|
|
out.append(f"\n Gap ratio comparison (top 5 gaps, normalized to largest):")
|
|||
|
|
out.append(f" Observed: {', '.join(f'{r:.3f}' for r in observed_ratios)}")
|
|||
|
|
out.append(f" Nuclear: {', '.join(f'{r:.3f}' for r in MAGIC_GAP_RATIOS)}")
|
|||
|
|
correlation = np.corrcoef(observed_ratios, MAGIC_GAP_RATIOS[:5])[0, 1]
|
|||
|
|
out.append(f" Pearson correlation: {correlation:.4f}")
|
|||
|
|
|
|||
|
|
# Per-omega gap structure
|
|||
|
|
out.append(f"\n Per-omega max gap:")
|
|||
|
|
omega_vals = sorted(df['omega'].unique())
|
|||
|
|
for omega in omega_vals:
|
|||
|
|
group = df[df['omega'] == omega]
|
|||
|
|
coh_sorted = np.sort(group['coherence'].values)
|
|||
|
|
g = np.diff(coh_sorted)
|
|||
|
|
if len(g) > 0:
|
|||
|
|
max_gap = g.max()
|
|||
|
|
mean_gap = g.mean()
|
|||
|
|
ratio = max_gap / mean_gap if mean_gap > 0 else 0
|
|||
|
|
out.append(f" Ω={omega:.1f}: max_gap={max_gap:.6f} mean_gap={mean_gap:.6f} "
|
|||
|
|
f"ratio={ratio:.2f}x")
|
|||
|
|
|
|||
|
|
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
# Analysis 4: Omega-Resolved Shell Occupancy
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
def analysis_shell_occupancy(df, out):
|
|||
|
|
out.append("")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
out.append("ANALYSIS 4: Omega-Resolved Shell Occupancy")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
|
|||
|
|
global_mean = df['coherence'].mean()
|
|||
|
|
global_std = df['coherence'].std()
|
|||
|
|
|
|||
|
|
# Define "shells" as coherence bands
|
|||
|
|
n_shells = 6
|
|||
|
|
coh_min = df['coherence'].min()
|
|||
|
|
coh_max = df['coherence'].max()
|
|||
|
|
shell_edges = np.linspace(coh_min, coh_max + 1e-9, n_shells + 1)
|
|||
|
|
|
|||
|
|
out.append(f"\nShell definition: {n_shells} equal-width coherence bands")
|
|||
|
|
out.append(f" Coherence range: [{coh_min:.6f}, {coh_max:.6f}]")
|
|||
|
|
out.append(f" Shell width: {(coh_max - coh_min) / n_shells:.6f}")
|
|||
|
|
out.append("")
|
|||
|
|
|
|||
|
|
omega_vals = sorted(df['omega'].unique())
|
|||
|
|
|
|||
|
|
# Build occupancy matrix: omega × shell
|
|||
|
|
occupancy = np.zeros((len(omega_vals), n_shells), dtype=int)
|
|||
|
|
for i, omega in enumerate(omega_vals):
|
|||
|
|
group = df[df['omega'] == omega]
|
|||
|
|
for j in range(n_shells):
|
|||
|
|
count = ((group['coherence'] >= shell_edges[j]) &
|
|||
|
|
(group['coherence'] < shell_edges[j+1])).sum()
|
|||
|
|
occupancy[i, j] = count
|
|||
|
|
|
|||
|
|
# Display occupancy matrix
|
|||
|
|
header = f" {'Omega':>6} " + " ".join(f"S{j+1:d}" for j in range(n_shells)) + " Total Pattern"
|
|||
|
|
out.append(header)
|
|||
|
|
for i, omega in enumerate(omega_vals):
|
|||
|
|
row = occupancy[i]
|
|||
|
|
total = row.sum()
|
|||
|
|
# Binary pattern: 1 if occupied, 0 if not
|
|||
|
|
pattern = "".join("█" if x > 0 else "·" for x in row)
|
|||
|
|
out.append(f" {omega:6.1f} " + " ".join(f"{x:2d}" for x in row) +
|
|||
|
|
f" {total:5d} {pattern}")
|
|||
|
|
|
|||
|
|
# Count unique occupancy patterns
|
|||
|
|
patterns = ["".join("1" if x > 0 else "0" for x in occupancy[i]) for i in range(len(omega_vals))]
|
|||
|
|
pattern_counts = Counter(patterns)
|
|||
|
|
out.append(f"\n Unique occupancy patterns: {len(pattern_counts)}")
|
|||
|
|
for pat, count in sorted(pattern_counts.items(), key=lambda x: -x[1]):
|
|||
|
|
visual = "".join("█" if c == "1" else "·" for c in pat)
|
|||
|
|
out.append(f" {visual} ({pat}): {count} omega values")
|
|||
|
|
|
|||
|
|
# Shell filling: total occupancy per shell across all omega
|
|||
|
|
shell_totals = occupancy.sum(axis=0)
|
|||
|
|
out.append(f"\n Total occupancy per shell:")
|
|||
|
|
for j in range(n_shells):
|
|||
|
|
bar = "█" * (shell_totals[j] // 2) if shell_totals[j] > 0 else ""
|
|||
|
|
out.append(f" S{j+1} [{shell_edges[j]:.5f} – {shell_edges[j+1]:.5f}]: "
|
|||
|
|
f"{shell_totals[j]:4d} {bar}")
|
|||
|
|
|
|||
|
|
# Compare to nuclear filling order
|
|||
|
|
if len(omega_vals) >= 5:
|
|||
|
|
# "Closed shell" = omega where all points fall in same shell
|
|||
|
|
closed = []
|
|||
|
|
for i, omega in enumerate(omega_vals):
|
|||
|
|
nonzero = np.count_nonzero(occupancy[i])
|
|||
|
|
if nonzero == 1:
|
|||
|
|
filled_shell = np.argmax(occupancy[i])
|
|||
|
|
closed.append((omega, filled_shell + 1))
|
|||
|
|
out.append(f"\n Closed-shell configurations (all points in one band):")
|
|||
|
|
if closed:
|
|||
|
|
for omega, shell in closed:
|
|||
|
|
out.append(f" Ω={omega:.1f} → Shell {shell}")
|
|||
|
|
else:
|
|||
|
|
out.append(f" None found (points spread across multiple bands)")
|
|||
|
|
|
|||
|
|
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
# Analysis 5: 2D Torus Mode Comparison
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
def analysis_torus_modes(df, out):
|
|||
|
|
out.append("")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
out.append("ANALYSIS 5: 2D Torus Mode Comparison")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
|
|||
|
|
# On a 2D torus (periodic lattice), modes are labeled (n, m)
|
|||
|
|
# with energy ~ n² + m². Degeneracy = # of (n,m) pairs giving same E.
|
|||
|
|
# This is the sum-of-two-squares function r₂(E).
|
|||
|
|
max_E = 50
|
|||
|
|
torus_degeneracy = {}
|
|||
|
|
for n in range(-int(np.sqrt(max_E)) - 1, int(np.sqrt(max_E)) + 2):
|
|||
|
|
for m in range(-int(np.sqrt(max_E)) - 1, int(np.sqrt(max_E)) + 2):
|
|||
|
|
E = n * n + m * m
|
|||
|
|
if 0 < E <= max_E:
|
|||
|
|
torus_degeneracy[E] = torus_degeneracy.get(E, 0) + 1
|
|||
|
|
|
|||
|
|
torus_energies = sorted(torus_degeneracy.keys())
|
|||
|
|
torus_degens = [torus_degeneracy[E] for E in torus_energies]
|
|||
|
|
|
|||
|
|
out.append(f"\nTheoretical 2D torus modes (E = n² + m², E ≤ {max_E}):")
|
|||
|
|
out.append(f" Representable energies: {len(torus_energies)}")
|
|||
|
|
out.append(f" Cumulative modes at each energy:")
|
|||
|
|
cumulative = np.cumsum(torus_degens)
|
|||
|
|
out.append(f" {'E':>4} {'Degen':>6} {'Cumul':>6} {'Magic?':>7}")
|
|||
|
|
for E, d, c in zip(torus_energies, torus_degens, cumulative):
|
|||
|
|
magic_hit = " <<<" if c in MAGIC_NUMBERS else ""
|
|||
|
|
out.append(f" {E:4d} {d:6d} {c:6d}{magic_hit}")
|
|||
|
|
|
|||
|
|
# Compare torus cumulative degeneracies to magic numbers
|
|||
|
|
magic_hits = []
|
|||
|
|
for mn in MAGIC_NUMBERS[:6]:
|
|||
|
|
if mn in cumulative.tolist():
|
|||
|
|
idx = cumulative.tolist().index(mn)
|
|||
|
|
magic_hits.append((mn, torus_energies[idx]))
|
|||
|
|
out.append(f"\n Torus shell closures matching magic numbers:")
|
|||
|
|
if magic_hits:
|
|||
|
|
for mn, E in magic_hits:
|
|||
|
|
out.append(f" Magic N={mn} occurs at torus energy E={E}")
|
|||
|
|
else:
|
|||
|
|
out.append(f" No exact matches")
|
|||
|
|
# Find nearest
|
|||
|
|
for mn in MAGIC_NUMBERS[:6]:
|
|||
|
|
nearest_idx = np.argmin(np.abs(cumulative - mn))
|
|||
|
|
out.append(f" Magic N={mn}: nearest cumulative = {cumulative[nearest_idx]} at E={torus_energies[nearest_idx]}")
|
|||
|
|
|
|||
|
|
# Now compare to observed data
|
|||
|
|
# Use coherence as proxy for "energy level"
|
|||
|
|
# Count modes in the observed spectrum and compare degeneracies
|
|||
|
|
omega_vals = sorted(df['omega'].unique())
|
|||
|
|
resolution = np.std(df['coherence'].values) * 0.1
|
|||
|
|
if resolution < 1e-6:
|
|||
|
|
resolution = 1e-4
|
|||
|
|
|
|||
|
|
# Per-omega mode degeneracy: count how many (khra, gixx) pairs
|
|||
|
|
# give the same coherence level (within resolution)
|
|||
|
|
out.append(f"\n Observed mode degeneracies per omega:")
|
|||
|
|
out.append(f" {'Omega':>6} {'Modes':>6} {'Max Degen':>10} {'Degen Pattern':>20}")
|
|||
|
|
for omega in omega_vals:
|
|||
|
|
group = df[df['omega'] == omega]
|
|||
|
|
coh_sorted = np.sort(group['coherence'].values)
|
|||
|
|
# Bin into distinct modes
|
|||
|
|
modes = []
|
|||
|
|
current_mode = [coh_sorted[0]]
|
|||
|
|
for c in coh_sorted[1:]:
|
|||
|
|
if c - current_mode[-1] > resolution:
|
|||
|
|
modes.append(len(current_mode))
|
|||
|
|
current_mode = [c]
|
|||
|
|
else:
|
|||
|
|
current_mode.append(c)
|
|||
|
|
modes.append(len(current_mode))
|
|||
|
|
# modes[] now holds the degeneracy of each mode
|
|||
|
|
max_degen = max(modes)
|
|||
|
|
pattern = ",".join(str(d) for d in modes[:8])
|
|||
|
|
if len(modes) > 8:
|
|||
|
|
pattern += "..."
|
|||
|
|
out.append(f" {omega:6.1f} {len(modes):6d} {max_degen:10d} {pattern:>20}")
|
|||
|
|
|
|||
|
|
# Correlation between observed degeneracy spectrum and torus degeneracies
|
|||
|
|
# Use the full dataset: histogram of degeneracies
|
|||
|
|
all_coh = np.sort(df['coherence'].values)
|
|||
|
|
global_modes = []
|
|||
|
|
current_mode = [all_coh[0]]
|
|||
|
|
for c in all_coh[1:]:
|
|||
|
|
if c - current_mode[-1] > resolution:
|
|||
|
|
global_modes.append(len(current_mode))
|
|||
|
|
current_mode = [c]
|
|||
|
|
else:
|
|||
|
|
current_mode.append(c)
|
|||
|
|
global_modes.append(len(current_mode))
|
|||
|
|
|
|||
|
|
obs_degen_hist = Counter(global_modes)
|
|||
|
|
torus_degen_hist = Counter(torus_degens)
|
|||
|
|
|
|||
|
|
out.append(f"\n Degeneracy histograms:")
|
|||
|
|
out.append(f" {'Degen':>6} {'Observed':>9} {'Torus':>6}")
|
|||
|
|
all_degens = sorted(set(list(obs_degen_hist.keys()) + list(torus_degen_hist.keys())))
|
|||
|
|
for d in all_degens[:15]:
|
|||
|
|
out.append(f" {d:6d} {obs_degen_hist.get(d, 0):9d} {torus_degen_hist.get(d, 0):6d}")
|
|||
|
|
|
|||
|
|
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
# Analysis 6: GUE Pair Correlation (Random Matrix Theory)
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
def analysis_gue_correlation(df, out):
|
|||
|
|
out.append("")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
out.append("ANALYSIS 6: GUE Pair Correlation (Random Matrix Theory)")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
|
|||
|
|
# Nearest-neighbor spacing distribution
|
|||
|
|
# For GUE (β=2): P(s) = (32/π²) s² exp(-4s²/π) (Wigner surmise)
|
|||
|
|
# For Poisson: P(s) = exp(-s)
|
|||
|
|
# Normalize spacings to mean = 1
|
|||
|
|
|
|||
|
|
coh_sorted = np.sort(df['coherence'].unique())
|
|||
|
|
spacings = np.diff(coh_sorted)
|
|||
|
|
|
|||
|
|
if len(spacings) < 5:
|
|||
|
|
out.append("\n Insufficient unique coherence values for spacing analysis.")
|
|||
|
|
return
|
|||
|
|
|
|||
|
|
mean_spacing = spacings.mean()
|
|||
|
|
if mean_spacing > 0:
|
|||
|
|
s_normalized = spacings / mean_spacing # normalize to <s> = 1
|
|||
|
|
else:
|
|||
|
|
out.append("\n Zero mean spacing — all values identical.")
|
|||
|
|
return
|
|||
|
|
|
|||
|
|
out.append(f"\nSpacing statistics (normalized to mean=1):")
|
|||
|
|
out.append(f" N unique levels: {len(coh_sorted)}")
|
|||
|
|
out.append(f" N spacings: {len(spacings)}")
|
|||
|
|
out.append(f" Raw mean spacing: {mean_spacing:.6e}")
|
|||
|
|
out.append(f" Normalized <s>: {s_normalized.mean():.4f}")
|
|||
|
|
out.append(f" Normalized var: {np.var(s_normalized):.4f}")
|
|||
|
|
out.append(f" Normalized <s²>: {np.mean(s_normalized**2):.4f}")
|
|||
|
|
|
|||
|
|
# GUE prediction: var(s) = (4 - π) * π / (2π²) ≈ 0.178
|
|||
|
|
# Poisson prediction: var(s) = 1.0
|
|||
|
|
gue_var = (4 - np.pi) * np.pi / (2 * np.pi**2)
|
|||
|
|
obs_var = np.var(s_normalized)
|
|||
|
|
out.append(f"\n Variance comparison:")
|
|||
|
|
out.append(f" Observed: {obs_var:.4f}")
|
|||
|
|
out.append(f" GUE (β=2): {gue_var:.4f}")
|
|||
|
|
out.append(f" Poisson: 1.0000")
|
|||
|
|
|
|||
|
|
if abs(obs_var - gue_var) < abs(obs_var - 1.0):
|
|||
|
|
out.append(f" → Closer to GUE (level repulsion present)")
|
|||
|
|
else:
|
|||
|
|
out.append(f" → Closer to Poisson (uncorrelated levels)")
|
|||
|
|
|
|||
|
|
# Histogram of normalized spacings
|
|||
|
|
n_bins = 20
|
|||
|
|
bin_edges = np.linspace(0, max(4.0, s_normalized.max()), n_bins + 1)
|
|||
|
|
hist, _ = np.histogram(s_normalized, bins=bin_edges, density=True)
|
|||
|
|
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
|
|||
|
|
|
|||
|
|
# Theoretical curves
|
|||
|
|
gue_pdf = (32.0 / (np.pi**2)) * bin_centers**2 * np.exp(-4.0 * bin_centers**2 / np.pi)
|
|||
|
|
poisson_pdf = np.exp(-bin_centers)
|
|||
|
|
|
|||
|
|
out.append(f"\n Spacing distribution P(s):")
|
|||
|
|
out.append(f" {'s':>6} {'Observed':>9} {'GUE':>7} {'Poisson':>8}")
|
|||
|
|
for i in range(n_bins):
|
|||
|
|
out.append(f" {bin_centers[i]:6.2f} {hist[i]:9.4f} {gue_pdf[i]:7.4f} {poisson_pdf[i]:8.4f}")
|
|||
|
|
|
|||
|
|
# Chi-squared goodness of fit (manual, no scipy)
|
|||
|
|
# Against GUE and Poisson
|
|||
|
|
chi2_gue = 0
|
|||
|
|
chi2_poisson = 0
|
|||
|
|
bins_used = 0
|
|||
|
|
for i in range(n_bins):
|
|||
|
|
if gue_pdf[i] > 0.01: # only use bins with sufficient expected density
|
|||
|
|
chi2_gue += (hist[i] - gue_pdf[i])**2 / gue_pdf[i]
|
|||
|
|
bins_used += 1
|
|||
|
|
if poisson_pdf[i] > 0.01:
|
|||
|
|
chi2_poisson += (hist[i] - poisson_pdf[i])**2 / poisson_pdf[i]
|
|||
|
|
|
|||
|
|
out.append(f"\n Goodness of fit (χ²-like, lower is better):")
|
|||
|
|
out.append(f" vs GUE: {chi2_gue:.4f} (over {bins_used} bins)")
|
|||
|
|
out.append(f" vs Poisson: {chi2_poisson:.4f}")
|
|||
|
|
if chi2_gue < chi2_poisson:
|
|||
|
|
out.append(f" → GUE is better fit")
|
|||
|
|
else:
|
|||
|
|
out.append(f" → Poisson is better fit")
|
|||
|
|
|
|||
|
|
# Number variance Σ²(L): count fluctuations in intervals of length L
|
|||
|
|
out.append(f"\n Number variance Σ²(L):")
|
|||
|
|
out.append(f" {'L':>6} {'Σ²(obs)':>9} {'GUE':>7} {'Poisson':>8}")
|
|||
|
|
for L in [0.5, 1.0, 1.5, 2.0, 3.0, 5.0]:
|
|||
|
|
# Count how many spacings fall in windows of size L*mean_spacing
|
|||
|
|
window = L * mean_spacing
|
|||
|
|
counts = []
|
|||
|
|
for start_idx in range(len(coh_sorted) - 1):
|
|||
|
|
start_val = coh_sorted[start_idx]
|
|||
|
|
# Count levels in [start_val, start_val + window)
|
|||
|
|
n_in_window = np.sum((coh_sorted >= start_val) & (coh_sorted < start_val + window))
|
|||
|
|
counts.append(n_in_window)
|
|||
|
|
counts = np.array(counts, dtype=float)
|
|||
|
|
sigma2_obs = np.var(counts) if len(counts) > 0 else 0
|
|||
|
|
|
|||
|
|
# GUE: Σ²(L) ≈ (2/π²)(ln(2πL) + γ + 1) for large L (γ = Euler-Mascheroni)
|
|||
|
|
gamma_em = 0.5772156649
|
|||
|
|
sigma2_gue = (2.0 / np.pi**2) * (np.log(2 * np.pi * L) + gamma_em + 1) if L > 0 else 0
|
|||
|
|
sigma2_poisson = L # Poisson: Σ²(L) = L
|
|||
|
|
|
|||
|
|
out.append(f" {L:6.1f} {sigma2_obs:9.4f} {sigma2_gue:7.4f} {sigma2_poisson:8.4f}")
|
|||
|
|
|
|||
|
|
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
# Main
|
|||
|
|
# ═══════════════════════════════════════════════════════════════════
|
|||
|
|
def main():
|
|||
|
|
if len(sys.argv) < 2:
|
|||
|
|
# Auto-find latest sweep CSV
|
|||
|
|
sweep_dir = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
|
|||
|
|
"sweep_results")
|
|||
|
|
csvs = sorted([f for f in os.listdir(sweep_dir) if f.startswith("em_direct_sweep") and f.endswith(".csv")])
|
|||
|
|
if not csvs:
|
|||
|
|
print("Usage: python3 nuclear_magic_analyzer.py <sweep_csv>")
|
|||
|
|
print(" No sweep CSVs found in sweep_results/")
|
|||
|
|
sys.exit(1)
|
|||
|
|
csv_path = os.path.join(sweep_dir, csvs[-1])
|
|||
|
|
print(f"Auto-selected latest sweep: {csvs[-1]}")
|
|||
|
|
else:
|
|||
|
|
csv_path = sys.argv[1]
|
|||
|
|
|
|||
|
|
if not os.path.exists(csv_path):
|
|||
|
|
print(f"ERROR: File not found: {csv_path}")
|
|||
|
|
sys.exit(1)
|
|||
|
|
|
|||
|
|
df = load_sweep(csv_path)
|
|||
|
|
print(f"Loaded {len(df)} data points from {os.path.basename(csv_path)}")
|
|||
|
|
print(f" Omega range: {df['omega'].min():.1f} – {df['omega'].max():.1f}")
|
|||
|
|
print(f" Coherence range: {df['coherence'].min():.6f} – {df['coherence'].max():.6f}")
|
|||
|
|
print()
|
|||
|
|
|
|||
|
|
out = []
|
|||
|
|
out.append("╔══════════════════════════════════════════════════════════════════════╗")
|
|||
|
|
out.append("║ NUCLEAR MAGIC NUMBER ANALYSIS — RESONANCE ENGINE ║")
|
|||
|
|
out.append("╚══════════════════════════════════════════════════════════════════════╝")
|
|||
|
|
out.append(f"Source: {os.path.basename(csv_path)}")
|
|||
|
|
out.append(f"Points: {len(df)}")
|
|||
|
|
out.append(f"Date: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
|
|||
|
|
out.append(f"Omega: {df['omega'].min():.1f} – {df['omega'].max():.1f} ({df['omega'].nunique()} steps)")
|
|||
|
|
out.append(f"Coherence: {df['coherence'].min():.6f} – {df['coherence'].max():.6f}")
|
|||
|
|
|
|||
|
|
# Run all six analyses
|
|||
|
|
peaks_df = analysis_coherence_peaks(df, out)
|
|||
|
|
mode_counts = analysis_mode_counting(df, out)
|
|||
|
|
analysis_gap_structure(df, out)
|
|||
|
|
analysis_shell_occupancy(df, out)
|
|||
|
|
analysis_torus_modes(df, out)
|
|||
|
|
analysis_gue_correlation(df, out)
|
|||
|
|
|
|||
|
|
# Summary
|
|||
|
|
out.append("")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
out.append("SUMMARY")
|
|||
|
|
out.append("=" * 70)
|
|||
|
|
|
|||
|
|
best_omega = peaks_df.loc[peaks_df['coherence'].idxmax()]
|
|||
|
|
out.append(f"\n Best coherence: {best_omega['coherence']:.6f} at "
|
|||
|
|
f"Ω={best_omega['omega']:.1f} K={best_omega['khra_amp']:.3f} G={best_omega['gixx_amp']:.4f}")
|
|||
|
|
|
|||
|
|
mode_arr = np.array(list(mode_counts.values()))
|
|||
|
|
out.append(f" Mode count range: {mode_arr.min()} – {mode_arr.max()}")
|
|||
|
|
|
|||
|
|
magic_matches = sum(1 for n in mode_counts.values() if n in MAGIC_NUMBERS)
|
|||
|
|
out.append(f" Omega slices matching a magic number: {magic_matches}/{len(mode_counts)}")
|
|||
|
|
|
|||
|
|
out.append(f"\n Nuclear magic numbers for reference: {MAGIC_NUMBERS}")
|
|||
|
|
out.append("")
|
|||
|
|
|
|||
|
|
# Print to stdout
|
|||
|
|
report = "\n".join(out)
|
|||
|
|
print(report)
|
|||
|
|
|
|||
|
|
# Save to file
|
|||
|
|
results_dir = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
|
|||
|
|
"results")
|
|||
|
|
os.makedirs(results_dir, exist_ok=True)
|
|||
|
|
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
|
|||
|
|
output_path = os.path.join(results_dir, f"nuclear_magic_analysis_{timestamp}.txt")
|
|||
|
|
with open(output_path, 'w', encoding='utf-8') as f:
|
|||
|
|
f.write(report)
|
|||
|
|
print(f"\nSaved to: {output_path}")
|
|||
|
|
|
|||
|
|
|
|||
|
|
if __name__ == "__main__":
|
|||
|
|
main()
|