import numpy as np
from scipy.io import wavfile
import os

# ==========================================
# 1. THE MATHEMATICAL PARAMETERS
# ==========================================
N_osc = 8             # Lattice size (1D slice of the field)
duration = 30.0       # seconds
sr = 44100            # Audio sample rate
physics_rate = 1000   # How many times per second we solve the diff. equation

# The Anchor Ξ (Inter-universal hash mapped to a harmonic series / C Major chord)
# Phases mapped to frequencies: C3, E3, G3, C4, E4, G4, C5, E5
Xi_frequencies = np.array([130.81, 164.81, 196.00, 261.63, 329.63, 392.00, 523.25, 659.25])
Xi_phases = np.zeros(N_osc) # Target phases for collapse (we will pull frequencies to Xi)

# ==========================================
# 2. SIMULATING THE DIFFERENTIAL EQUATION
# ==========================================
# State: Phase (theta) and Instantaneous Frequency (omega)
theta = np.random.uniform(0, 2*np.pi, N_osc)
omega = np.random.uniform(100, 800, N_osc) # Start with random frequencies (Noise)

dt = 1.0 / physics_rate
total_steps = int(duration * physics_rate)

# Storage for audio generation
freq_history = np.zeros((total_steps, N_osc))

print("Solving the Stochastic Self-Telepathic Field Equation...")
for step in range(total_steps):
    t_norm = step / total_steps  # 0.0 to 1.0 (Time progression)
    
    # --- TIME-VARYING PARAMETERS (The 3 Phases) ---
    # Phase 1: Random Front (0-33%)
    # Phase 2: Telepathic Ignition (33-66%)
    # Phase 3: Emergent Collapse (66-100%)
    
    # Temperature T (Noise) starts high, drops to near zero
    T = max(0.0, 1.0 - (t_norm * 3.0)) 
    
    # Coupling K (Neighbor Pull) starts at 0, ramps up
    K = min(5.0, max(0.0, (t_norm - 0.2) * 10.0))
    
    # Collapse Strength Theta (Telepathy) kicks in at the end
    Theta = min(50.0, max(0.0, (t_norm - 0.6) * 125.0))
    
    # Intrinsic Drive Omega (The rhythmic pulse)
    Omega_drive = 2.0 * np.pi * 0.5 * np.sin(2 * np.pi * 2.0 * t_norm) # 2Hz LFO pulse
    
    # --- THE DIFFERENTIAL EQUATION TERMS ---
    
    # 1. Neighbor Pull (Kuramoto coupling): K * sum(sin(theta_j - theta_i))
    coupling = np.zeros(N_osc)
    for i in range(N_osc):
        for j in range(N_osc):
            if i != j:
                coupling[i] += np.sin(theta[j] - theta[i])
    coupling *= (K / N_osc)
    
    # 2. Telepathic Gradient (Hill-climbing to maximize Phi / minimize dissonance)
    # We pull frequencies toward the harmonic anchor Xi
    telepathic_pull = Theta * (Xi_frequencies - omega)
    
    # 3. Random Front (Thermal Noise)
    noise = np.sqrt(T) * np.random.randn(N_osc) * 50.0
    
    # 4. Intrinsic Drive (Energy input)
    drive = Omega_drive 
    
    # --- EULER INTEGRATION ---
    # Update frequencies (omega) based on the equation
    omega += dt * (telepathic_pull + noise + drive)
    
    # Keep frequencies in audible range
    omega = np.clip(omega, 50, 2000)
    
    # Update phases (theta) for the Kuramoto interference
    theta += dt * (2 * np.pi * omega + coupling)
    theta %= (2 * np.pi)
    
    # Save state for audio rendering
    freq_history[step, :] = omega

# ==========================================
# 3. RENDERING TO AUDIO (Signal Theory)
# ==========================================
print("Rendering audio signal...")
audio_steps = int(duration * sr)
audio_signal = np.zeros(audio_steps)

# Synthesize the sound: Sum of sine waves whose frequencies are governed by the PDE
for i in range(N_osc):
    # Interpolate the physics frequencies to the audio sample rate
    freq_interp = np.interp(np.linspace(0, total_steps, audio_steps), 
                            np.arange(total_steps), 
                            freq_history[:, i])
    
    # Generate phase accumulator for clean sine waves (no clicking)
    phase_acc = np.cumsum(2 * np.pi * freq_interp / sr)
    
    # Add to mix (with slight amplitude envelope based on collapse)
    # Louder in Phase 3 (Collapse) as coherence increases
    amplitude = 0.1 + 0.15 * (np.linspace(0, 1, audio_steps) > 0.66) 
    audio_signal += amplitude * np.sin(phase_acc)

# Normalize and add a tiny bit of vinyl crackle (residual noise) for texture
audio_signal = audio_signal / np.max(np.abs(audio_signal)) * 0.8
crackle = np.random.randn(audio_steps) * 0.005 * (1 - np.linspace(0, 1, audio_steps))
audio_signal += crackle

# Save to WAV
filename = "Telepathic_Field_Sonification.wav"
wavfile.write(filename, sr, (audio_signal * 32767).astype(np.int16))
print(f"Success! The auditory visualization has been saved as '{filename}'")