interferencja_fal.py

import numpy as np
import threading
import tkinter as tk
from tkinter import ttk

import matplotlib
matplotlib.use("TkAgg")
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

class SimState:
    def __init__(self):
        self.Lx, self.Ly = 8.0, 4.0
        self.nx, self.ny = 401, 201
        self.x = np.linspace(0, self.Lx, self.nx)
        self.y = np.linspace(0, self.Ly, self.ny)
        self.dx = self.x[1] - self.x[0]
        self.dy = self.y[1] - self.y[0]

        self.c = 1.0
        self.gamma = 0.01
        self.t = 0.0
        self.running = True

        self.f1, self.f2 = 0.65, 0.85
        self.A1, self.A2 = 1.0, 1.0
        self.phi1, self.phi2 = 0.0, np.pi/4
        self.d = 1.2
        self.bc = 'neumann'

        self.u_nm1 = np.zeros((self.ny, self.nx))
        self.u_n   = np.zeros((self.ny, self.nx))
        self.u_np1 = np.zeros((self.ny, self.nx))

        sponge_width = int(0.08 * self.nx)
        self.sponge = np.ones((self.ny, self.nx))
        if sponge_width > 0:
            left = np.linspace(1.0, 0.0, sponge_width)**2
            right = left[::-1]
            self.sponge[:, :sponge_width] *= left
            self.sponge[:, -sponge_width:] *= right

        self.probe_x = self.Lx / 2
        self.probe_y = self.Ly / 2
        self.ix_probe, self.iy_probe = self.nearest_idx(self.probe_x, self.probe_y)

        self.fft_buffer_len = 2048
        self.probe_buffer = np.zeros(self.fft_buffer_len)
        self.probe_ptr = 0
        self.probe_times = []

        self.lock = threading.Lock()

    def compute_dt(self):
        return 0.9 / (self.c * np.sqrt((1/self.dx**2) + (1/self.dy**2)))

    def source_positions(self):
        x1 = self.Lx/2 - self.d/2
        x2 = self.Lx/2 + self.d/2
        y0 = self.Ly/2
        return (x1, y0), (x2, y0)

    def nearest_idx(self, x0, y0):
        ix = int(np.clip(np.round((x0 - self.x[0]) / self.dx), 0, self.nx-1))
        iy = int(np.clip(np.round((y0 - self.y[0]) / self.dy), 0, self.ny-1))
        return ix, iy


def apply_vertical_bc(state, U):
    if state.bc == 'neumann':
        U[0, :]  = U[1, :]
        U[-1, :] = U[-2, :]
    else:
        U[0, :]  = 0.0
        U[-1, :] = 0.0

def add_sources(state, U, dt_local, t_now):
    (sx1, sy1), (sx2, sy2) = state.source_positions()
    ix1, iy1 = state.nearest_idx(sx1, sy1)
    ix2, iy2 = state.nearest_idx(sx2, sy2)
    s1 = state.A1 * np.sin(2*np.pi*state.f1 * t_now + state.phi1)
    s2 = state.A2 * np.sin(2*np.pi*state.f2 * t_now + state.phi2)
    U[iy1, ix1] += (dt_local**2) * s1
    U[iy2, ix2] += (dt_local**2) * s2

def launch_animation_window(state: SimState):
    fig, ax = plt.subplots(figsize=(10, 5))
    im = ax.imshow(state.u_n, extent=[0, state.Lx, 0, state.Ly], origin='lower',
                   cmap='RdBu_r', vmin=-1.5, vmax=1.5, interpolation='bilinear')
    ax.set_title('Pole falowe u(x, y, t)')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_aspect('equal')

    # Markery źródeł
    (sx1, sy1), (sx2, sy2) = state.source_positions()
    s1 = ax.scatter([sx1], [sy1], c='k', marker='x', s=60, label='Źródło 1')
    s2 = ax.scatter([sx2], [sy2], c='k', marker='x', s=60, label='Źródło 2')

    # Marker sondy
    probe_marker = ax.scatter([state.probe_x], [state.probe_y], c='lime', marker='o', s=60, label='Sonda')
    ax.legend(loc='upper right', fontsize=8, framealpha=0.8)

    # Przesuwanie sondy kliknięciem
    def on_click(event):
        if event.inaxes != ax:
            return
        with state.lock:
            state.probe_x = float(np.clip(event.xdata, 0.0, state.Lx))
            state.probe_y = float(np.clip(event.ydata, 0.0, state.Ly))
            state.ix_probe, state.iy_probe = state.nearest_idx(state.probe_x, state.probe_y)
            # Zerowanie bufora po zmianie sondy, by FFT i spektrogram miały spójne dane
            state.probe_buffer.fill(0.0)
            state.probe_ptr = 0
            state.probe_times.clear()
    fig.canvas.mpl_connect('button_press_event', on_click)

    def update(_frame):
        with state.lock:
            if not state.running:
                return (im,)

            dt = state.compute_dt()
            lap = (
                (np.roll(state.u_n, +1, axis=1) - 2*state.u_n + np.roll(state.u_n, -1, axis=1)) / state.dx**2 +
                (np.roll(state.u_n, +1, axis=0) - 2*state.u_n + np.roll(state.u_n, -1, axis=0)) / state.dy**2
            )
            state.u_np1 = (2.0 - state.gamma*dt) * state.u_n - (1.0 - state.gamma*dt) * state.u_nm1 + (state.c*dt)**2 * lap

            apply_vertical_bc(state, state.u_np1)
            add_sources(state, state.u_np1, dt, state.t)

            # Boczne tłumienie
            state.u_np1 *= state.sponge

            # Przesunięcie czasowe
            state.u_nm1, state.u_n = state.u_n, state.u_np1

            # Próbka do sondy (po aktualizacji stanu)
            val = state.u_n[state.iy_probe, state.ix_probe]
            state.probe_buffer[state.probe_ptr % state.fft_buffer_len] = val
            state.probe_ptr += 1
            state.probe_times.append(state.t)
            if len(state.probe_times) > state.fft_buffer_len:
                state.probe_times = state.probe_times[-state.fft_buffer_len:]

            # Aktualizacja krzyżyków źródeł (rozstaw d może się zmieniać)
            (sx1, sy1), (sx2, sy2) = state.source_positions()
            s1.set_offsets([[sx1, sy1]])
            s2.set_offsets([[sx2, sy2]])

            # Aktualizacja markera sondy
            probe_marker.set_offsets([[state.probe_x, state.probe_y]])

            # Auto-skalowanie kolorów
            vmax = np.percentile(np.abs(state.u_n), 99.5) + 1e-6
            im.set_clim(-vmax, vmax)
            im.set_data(state.u_n)

            state.t += dt

        return (im,)

    ani = FuncAnimation(fig, update, interval=15, blit=False)
    return fig, ani



def launch_analysis_window(state: SimState):
    fig, (ax_probe, ax_fft, ax_spec) = plt.subplots(3, 1, figsize=(8, 9))
    fig.suptitle("Analiza sygnału w sondzie")

    # Wykres czasowy
    line_probe, = ax_probe.plot([], [], lw=1.2)
    ax_probe.set_title("Sygnał w punkcie pomiarowym")
    ax_probe.set_xlabel("Czas [s]")
    ax_probe.set_ylabel("u")

    # Widmo FFT
    line_fft, = ax_fft.plot([], [], lw=1.2)
    ax_fft.set_title("Widmo amplitudy (FFT)")
    ax_fft.set_xlabel("Częstotliwość [Hz]")
    ax_fft.set_ylabel("Amplituda")

    # Spektrogram
    ax_spec.set_title("Spektrogram")
    ax_spec.set_xlabel("Czas [s]")
    ax_spec.set_ylabel("Częstotliwość [Hz]")

    def update(_frame):
        with state.lock:
            times = np.array(state.probe_times, dtype=float)
            if times.size >= 8:
                # Pobierz zsynkronizowane wartości z bufora
                idxs = np.arange(state.probe_ptr - times.size, state.probe_ptr)
                vals = state.probe_buffer.take(idxs, mode='wrap')

                # Wykres czasowy
                line_probe.set_data(times, vals)
                ax_probe.set_xlim(times[0], times[-1])
                vspan = float(np.max(np.abs(vals)))
                vspan = max(1e-6, vspan)
                ax_probe.set_ylim(-1.2*vspan, 1.2*vspan)

                # FFT z okna danych
                dt = state.compute_dt()
                win = np.hanning(len(vals))
                sigw = vals * win
                freqs = np.fft.rfftfreq(len(sigw), d=dt)
                amps = (2.0 / np.sum(win)) * np.abs(np.fft.rfft(sigw))
                line_fft.set_data(freqs, amps)
                ax_fft.set_xlim(0, max(2.5, max(state.f1, state.f2) * 3.0))
                ymax = float(np.percentile(amps, 99.5)) * 1.2 if np.any(amps > 0) else 1.0
                ax_fft.set_ylim(0, ymax)

                # Spektrogram (odświeżany rzadziej dla płynności)
                if state.probe_ptr % 20 == 0:
                    ax_spec.cla()
                    ax_spec.set_title("Spektrogram")
                    ax_spec.set_xlabel("Czas [s]")
                    ax_spec.set_ylabel("Częstotliwość [Hz]")
                    Fs = 1.0 / dt
                    # Użyj tego samego okna „vals” dla spójności
                    ax_spec.specgram(vals, NFFT=256, Fs=Fs, noverlap=192, cmap='magma')

        return (line_probe, line_fft)

    ani = FuncAnimation(fig, update, interval=300, blit=False)
    return fig, ani



def launch_control_window(state: SimState):
    root = tk.Tk()
    root.title("Sterowanie symulacją fal (Tkinter)")

    def add_slider(label, from_, to, getter, setter, fmt="{:.3f}"):
        frame = ttk.Frame(root)
        frame.pack(fill='x', padx=8, pady=4)

        ttk.Label(frame, text=label, width=22).pack(side='left')

        val_var = tk.DoubleVar(value=getter())

        # Etykieta wartości po prawej
        val_label = ttk.Label(frame, text=fmt.format(val_var.get()), width=10, anchor='e')
        val_label.pack(side='right')

        def on_slide(v):
            v = float(v)
            with state.lock:
                setter(v)
            val_label.config(text=fmt.format(v))

        scale = ttk.Scale(frame, from_=from_, to=to, orient='horizontal', command=on_slide)
        scale.set(val_var.get())
        scale.pack(side='left', fill='x', expand=True, padx=8)

        return scale, val_label

    def add_button(text, cmd):
        def wrapped():
            with state.lock:
                cmd()
        btn = ttk.Button(root, text=text, command=wrapped)
        btn.pack(fill='x', padx=8, pady=4)
        return btn

    def add_radio(label, options, getter, setter):
        frame = ttk.Frame(root)
        frame.pack(fill='x', padx=8, pady=4)
        ttk.Label(frame, text=label, width=22).pack(side='left')
        var = tk.StringVar(value=getter())
        def on_change():
            with state.lock:
                setter(var.get())
        for opt in options:
            rb = ttk.Radiobutton(frame, text=opt.capitalize(), value=opt, variable=var, command=on_change)
            rb.pack(side='left', padx=4)
        return frame

    # Slidery parametrów z etykietami wartości
    add_slider("Częstotliwość f1 [Hz]", 0.05, 2.0, lambda: state.f1, lambda v: setattr(state, 'f1', v), fmt="{:.3f}")
    add_slider("Częstotliwość f2 [Hz]", 0.05, 2.0, lambda: state.f2, lambda v: setattr(state, 'f2', v), fmt="{:.3f}")
    add_slider("Amplituda A1", 0.0, 2.0, lambda: state.A1, lambda v: setattr(state, 'A1', v), fmt="{:.3f}")
    add_slider("Amplituda A2", 0.0, 2.0, lambda: state.A2, lambda v: setattr(state, 'A2', v), fmt="{:.3f}")
    add_slider("Faza phi1 [rad]", -np.pi, np.pi, lambda: state.phi1, lambda v: setattr(state, 'phi1', v), fmt="{:.2f}")
    add_slider("Faza phi2 [rad]", -np.pi, np.pi, lambda: state.phi2, lambda v: setattr(state, 'phi2', v), fmt="{:.2f}")
    add_slider("Rozstaw d", 0.1, 0.9*state.Lx, lambda: state.d, lambda v: setattr(state, 'd', v), fmt="{:.3f}")
    add_slider("Prędkość c", 0.2, 3.0, lambda: state.c, lambda v: setattr(state, 'c', v), fmt="{:.3f}")
    add_slider("Tłumienie gamma", 0.0, 0.1, lambda: state.gamma, lambda v: setattr(state, 'gamma', v), fmt="{:.4f}")

    add_radio("Ściany góra/dół", ['neumann', 'dirichlet'], lambda: state.bc, lambda v: setattr(state, 'bc', v))

    # Przyciski sterujące
    add_button("Start / Stop", lambda: setattr(state, 'running', not state.running))

    def do_reset():
        state.u_nm1.fill(0.0)
        state.u_n.fill(0.0)
        state.u_np1.fill(0.0)
        state.t = 0.0
        state.probe_buffer.fill(0.0)
        state.probe_ptr = 0
        state.probe_times.clear()
    add_button("Reset", do_reset)

    root.geometry("520x560")
    root.mainloop()



if __name__ == "__main__":
    state = SimState()

    # Panel sterowania w wątku pobocznym
    threading.Thread(target=launch_control_window, args=(state,), daemon=True).start()

    # Okna Matplotlib w głównym wątku
    fig_anim, ani_anim = launch_animation_window(state)
    fig_an, ani_an = launch_analysis_window(state)

    plt.show()