Za kulisami symulacji #1 — Fale EM: jak zbudowaliśmy…

Punkt wyjścia: czym jest fala EM?

James Clerk Maxwell zunifikował elektryczność i magnetyzm w 1865 roku za pomocą czterech równań. Kluczowym spostrzeżeniem dla tej symulacji jest wynikające z nich równanie falowe: zmieniające się pole elektryczne indukuje pole magnetyczne, i na odwrót. Oba pola gonią się nawzajem w przestrzeni z prędkością światła.

W 2D pracujemy w trybie fali poprzecznej magnetycznej (TM): pole elektryczne ma tylko składową z, Ez, a pole magnetyczne ma składowe x i y, Hx i Hy. Redukuje to równania rotacji Maxwella do trzech równań aktualizujących — wystarczająco prostych, by obliczać je dla każdej komórki siatki w każdej klatce.

FDTD: siatka Yee

Metoda różnic skończonych w dziedzinie czasu (FDTD), wprowadzona przez Kane'a Yee w 1966 roku, dyskretyzuje równania Maxwella na przesuniętej siatce czasoprzestrzennej. Składowe pola elektrycznego i magnetycznego są przesunięte o pół kroku przestrzennego i czasowego, co automatycznie spełnia prawo Gaussa (∇·B = 0) bez żadnego kroku korekcyjnego.

// Równania aktualizacji FDTD trybu TM (2D, siatka jednorodna)
// E i H są przeplecione: najpierw aktualizuje się H, potem E

// Aktualizacja Hx i Hy w kroku czasowym n+½
Hx[i][j] -= (Cz / dy) * (Ez[i][j+1] - Ez[i][j])
Hy[i][j] += (Cz / dx) * (Ez[i+1][j] - Ez[i][j])

// Aktualizacja Ez w kroku czasowym n+1
Ez[i][j] += Cz * (
    (Hy[i][j] - Hy[i-1][j]) / dx -
    (Hx[i][j] - Hx[i][j-1]) / dy
)

// Cz = dt / (ε₀ · μ₀) — liczba Couranta

Warunek stabilności Couranta wymaga dt ≤ dx / (c · √2) dla jednorodnej siatki 2D. Ustawiamy dt = 0,9 · dx / (c · √2) — 90% granicy stabilności dla niewielkiego marginesu bezpieczeństwa. Przy siatce 200×200 i dx = 0,01 m krok czasowy wynosi około 21 ps na aktualizację.

Oś czasu budowy

Dzień 1 — Goła FDTD na zwykłym canvasie

Zaczęliśmy od tablicy JavaScript 100×100, czystego CPU, pojedynczego źródła punktowego. Fala propagowała się poprawnie — kręgi rozchodzące się na zewnątrz. Ale odbijała się od krawędzi idealnie, więc cała domena wypełniała się falami stojącymi w ciągu kilku sekund. Wyraźnie potrzebne były granice absorbujące.

Dzień 2 — Perfectly Matched Layers (PML)

PML to standardowa absorbująca granica FDTD: stratna, gąbczasta strefa wokół domeny, w której pola są tłumione bez odbicia. Implementacja podwoiła złożoność kodu, ale niemal idealnie wyeliminowała odbicia. Domena zaczęła sprawiać wrażenie „otwartej" — fale wychodzą czysto.

Dzień 3 — Przejście na WebGL

Implementacja CPU działała z ~30 FPS na siatce 100×100. Rozszerzenie do 200×200 (minimalna siatka użyteczna do wizualizacji) obniżyło to do 8 FPS. Równania aktualizacji FDTD są aż nadmiarowo równoległe — każdą komórkę można aktualizować niezależnie. Przenieśliśmy pętlę aktualizacji do fragment shadera WebGL zapisującego do pary tekstur zmiennoprzecinkowych (schemat ping-pong). Skoczyliśmy do 60 FPS przy 400×400.

Dzień 4 — Renderowanie i źródła

Zmapowaliśmy pole Ez na niebiesko-czerwoną skalę kolorów (od ujemnych do dodatnich), dodaliśmy sinusoidalne źródło punktowe modulowane obwiednią Gaussa, aby uniknąć początkowego „kliknięcia", i zaimplementowaliśmy podwójną szczelinę dla słynnego w klasie wzoru dyfrakcyjnego. Dodaliśmy źródło dipolowe, źródło fali płaskiej oraz przesuwalny obszar dielektryczny, który użytkownik może umieścić na canvasie.

Dzień 5 — Nakładka pola wektorowego

Nałożenie Hx/Hy jako strzałek za pomocą instancjonowanych odcinków linii w WebGL. Strzałki są podpróbkowane do siatki 20×20, skalowane wielkością pola i obracane zgodnie z kierunkiem pola. Połączenie mapy kolorów Ez i strzałek pola H daje namacalne wyczucie tego, jak oba pola pozostają w kwadraturze.

Najtrudniejsza część: tekstury zmiennoprzecinkowe

WebGL1 nie obsługuje odczytu z tekstury zmiennoprzecinkowej, do której właśnie zapisano — rozszerzenie „render to float texture" (OES_texture_float) istnieje, ale nie jest powszechnie zapisywalne. Rozwiązaniem jest bufor ping-pong: dwa bufory ramki, zamieniające się rolami co klatkę. Shader czyta z bufora A i zapisuje do bufora B; w następnej klatce czyta z B i zapisuje do A.

W Safari na iOS to rozszerzenie w ogóle nie było dostępne, co wymagało fallbacku do 16-bitowych tekstur half-float o zmniejszonym zakresie dynamicznym. Wartości pola są normalizowane przed zapisem i przeskalowywane przy odczycie — niewielka utrata precyzji, ale niewidoczna przy rozdzielczości kolorów ekranu.

Tekstura przenikalności dielektrycznej

Aby umożliwić użytkownikom umieszczanie materiałów (szkła, metalu itd.) w domenie, współczynniki aktualizacji FDTD muszą zmieniać się przestrzennie. Zapisujemy przenikalność ε_r do osobnej tekstury R8 i próbkujemy ją w shaderze aktualizacji. Naciśnięcie przycisku na pasku narzędzi włącza „tryb pędzla" — kolejne zdarzenia myszy malują wartości przenikalności w teksturze materiału za pomocą małego, dodatkowego bufora ramki.

// W shaderze GLSL — odczyt tekstury materiału
float eps_r = texture2D(u_material, vTexCoord).r * 10.0 + 1.0;
// Zakres: [1.0, 11.0] — od powietrza po dielektryk podobny do szkła

float Cz_local = dt / (eps_r * eps0 * mu0 * dx * dx);
// Ten sam wzór aktualizacji, ale Cz zależy od materiału

Zastosowanie w klasie: co działa, a co nie

Po premierze symulacji otrzymaliśmy informację zwrotną od kilku nauczycieli fizyki. Funkcje, które okazały się najcenniejsze w klasie:

Preset dyfrakcji na podwójnej szczelinie — uczniowie mogą bezpośrednio zaobserwować zasadę Huygensa: każda szczelina staje się źródłem wtórnym, a wzór interferencyjny wyłania się z fizyki FDTD bez żadnego ręcznego dostrajania.
Suwak prędkości światła — skalowanie c w dół 10-krotnie sprawia, że propagacja fali jest wystarczająco wolna, by śledzić ją wzrokiem. Czyni to zjawisko namacalnym.
Przełącznik trybu polaryzacji — przełączanie między trybami TE i TM pokazuje, jak dwie polaryzacje światła są po prostu obróconymi wersjami tej samej struktury pola.

Co nie zadziałało: początkowe narzędzie do swobodnego rysowania dielektryka myliło uczniów, ponieważ „rysowanie szkła" wyglądało tak samo jak „rysowanie powietrza", zanim fala do niego dotarła. Dodaliśmy trwałą nakładkę kolorystyczną dla umieszczonych materiałów (bursztynowy odcień = dielektryk, szary = przewodnik), co bardzo pomogło.

Wypróbuj sam: Otwórz symulację Fal EM, kliknij preset „Podwójna szczelina" i przeciągnij suwak odległości między szczelinami. Zobaczysz, jak wzór dyfrakcyjny zmienia się w czasie rzeczywistym — ta sama fizyka co w eksperymencie Thomasa Younga z 1801 roku, działająca w twojej przeglądarce.

Kolejne kroki dla symulacji

2D FDTD w trybie TM to solidne narzędzie dydaktyczne, ale istnieją naturalne rozszerzenia, które chcemy zbudować:

Tryb 3D — wszystkie sześć składowych pola na trójwymiarowej siatce Yee. Wymaga obliczeń WebGPU compute dla przebiegu aktualizacji; dane nie mieszczą się w teksturze 2D przy użytecznych rozmiarach siatki.
Materiały dyspersyjne — model Drudego dla metali (przenikalność zależna od częstotliwości), umożliwiający wizualizację rezonansu plazmonowego.
Transformacja pola bliskiego do dalekiego — obliczanie charakterystyki promieniowania w polu dalekim (charakterystyka anteny) z pola bliskiego zapisanego na „skrzynce Huygensa" wokół źródła.