"""
LFM Double-Slit: RESONANT DETECTOR
==================================

The detector is tuned to RESONATE with the incoming wave.
This should maximize energy transfer = actual detection.

Key: detector natural frequency ≈ wave frequency

Author: Greg Partin
Date: February 4, 2026
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy.signal import find_peaks

print("="*70)
print("LFM DOUBLE-SLIT: RESONANT DETECTOR")
print("Detector tuned to incoming wave frequency")
print("="*70)

# Grid
N = 300
L = 30.0
dx = L / N
x = np.linspace(-L/2, L/2, N)
y = np.linspace(-L/2, L/2, N)
X, Y = np.meshgrid(x, y, indexing='ij')

c = 1.0
dt = 0.02

chi_0 = 1.0
kappa = 1.0  # Strong E-chi coupling
E0_sq = 0.0

# Barrier
barrier_x = 0.0
barrier_width = 0.3
barrier_chi = 10.0
slit_separation = 4.0
slit_width = 1.5

barrier_region = (np.abs(X - barrier_x) < barrier_width/2)
slit_1 = (np.abs(Y - slit_separation/2) < slit_width/2)
slit_2 = (np.abs(Y + slit_separation/2) < slit_width/2)
barrier_mask = barrier_region & ~slit_1 & ~slit_2

# Source
source_x = -10.0
source_freq = 3.0
R_source = np.sqrt((X - source_x)**2 + Y**2)
source_envelope = np.exp(-R_source**2 / 1.0)

# Screen
screen_x = 10.0
screen_idx = int((screen_x + L/2) / dx)

# RESONANT DETECTOR
# Natural frequency of E in chi well: omega = chi (in our units)
# For resonance: chi_detector ≈ source_freq
detector_x = barrier_x
detector_y = -slit_separation / 2
detector_radius = slit_width * 0.4  # Smaller, concentrated
detector_chi = source_freq * 0.8  # Tuned near source frequency!

R_detector = np.sqrt((X - detector_x)**2 + (Y - detector_y)**2)
detector_mask = R_detector < detector_radius
detector_shape = np.exp(-R_detector**2 / (detector_radius**2))

print(f"Source frequency: {source_freq}")
print(f"Detector chi (= natural freq): {detector_chi}")
print(f"Detuning: {abs(source_freq - detector_chi)/source_freq*100:.0f}%")

# Helpers
def laplacian(f, dx):
    lap = np.zeros_like(f)
    lap[1:-1, 1:-1] = (f[2:, 1:-1] + f[:-2, 1:-1] + 
                       f[1:-1, 2:] + f[1:-1, :-2] - 
                       4*f[1:-1, 1:-1]) / dx**2
    return lap

def absorb_boundary(f, w=15, s=0.05):
    for i in range(w):
        d = s * ((w-i)/w)**2
        f[i,:] *= (1-d); f[-(i+1),:] *= (1-d)
        f[:,i] *= (1-d); f[:,-(i+1)] *= (1-d)
    return f

# Simulation
def run(use_detector, n_steps, name):
    print(f"\n--- {name} ---")
    
    E = np.zeros((N, N))
    E_prev = E.copy()
    
    chi = np.ones((N, N)) * chi_0
    chi[barrier_mask] = barrier_chi
    
    if use_detector:
        # Set chi IN the detector to the resonant value
        chi[detector_mask] = detector_chi
        
        # Small initial E disturbance (detector ready to resonate)
        E += 0.1 * detector_shape
        E_prev = E.copy()
        
        print(f"  Detector chi = {detector_chi:.2f} (resonant with wave)")
    
    chi_prev = chi.copy()
    chi_fixed = chi.copy()  # Save fixed chi structure
    
    frames = []
    screen_history = []
    det_E = []
    det_E_signed = []  # Track actual E, not E²
    
    for step in range(n_steps):
        t = step * dt
        
        drive = 1.0 * np.sin(source_freq * t)
        
        # GOV-01
        lap_E = laplacian(E, dx)
        E_new = 2*E - E_prev + dt**2 * (c**2 * lap_E - chi**2 * E)
        E_new += dt**2 * drive * source_envelope * 10.0
        
        # GOV-02 with coupling (but maintain detector structure)
        lap_chi = laplacian(chi, dx)
        chi_new = 2*chi - chi_prev + dt**2 * (c**2 * lap_chi - kappa * (E**2 - E0_sq))
        chi_new = np.clip(chi_new, 0.05, barrier_chi + 1)
        
        # Keep barrier and detector structure fixed
        chi_new[barrier_mask] = barrier_chi
        if use_detector:
            chi_new[detector_mask] = detector_chi
        
        E_new = absorb_boundary(E_new)
        
        E_prev, E = E, E_new
        chi_prev, chi = chi, chi_new
        
        screen_history.append(E[screen_idx, :].copy())
        
        if use_detector:
            det_E.append(np.sum(E[detector_mask]**2))
            det_E_signed.append(np.mean(E[detector_mask]))
        
        if step % 50 == 0:
            frames.append({'E': E.copy(), 'chi': chi.copy(), 't': t})
        
        if step % 500 == 0:
            det_info = f", det_E={det_E[-1]:.2f}" if det_E else ""
            print(f"  Step {step}: E_max={np.max(np.abs(E)):.2f}{det_info}")
    
    pattern = np.sum(np.array(screen_history)**2, axis=0)
    return frames, pattern, det_E, det_E_signed

n_steps = 3000

frames_no, pattern_no, _, _ = run(False, n_steps, "NO DETECTOR")
frames_det, pattern_det, det_E, det_E_signed = run(True, n_steps, "RESONANT DETECTOR")

# Analysis
print("\n" + "="*70)
print("ANALYSIS")
print("="*70)

pmax = max(np.max(pattern_no), np.max(pattern_det))
p_no = pattern_no / pmax
p_det = pattern_det / pmax

peaks_no, _ = find_peaks(p_no, height=0.02, distance=5, prominence=0.01)
peaks_det, _ = find_peaks(p_det, height=0.02, distance=5, prominence=0.01)

def vis(p, peaks):
    if len(peaks) < 2: return 0
    pk = p[peaks]
    valleys = [np.min(p[peaks[i]:peaks[i+1]]) for i in range(len(peaks)-1)]
    if not valleys: return 0
    Imax, Imin = np.mean(pk), np.mean(valleys)
    return (Imax - Imin)/(Imax + Imin) if (Imax + Imin) > 0 else 0

v_no = vis(p_no, peaks_no)
v_det = vis(p_det, peaks_det)

print(f"\nWithout detector: {len(peaks_no)} peaks, V = {v_no:.3f}")
print(f"With detector:    {len(peaks_det)} peaks, V = {v_det:.3f}")

if det_E:
    init = np.mean(det_E[:10])
    peak = max(det_E)
    print(f"\nDetector response:")
    print(f"  Initial E²: {init:.2f}")
    print(f"  Peak E²:    {peak:.2f}")
    print(f"  Excitation: {peak/max(init,0.01):.1f}x")
    
    # Check for resonant oscillation
    if len(det_E_signed) > 100:
        oscillation_amp = np.std(det_E_signed[len(det_E_signed)//2:])
        print(f"  Oscillation amplitude: {oscillation_amp:.3f}")

# Pattern asymmetry (key test!)
# If detector absorbed from bottom slit, pattern should shift
mid = N // 2
top_intensity = np.sum(pattern_det[mid:])
bot_intensity = np.sum(pattern_det[:mid])
asymmetry = (top_intensity - bot_intensity) / (top_intensity + bot_intensity)

top_intensity_no = np.sum(pattern_no[mid:])
bot_intensity_no = np.sum(pattern_no[:mid])
asymmetry_no = (top_intensity_no - bot_intensity_no) / (top_intensity_no + bot_intensity_no)

print(f"\nPattern asymmetry (top vs bottom of screen):")
print(f"  Without detector: {asymmetry_no:.3f}")
print(f"  With detector:    {asymmetry:.3f}")

# Plot
fig, axes = plt.subplots(2, 3, figsize=(14, 9))

ax = axes[0, 0]
frame = frames_no[len(frames_no)*2//3]
ax.imshow(frame['E'].T, origin='lower', extent=[-L/2,L/2,-L/2,L/2],
          cmap='RdBu', vmin=-0.3, vmax=0.3)
ax.axvline(barrier_x, color='k', lw=1)
ax.axvline(screen_x, color='g', ls='--')
ax.set_title('NO Detector')

ax = axes[0, 1]
frame = frames_det[len(frames_det)*2//3]
ax.imshow(frame['E'].T, origin='lower', extent=[-L/2,L/2,-L/2,L/2],
          cmap='RdBu', vmin=-0.3, vmax=0.3)
ax.axvline(barrier_x, color='k', lw=1)
ax.axvline(screen_x, color='g', ls='--')
c = plt.Circle((detector_x, detector_y), detector_radius, fill=False, color='yellow', lw=2)
ax.add_patch(c)
ax.set_title('WITH Resonant Detector')

ax = axes[0, 2]
ax.imshow(frame['chi'].T, origin='lower', extent=[-L/2,L/2,-L/2,L/2],
          cmap='viridis', vmin=0, vmax=barrier_chi)
c = plt.Circle((detector_x, detector_y), detector_radius, fill=False, color='white', lw=2)
ax.add_patch(c)
ax.set_title(f'Chi field (detector chi={detector_chi:.1f})')

ax = axes[1, 0]
ax.plot(y, p_no, 'b-', lw=2, label=f'No det (V={v_no:.2f})')
ax.plot(y, p_det, 'r--', lw=2, label=f'With det (V={v_det:.2f})')
ax.legend()
ax.set_title('Screen Patterns')
ax.set_xlim(-10, 10)
ax.grid(alpha=0.3)

ax = axes[1, 1]
# Plot difference
diff = p_det - p_no
ax.plot(y, diff, 'purple', lw=2)
ax.axhline(0, color='k', ls='--')
ax.set_title('Difference (with - without detector)')
ax.set_xlim(-10, 10)
ax.grid(alpha=0.3)
ax.set_ylabel('Intensity difference')

ax = axes[1, 2]
if det_E_signed:
    t_arr = np.arange(len(det_E_signed))*dt
    ax.plot(t_arr, det_E_signed, 'g-', lw=0.5)
    ax.set_title('Detector E (showing oscillation)')
    ax.set_xlabel('Time')
    ax.set_ylabel('E at detector')
    ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('figures/resonant_detector_result.png', dpi=150)
plt.show()

print(f"\nSaved: figures/resonant_detector_result.png")

# GIFs
print("\nCreating GIFs...")

def make_gif(frames, fname, title, show_det=False):
    fig, ax = plt.subplots(figsize=(7,7))
    def update(i):
        ax.clear()
        ax.imshow(frames[i]['E'].T, origin='lower', 
                  extent=[-L/2,L/2,-L/2,L/2], cmap='RdBu', vmin=-0.3, vmax=0.3)
        ax.axvline(barrier_x, color='k', lw=1)
        ax.axvline(screen_x, color='g', ls='--')
        if show_det:
            c = plt.Circle((detector_x, detector_y), detector_radius, 
                          fill=False, color='yellow', lw=2)
            ax.add_patch(c)
        ax.set_title(f"{title}\nt={frames[i]['t']:.1f}")
        return [ax]
    anim = FuncAnimation(fig, update, frames=len(frames), interval=80)
    anim.save(fname, writer=PillowWriter(fps=12))
    plt.close()
    print(f"  {fname}")

make_gif(frames_no, 'figures/double_slit_no_detector.gif', 'NO Detector', False)
make_gif(frames_det, 'figures/double_slit_with_detector.gif', 'WITH Resonant Detector', True)

# Conclusion
print("\n" + "="*70)
print("CONCLUSION")
print("="*70)

print(f"""
RESONANT DETECTOR EXPERIMENT (pure GOV-01 + GOV-02):

Interference visibility:
  Without detector: {v_no:.3f}
  With detector:    {v_det:.3f}
  Change: {100*(v_det-v_no)/max(v_no,0.01):.1f}%

Pattern asymmetry:
  Without: {asymmetry_no:.3f}
  With:    {asymmetry:.3f}
""")

excited = max(det_E) > np.mean(det_E[:10]) * 5 if det_E else False
vis_reduced = v_det < v_no * 0.7
asymmetric = abs(asymmetry - asymmetry_no) > 0.05

if vis_reduced:
    print("✓ Visibility REDUCED - detector affected coherence")
else:
    print("✗ Visibility PRESERVED - detector did not break coherence")

if excited:
    print("✓ Detector WAS excited - energy transferred")
else:
    print("✗ Detector NOT excited - wave passed through")

if asymmetric:
    print("✓ Pattern asymmetric - one slit affected more")
else:
    print("✗ Pattern symmetric - both slits equal")

print("""
This is the HONEST result from pure substrate physics.
No hardcoded absorption, collapse, or measurement rules.
""")
