%reload_ext autoreload
%autoreload 2

from visualization_1D import *

# ============================================================================
# EXAMPLES GALLERY
# ============================================================================

def example_harmonic_oscillator():
    """Example: Harmonic oscillator (integrable system)"""
    x, xi = sp.symbols('x xi', real=True)
    H = xi**2 + x**2
    
    fig, geos, orbits, spec = visualize_symbol(
        H,
        x_range=(-3, 3),
        xi_range=(-3, 3),
        geodesics_params=[
            (0, 0.5, np.pi, 'red'),
            (0, 1.0, np.pi, 'blue'),
            (0, 1.5, np.pi, 'green'),
            (0, 2.0, np.pi, 'purple'),
        ],
        E_range=(0.5, 8),
        hbar=1.0,
        resolution=80
    )
    
    plt.savefig('harmonic_oscillator.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    return fig, geos, orbits, spec


def example_fractional_operator():
    """Example: Fractional operator α=1.5"""
    x, xi = sp.symbols('x xi', real=True)
    alpha = 1.5
    H = xi**alpha
    
    fig, geos, orbits, spec = visualize_symbol(
        H,
        x_range=(-5, 30),
        xi_range=(0.1, 4),
        geodesics_params=[
            (0, 0.5, 5, 'red'),
            (0, 1.0, 5, 'blue'),
            (0, 2.0, 5, 'green'),
            (0, 3.0, 5, 'purple'),
        ],
        E_range=(0.5, 8),
        hbar=0.5,
        resolution=80
    )
    
    plt.savefig('fractional_operator.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    return fig, geos, orbits, spec


def example_double_well():
    """Example: Double well potential (tunneling)"""
    x, xi = sp.symbols('x xi', real=True)
    H = xi**2 + x**4 - 2*x**2
    
    fig, geos, orbits, spec = visualize_symbol(
        H,
        x_range=(-2, 2),
        xi_range=(-3, 3),
        geodesics_params=[
            (-1, 0.5, 3, 'red'),
            (-1, 1.5, 3, 'blue'),
            (0, 1.0, 3, 'green'),
            (1, -0.5, 3, 'orange'),
        ],
        E_range=(-1.0, 5),
        hbar=0.5,
        resolution=80
    )
    
    plt.savefig('double_well.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    return fig, geos, orbits, spec


def example_anharmonic():
    """Example: Anharmonic oscillator"""
    x, xi = sp.symbols('x xi', real=True)
    H = xi**2 + 0.1*x**4
    
    fig, geos, orbits, spec = visualize_symbol(
        H,
        x_range=(-4, 4),
        xi_range=(-4, 4),
        geodesics_params=[
            (0, 0.5, 5, 'red'),
            (0, 1.0, 5, 'blue'),
            (0, 2.0, 5, 'green'),
        ],
        E_range=(0.5, 8),
        hbar=1.0,
        resolution=80
    )
    
    plt.savefig('anharmonic_oscillator.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    return fig, geos, orbits, spec

# ============================================================================
# MAIN
# ============================================================================

if __name__ == "__main__":
    print_theory()
    
    print("\n" + "="*70)
    print("RUNNING EXAMPLES")
    print("="*70)
    
    # Run examples
    print("\n[1/5] Harmonic oscillator...")
    example_harmonic_oscillator()
    
    print("\n[2/5] Fractional operator...")
    example_fractional_operator()
    
    print("\n[3/4] Double well...")
    example_double_well()
    
    print("\n[4/4] Anharmonic oscillator...")
    example_anharmonic()
    
    print("\n✓ All visualizations complete!")

╔══════════════════════════════════════════════════════════════════════╗
║           GEOMETRIC VISUALIZATION OF SYMBOLS                          ║
║        Pseudodifferential Operators & Spectral Analysis               ║
╚══════════════════════════════════════════════════════════════════════╝

GEOMETRIC STRUCTURE
───────────────────

1. HAMILTONIAN SURFACE
   H: T*M → ℝ defines a hypersurface in phase space
   Geodesics are integral curves of the Hamiltonian vector field

2. SYMPLECTIC FOLIATION
   Level sets {H = const} foliate the phase space
   Each leaf is a symplectic manifold

3. HAMILTONIAN FLOW
   Generated by X_H = (∂H/∂ξ, -∂H/∂x)
   Preserves symplectic structure

4. CAUSTICS
   Singularities where Jacobian ∂x/∂ξ₀ vanishes
   Focal points where rays concentrate

SPECTRAL THEORY
───────────────

5. WEYL'S LAW (1911)
   N(E) ~ (2πℏ)^(-d) Vol{H ≤ E}
   Relates spectrum to phase space volume

6. EBK QUANTIZATION (1917)
   ∮_γ p dq = 2πℏ(n + α)
   Periodic orbits quantize the spectrum

7. GUTZWILLER TRACE FORMULA (1971)
   Tr[e^(-iHt/ℏ)] = Σ_γ A_γ e^(iS_γ/ℏ)
   Relates spectrum to periodic orbits

8. BERRY-TABOR (1977)
   Integrable: P(s) ~ e^(-s) (Poisson)
   Level spacings are uncorrelated

9. BGS CONJECTURE (1984)
   Chaotic: P(s) ~ (πs/2) e^(-πs²/4) (Wigner)
   Level repulsion from chaos

10. COLIN DE VERDIÈRE (1985)
    Spectrum determines geometry (partially)
    "Can one hear the shape of a drum?"

11. FERMI'S GOLDEN RULE
    Γ = (2π/ℏ) |V|² ρ(E)
    Transition rate ~ density of states

KEY INSIGHTS
────────────

- Geometry encodes spectrum
- Periodic orbits = spectral oscillations
- Caustics = singularities of wave propagation
- Integrability ↔ Poisson statistics
- Chaos ↔ Wigner statistics

USAGE
─────

from symbol_visualization import visualize_symbol
import sympy as sp

x, xi = sp.symbols('x xi', real=True)
H = xi**2 + x**2  # Your Hamiltonian

fig, geodesics, orbits, spectrum = visualize_symbol(
    symbol=H,
    x_range=(-3, 3),
    xi_range=(-3, 3),
    geodesics_params=[
        (x0, xi0, t_max, color),
        ...
    ],
    E_range=(E_min, E_max),  # Optional: for spectral analysis
    hbar=1.0
)

plt.show()
    

======================================================================
RUNNING EXAMPLES
======================================================================

[1/5] Harmonic oscillator...

=== Spectral Analysis Examples ===
Weyl law estimate of states below E=62.64: N(E) ~ 9.97
Berry-Tabor density of states at E=4.25: ρ(E) ~ 0.000

[2/5] Fractional operator...

[3/4] Double well...

=== Spectral Analysis Examples ===
Weyl law estimate of states below E=15.66: N(E) ~ 4.99
Berry-Tabor density of states at E=2.00: ρ(E) ~ 0.000

[4/4] Anharmonic oscillator...

=== Spectral Analysis Examples ===
Weyl law estimate of states below E=62.64: N(E) ~ 9.97
Berry-Tabor density of states at E=4.25: ρ(E) ~ 0.000

✓ All visualizations complete!