Počítačové simulace vybraných fyzikálních jevů

Přírodovědecká fakulta Univerzity Palackého v Olomouci Katedra optiky

Počítačové simulace vybraných fyzikálních jevů Bakalářská práce

Karel Lemr

Olomouc

duben 2006

Prohlašuji, že jsem tuto práci zpracoval samostatně a pouze s použitím uvedené literatury.

Děkuji prof. RNDr. Jiřímu Bajerovi, CSc. za cenné rady a konzultace.

Obsah 1 Úvod

1

2 Počítačová simulace difrakčních obrazců 2.1 Teoretický základ . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Algoritmické řešení apletu . . . . . . . . . . . . . . . . . . . . 2.2.1 Významné proměnné použité v algoritmu . . . . . . . 2.2.2 Základní algoritmus . . . . . . . . . . . . . . . . . . . 2.2.3 Optimalizovaný algoritmus . . . . . . . . . . . . . . . 2.2.4 Vykreslení difrakčního obrazce . . . . . . . . . . . . . 2.2.5 Porovnání časové náročnosti představených algoritmů 2.3 Návod k použití apletu . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

2 . 2 . 4 . 4 . 5 . 6 . 10 . 11 . 11

3 Počítačová simulace šíření světelného paprsku 15 3.1 Teoretický základ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Algoritmické řešení apletu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.3 Návod k použití apletu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4 Počítačová simulace Lissajousových 4.1 Teoretický základ . . . . . . . . . . 4.2 Algoritmické řešení apletu . . . . . 4.3 Návod k použití apletu . . . . . . . 5 Závěr

křivek 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 25

Počítačové simulace vybraných fyzikálních jevů

1

11. června 2006

Úvod

V mnoha přírodních vědách se objevují problémy, které je z různých důvodů vhodné řešit pomocí simulace. Tyto důvody mohou být časové, ekonomické, ale také například bezpečnostní. V minulosti proto vznikaly důmyslné přístroje, které by k řešení alespoň z části napomohly. S příchodem moderní výpočetní techniky se vědcům otevřely nové možnosti, jak řešit celou řadu problémů za pomoci počítačových simulací. Simulace je nejen mocným nástrojem pro ověření platnosti teoreticky odvozených vztahů, ale je také vynikající pomůckou při výuce. Znázornění daného fyzikálního jevu prostřednictvím simulace usnadňuje jeho pochopení a osvojení si pravidel, kterými se řídí. Zmíněné důvody mne motivovaly k tomu, abych svou pozornost věnoval právě počítačové simulaci. Jako student fyziky jsem si položil za cíl vytvořit počítačové programy simulující vybrané fyzikální jevy. Před zahájením práce jsem se musel rozhodnout, který programovací jazyk bude nejvhodnější. Po pečlivé úvaze jsem nakonec vybral jazyk Java vyvíjený firmou Sun Microsystems (viz obrázek 1). Mezi

Obrázek 1: Logo programovacího jazyka Java jeho hlavní přednosti patří nezávislost na platformě (operačním systému), kterou uživatel na svém počítači provozuje. Další významnou předností tohoto programovacího jazyka je dostupnost vývojových nástrojů zdarma přímo na webových stránkách výrobce (http://java.sun.com). Třetí poměrně zásadní předností Javy je integrovatelnost aplikací přímo do WWW stránek formou tzv. Java apletů. Uživatel není nucen program stahovat ani instalovat, ale stačí, když pomocí webového prohlížeče navštíví stránku obsahující program. Formu Java apletů vložených do WWW stránek jsem zvolil pro všechny své simulační programy předkládané v této práci. Každé simulaci je věnována samostatná kapitola členěná na tři sekce. V první sekci, nazvané teoretický základ, nastíním fyzikální pozadí studovaného jevu a uvedu matematický aparát nezbytný k sestavení vhodného počítačového algoritmu. Druhá sekce se věnuje popisu samotného algoritmu a seznamuje čtenáře s jeho chodem. Ve třetí sekci popisuji návod k použití apletu, abych uživateli představil danou simulaci a způsoby jejího ovládání. Jako náměty pro počítačové simulace slouží tři fyzikální jevy. V první části této práce se budu věnovat simulaci difrakce světla, ve druhé části se zaměřím na simulaci chodu světelných paprsků a ve třetí části představím simulaci Lissajousových křivek. Kromě textu, který právě čtete, tvoří mou práci i přiložený disk CD, na kterém naleznete nejen aplety samotné, ale také kompletní zdrojové kódy, které mohou sloužit k hlubšímu prostudování použitých algoritmů. Obsah CD je také k dispozici na webových stránkách http://optics.upol.cz/lemr/java/.

1


2 2.1

11. června 2006

Počítačová simulace difrakčních obrazců Teoretický základ

Narazí-li vlna při svém šíření na překážku, začne se ohýbat i do míst geometrického stínu. Popisovanému jevu říkáme difrakce a je jedním z typických projevů vlnění. Experimenty ukazují, že difrakci je možné pozorovat i u světla. Tento fakt podporuje chápání světla jako elektromagnetického vlnění [1]. Studiem difrakce se zabýval italský fyzik Francesco Maria Grimaldi (2. 4. 1618 – 28. 12. 1663). Ve své knize, publikované v roce 1665 (tedy dva roky po autorově smrti), popisuje difrakci světla na úzké štěrbině. Nezávisle na Grimaldiho objevech vydal roku 1690 dánský fyzik Christian Huygens (14. 4. 1629 – 8. 6. 1695) práci zabývající se šířením světla a nesoucí název Traité de la lumi`ere. Huygens v ní ukázal, že v každém okamžiku lze body, do kterých čelo vlny dospěje, považovat za elementární zdroje vlnění. Výsledný tvar vlnoplochy v následujícím okamžiku je poté určen superpozicí kulových vlnoploch od těchto elementárních zdrojů. Touto metodou dokázal Huygens vysvětlit například odraz světla nebo lom světla. Roku 1818 se Augustinu Jeanu Fresnelovi (10. 5. 1788 – 14. 7. 1827) podařilo s využitím Huygensova principu a interference popsat také difrakční jevy. Proto se v této souvislosti v literatuře často hovoří o Huygens-Fresnelově principu. Přesnou matematickou podobu dal tomuto principu Gustav Robert Kirchhoff (12. 3. 1824 – 17. 10. 1887). Odvodil integrální vztah " ! # 2πis 2πis 1 { ∂ e λ e λ ∂U U (P ) = U − dS (1) 4π ∂n s s ∂n S

známý dnes jako Kirchhoffův integrál, ve kterém ukázal, že komplexní amplitudu vlny U (P ) v libovolně zvoleném bodě P je možné získat ze znalosti rozložení komplexní amplitudy a její derivace na uzavřené ploše obklopující tento bod. Ve vztahu (1) označuje λ vlnovou délku záření a s vzdálenost bodu P od elementu dS na ploše S (viz obrázek 1). Pomocí tohoto vztahu lze přímo počítat difrakci monochromatických vln. V případě, kdy studujeme vlnění polychromatické, je nejprve třeba rozložit vlnu do spektra a poté řešit difrakci pro každou složku spektra zvlášť. V dalších úvahách se ovšem omezíme pouze na monochromatické záření. ~n dS

s

P S

Obrázek 1: Pomocí Kirchhoffova integrálu je možné určit komplexní amplitudu vlny v bodě P ze znalosti rozložení komplexní amplitudy a její derivace na ploše, která bod P uzavírá.

Ukazuje se, že přesné analytické provedení Kirchhoffova integrálu je možné pouze v několika vybraných případech. Jedním z mála takovýchto případů je difrakce světla na hraně elektricky dokonale vodivé poloroviny. Řešení tohoto problému představil roku 1896 německý fyzik Arnold Sommerfeld (5. 12. 1868 – 26. 4. 1951). Ve většině difrakčních úloh je však nutné využívat různých aproximací [2, 3, 4, 5, 6]. Často používanou aproximaci navrhl další německý fyzik Joseph Fraunhofer (6. 3. 1787 – 7. 6. 1826). Předpokládejme (viz obrázek 2) zdroj monochromatického záření Z, aperturu A a stínítko S. Na apertuře zavedeme souřadné osy ξ a η, na stínítku pak osy x a y. Pro Fraunhoferovu aproximaci je 2


11. června 2006

charakteristické pokládat zdroj i stínítko do vzdáleností mnohonásobně převyšujících příčné rozměry apertury (teoreticky do nekonečna). Vlny dopadající ze zdroje na aperturu i vlny, které z apertury dopadnou na stínítko, lze při této představě považovat za vlny rovinné [2, 3, 6, 7, 8]. η

y x

ξ s

r

Z A

S

Obrázek 2: Schéma realizace Fraunhoferovy difrakce. Dalším důležitým krokem je popis apertury. Zde zavádíme obecně komplexní funkci f (ξ, η). Naše úvahy však omezíme pouze na apertury tvořené průsvitnými otvory v neprůsvitném materiálu. Neuvažujeme apertury modulující fázi vlny. Funkce f (ξ, η) tak nabude tvaru 1 v místech otvorů, f (ξ, η) = (2) 0 jinde. Nyní již můžeme formulovat vztah pro Fraunhoferovu difrakci definující komplexní amplitudu U (x, y) v libovolném bodě stínítka { 2πi U (x, y) = U0 f (ξ, η)e λr (xξ+yη) dξdη. (3) A

Integraci provádíme přes celou aperturu a U0 označuje konstantu spojenou s vlastnostmi zdroje. Po zavedení substituce x y νx = a νy = (4) λr λr lze na Fraunhoferovu difrakci nahlížet jako na dvourozměrnou integrální Fourierovu transformaci funkce f (ξ, η). Tento přístup nám umožní využít při výpočtu difrakce známých matematických vztahů pro Fourierovu transformaci [2, 7, 9]. Je-li matematický popis funkce f (ξ, η) příliš komplikovaný, nahrazujeme analytický výpočet integrálu numerickými metodami. Proto přejdeme od integrální podoby Fourierovy transformace k podobě diskrétní. Komplexní amplitudu U (νx , νy ) pak lze určit následujícím vztahem U (νx , νy ) = U0

N −1 N −1 X X

0

f (ξ, η)e2πi(νx ξ+νy η)/N ,

(5)

ξ=0 η=0

kde jsme aperturu rozložili v každé ose na N intervalů, přes které sčítáme. Pro účely řešení difrakce za použití výpočetní techniky bude nezbytné také podobným způsobem rozdělit i stínítko na konečný počet bodů, které bude charakterizovat dvojice indexů νx a νy nabývajících hodnot 0, 1, 2, . . . , N 0 − 1. Simulace popisovaná v této práci pracuje s hodnotami N = 50 a N 0 = 200. Programovací jazyk Java není schopen pracovat s komplexními čísly, proto je nutné rozdělit výpočet komplexní amplitudy na reálnou a imaginární složku. S využitím Eulerových vztahů dostaneme <[U (νx , νy )] = U1

N −1 N −1 X X

f (ξ, η) cos 2π(νx ξ + νy η)/N 0

ξ=0 η=0

3

(6)


11. června 2006

pro reálnou část komplexní amplitudy a =[U (νx , νy )] = U2

N −1 N −1 X X

f (ξ, η) sin 2π(νx ξ + νy η)/N 0

(7)

ξ=0 η=0

pro imaginární část komplexní amplitudy. Vztahy (6) a (7) jsou platné pouze za předpokladu, že funkce f (ξ, η) je reálnou funkcí. Jak již bylo zmíněno, omezujeme se pouze na situace, kdy lze funkci f (ξ, η) vyjádřit vztahem (2). Tím je platnost vzorců (6) a (7) zaručena [10]. Pozorovanou veličinou na stínítku nebude komplexní amplituda, ale intenzita světla I(νx , νy ). Tu ze znalosti komplexní amplitudy dopočítáme pomocí vztahu I(νx , νy ) = |U (νx , νy )|2 = <[U (νx , νy )]2 + =[U (νx , νy )]2 .

2.2 2.2.1

(8)

Algoritmické řešení apletu Významné proměnné použité v algoritmu

Než začneme vlastním rozborem výpočetních algoritmů, věnujme pozornost několika významným proměnným, které jsou v programu použity. Zásadní roli hraje proměnná s názvem velikostBodu, která určuje vlastnosti stínítka. Jak již bylo zmíněno v teoretické části, stínítko je pro účely počítačové simulace rozděleno na množinu diskrétních bodů. Pro každý bod je výpočet intenzity prováděn samostatně. Minimální velikost jednoho bodu je 1×1 px a oblast stínítka určená pro vykreslení difrakčního obrazce má rozměry 200×200 px. Odtud vidíme, že simulace může obsahovat až 40 000 nezávislých bodů. Výpočet intenzity je pro takto vysoký počet prvků časově náročný. Pro zrychlení výpočtu je možné zvětšit velikost zobrazovaných bodů na celistvý násobek rozměru pixelu. K tomu je právě využita proměnná velikostBodu, neboť určuje výšku a šířku jednoho bodu v pixelech1 . Je zřejmé, že s rostoucí hodnotou velikostBodu bude klesat počet bodů na stínítku a tím i časová náročnost výpočtu difrakčního obrazce. S tím ovšem souvisí snížení celkového rozlišení. Proměnná velikostBodu určuje poměr mezi časovou náročností a kvalitou zobrazení. Uživatel má možnost volit hodnotu této proměnné tak, aby mohl dosáhnout požadovaného poměru. Další velmi důležitou proměnnou pro celý algoritmus je pole typu boolean s názvem pruchod[][]. Toto dvourozměrné pole může v každé své složce nabývat hodnot true nebo false a slouží k popisu apertury. Hodnoty true odpovídají otevřeným otvorům a hodnoty false otvorům uzavřeným. Pole pruchod[][] je programátorským vyjádřením funkce f (ξ, η) ve smyslu její definice (2). Hodnoty true pole pruchod[][] odpovídají hodnotám 1 funkce f (ξ, η). Podobně hodnoty false odpovídají hodnotám 0. Trojice dvourozměrných polí typu double s názvy re_amplituda[][], im_amplituda[][] a intenzita[][] určují postupně reálnou složku amplitudy, imaginární složku amplitudy a intenzitu jednotlivých bodů stínítka. V průběhu výpočtu jsou hodnoty těchto proměnných spočteny pro každý jednotlivý bod. Dále je do výpočtů zavedena proměnná s názvem intenzita_prum (typu double), pomocí které je určena průměrná intenzita na stínítku. Tato informace je dále využita ke škálování intenzity tak, aby bylo dosaženo obrazu s co možná nejoptimálnějším jasem. Uživatel má k dispozici volbu vlnové délky použitého záření. Ještě před výpočtem difrakčního obrazce je této vlnové délce třeba přiřadit barvu v systému RGB2 . Tento výpočet barvy z vlnové délky zajišťuje jednoduchý algoritmus využívající tabulku barev a vlnových délek [11]. Tabulka obsahuje vlnové délky pouze několika vybraných barev a převodní algoritmus je třeba považovat pouze za orientační. Informace o barvě jsou pro další použití uloženy v poli typu integer nesoucím název RGB[]. 1

Celkový počet pixelů v jednom bodu je tedy velikostBodu2 . Programovací jazyk Java popisuje barvy prostřednictvím systému RGB, kdy každé barvě vyhradí jeden byte informace. Celkový počet různých barev je 16 777 216. 2

4


2.2.2

11. června 2006

Základní algoritmus

Program předkládaný v této práci má implementovány dva algoritmy pro výpočet difrakčních obrazců, mezi kterými si může uživatel vybrat. První algoritmus (zdrojový kód 1) vychází přímo ze vztahů (6), (7) a (8). Provádí neoptimalizovanou diskrétní Fourierovu transformaci (zkráceně DFT). Zdrojový kód 1: Základní algoritmus 1 2 3 4 5 6 7 8 9 10

( i n t x = 0; x < 200; x = x + { f o r ( i n t y = 0; y < 200; y = y + { re_amplituda [ x ][ y ] = 0; im_amplituda [ x ][ y ] = 0; f o r ( i n t j = 0; j < 50; { f o r ( i n t k = 0; k < 50; {

for

velikostBodu ) velikostBodu )

j ++) k ++)

11

if

12 13

14

15

( pruchod [ j ][ k ] == t r u e ) { re_amplituda [ x ][ y ] = re_amplituda [ x ][ y ] + Math . cos (2* Math . PI *(( x -100) * j + ( y -100) * k ) / z /200*555/ lambda ) ; im_amplituda [ x ][ y ] = im_amplituda [ x ][ y ] + Math . sin (2* Math . PI *(( x -100) * j + ( y -100) * k ) / z /200*555/ lambda ) ; }

} }

16 17 18 19

20 21 22 23 24 25

26

intenzita [ x ][ y ] = ( Math . pow ( re_amplituda [ x ][ y ] ,2) + Math . pow ( im_amplituda [ x ][ y ] ,2) ) /256; intenzita_prum = intenzita_prum + intenzita [ x ][ y ]; i f ( x == y ) { g . setColor ( Color . black ) ; g . fillRect (275 ,75 ,325 ,125) ; g . setColor ( Color . white ) ; g . drawString ( String . valueOf ( Math . round (( d o u b l e ) x /2) ) + " % " ,295 ,95) ; }

27 28 29

} } První dva cykly (řádek 1 a 3) zajistí, aby algoritmus prošel všechny body na stínítku. Těmito dvěma cykly nás provádí indexy x a y, které odpovídají proměnným x a y znázorněným na obrázku 2. Všimněme si, že použité indexy x a y probíhají od nuly do 199 po skocích o velikosti velikostBodu. To je v souladu s výše zmíněnou funkcí proměnné velikostBodu. Skutečně vidíme, jak s její rostoucí hodnotou klesá počet proběhnutí obou cyklů a tím i časová náročnost výpočtu. Pro fixní hodnoty indexů x a y, které určují konkrétní bod na stínítku, je nyní třeba spočítat příspěvky komplexních amplitud od všech otvorů apertury. K tomu slouží další dva cykly (řádky 7 a 9), ve kterých sčítáme přes indexy j resp. k odpovídající proměnným ξ resp. η z obrázku 2. V cyklech je nejprve podmínkou na řádku 12 ověřeno, zda je otvor na pozici (j, k) otevřen. Pokud tomu tak je, 5


11. června 2006

přičte se k reálné i imaginární části komplexní amplitudy příspěvek od tohoto otvoru (řádky 13 a 14). Při přičtení se algoritmus opírá o vztahy (6) a (7), ve kterých je ovšem nutné provést několik úprav. Teoretické vztahy předpokládají, že optická osa prochází počátky souřadných systémů na apertuře i na stínítku. V programovacím jazyce Java však počátky indexování dvourozměrných polí neleží v jejich středu, ale v jejich levém horním rohu. S přihlédnutím k rozměrům stínítka 200×200 px a apertury 50×50 px zavedeme následující korekce x y j k

−→ −→ −→ −→

x − 100, y − 100, j − 25, k − 25.

(9)

Využijeme-li vztahu (5), můžeme s proměnnými použitými v programu psát U (x, y) = U0

49 X 49 X

f (j, k)e2πi(x−100)(j−25)/200 e2πi(y−100)(k−25)/200

(10)

j=0 k=0 a odtud U (x, y) = e

−πi(x+y−200)/4

U0

49 X 49 X

f (j, k)e2πi(x−100)j/200 e2πi(y−100)k/200 .

(11)

j=0 k=0 První člen na pravé straně ve vztahu (11) byl vytknut před znaménka sumy, protože nezávisí ani na j ani na k. Představuje pouze fázový posun závislý na poloze bodu (x, y). Je zřejmé, že se při výpočtu intenzity neuplatní a je tedy možné jej při konstrukci algoritmu vypustit. Dále si všimněme, že se v argumentech funkcí Math.cos a Math.sin na řádcích 13 a 14 vyskytuje ve jmenovateli proměnná z. Ta zajišťuje funkci lupy, kterou si může uživatel zvětšit prostřední část difrakčního obrazce. Její matematický význam lze chápat jako škálování souřadných os x a y. S rostoucí hodnotou z se bude zvětšovat měřítko na těchto osách. Výsledkem bude zvětšený výřez obrazu poblíž optické osy. V souladu s teorií se na řádcích 13 a 14 v argumentech funkcí Math.cos a Math.sin objevuje ve jmenovateli hodnota vlnové délky lambda. Navíc je zde ovšem korekční člen v čitateli s hodnotou 555. Ten zajišťuje pro rozměry použité v programu přijatelné škálování vlnové délky. Vidíme, že pro vlnovou délku 555 nm (lambda = 555) nedojde k žádné změně difrakčního obrazce. Pro větší vlnové délky bude difrakční obrazec více roztáhlý a pro menší vlnové délky bude roztáhlý méně. Tento efekt je v souladu s experimentálně pozorovanou skutečností. Po ukončení cyklů započatých na řádcích 7 a 9 je z reálné a komplexní složky amplitudy podle vztahu (8) vypočtena celková intenzita světla. Následně je tato intenzita přičtena k proměnné intenzita_prum za účelem výpočtu průměrné intenzity. Na řádcích 21 – 26 je zařazen krátký kód sloužící k zobrazení aktuálního průběhu výpočtu. Stav probíhajícího výpočtu se uživateli znázorňuje jako hodnota v procentech. 2.2.3

Optimalizovaný algoritmus

Druhý algoritmus (zdrojový kód 2), se kterým program pracuje, se inspiruje výhodami tzv. rychlé Fourierovy transformace (FFT) [12, 13, 14]. Algoritmus se však musí vypořádat nejen s dvourozměnou transformací, ale také se škálováním, neboť stínítko a apertura nemají stejné rozměry a není možné využít běžných algoritmů pro rychlou Fourierovu transformaci 1:1. Ukažme si nyní, jak lze optimalizovat výpočet difrakce zadané vztahem (11). Zavedeme-li proměnnou M = 25, můžeme psát U (x, y) = e

πi(x+y−8M )/4

U0

2M −1 2M −1 X X

f (j, k)e2πi(x−4M )j/8M e2πi(y−4M )k/8M ,

j=0 k=0 6

(12)


11. června 2006

tedy U (x, y) ∝

2M −1 2M −1 X X

f (j, k)e2πi(x−4M )j/8M e2πi(y−4M )k/8M .

(13)

j=0 k=0 Po roznásobení závorek v argumentech exponenciálních funkcí a využití identity eπiα ≡ (−1)α ,

pro ∀α ∈ Z

(14)

lze vztah (13) zapsat jako U (x, y) ∝

2M −1 X

k πiyk/4M

(−1) e

k=0

2M −1 X

(−1)j f (j, k)eπixj/4M .

(15)

j=0

Sumy v předchozím vztahu lze rozdělit na sumy sčítající pouze přes sudé a pouze přes liché hodnoty indexů j a k. Touto úpravou obdržíme vyjádření U (x, y)      M M −1 M −1 −1 −1  M X X X  X  πix(2j+1)/4M  πiy(2k+1)/4M  πixj/2M πiyk/2M − f (2j + 1, k)e − e f (2j, k)e U (x, y) ∝  e .       j=0  k=0  j =0 k =0 {z } | {z } | {z } | {z } | G1

G2

F1

F2

(16) Po roznásobení závorek můžeme psát3 F1 ~ G 1 =

PM −1 PM −1 πi(xj+yk)/2M j=0 k=0 f (2j, 2k)e

F1 ~ G2 = eπiy/4M

PM −1 PM −1 πi(xj+yk)/2M j=0 k=0 f (2j, 2k + 1)e

eπix/4M

PM −1 PM −1 πi(xj+yk)/2M j=0 k=0 f (2j + 1, 2k)e

(17) F2 ~ G 1 =

F2 ~ G2 = eπi(x+y)/4M

PM −1 PM −1 πi(xj+yk)/2M . j=0 k=0 f (2j + 1, 2k + 1)e

S využitím vztahu (16) a rovnic (17), je možné vyjádřit komplexní amplitudu na stínítku následujícím způsobem U (x, y) ∝ F1 ~ G1 − F1 ~ G2 − F2 ~ G1 + F2 ~ G2 U (x + 4M, y) ∝ F1 ~ G1 − F1 ~ G2 + F2 ~ G1 − F2 ~ G2 (18) U (x, y + 4M ) ∝ F1 ~ G1 + F1 ~ G2 − F2 ~ G1 − F2 ~ G2 U (x + 4M, y + 4M ) ∝ F1 ~ G1 + F1 ~ G2 + F2 ~ G1 + F2 ~ G2 . Díky tomuto vyjádření můžeme snížit počet cyklů, kterými probíhají indexy x a y z 200 na 100 a počet cyklů proměnných j a k z 50 na 25. Tato úprava představuje významné urychlení výpočtů. Její algoritmická implementace je uvedena ve zdrojovém kódu 2. Zdrojový kód 2: Optimalizovaný algoritmus 1

for

2 3

( i n t x = 0; x < 100; x = x + velikostBodu ) {

Symbol ~ v následujících vztazích neoznačuje klasický součin, ale součin ve smyslu pravých stran výrazů.

7


3 4 5 6 7 8 9 10 11 12 13

for

( i n t y = 0; y { F1G1re = F1G1im = F2G1re = F2G1im = F1G2re = F1G2im = F2G2re = F2G2im = for ( int

14 15 16 17

18 19 20

for

11. června 2006

< 100; y = y + velikostBodu ) 0; 0; 0; 0; 0; 0; 0; 0; k = 0; k < 25; k ++)

{ ( i n t j = 0; j < 25; j ++) { i f ( pruchod [2* j ][2* k ] || pruchod [2* j +1][2* k ] || pruchod [2* j ][2* k +1] || pruchod [2* j +1][2* k +1]) { FxGxre = Math . cos ( Math . PI *( j * x + k * y ) /50) ; FxGxim = Math . sin ( Math . PI *( j * x + k * y ) /50) ; }

21 22

if

( pruchod [2* j ][2* k ] == t r u e ) { F1G1re = F1G1re + FxGxre ; F1G1im = F1G1im + FxGxim ; }

if

( pruchod [2* j ][2* k +1] == t r u e ) { F1G2re = F1G2re - FxGxre ; F1G2im = F1G2im - FxGxim ; }

if

( pruchod [2* j +1][2* k ] == t r u e ) { F2G1re = F2G1re - FxGxre ; F2G1im = F2G1im - FxGxim ; }

if

( pruchod [2* j +1][2* k +1] == t r u e ) { F2G2re = F2G2re + FxGxre ; F2G2im = F2G2im + FxGxim ; }

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

} }

44 45

46

F1G2reB = F1G2re * Math . cos ( Math . PI * y /100) - F1G2im * Math . sin ( Math . PI * y /100) ; F1G2imB = F1G2re * Math . sin ( Math . PI * y /100) + F1G2im * Math . cos ( Math . PI * y /100) ;

47 48

F2G1reB = F2G1re * Math . cos ( Math . PI * x /100) - F2G1im * Math . sin ( Math . PI * x /100) ;

8


49

11. června 2006

F2G1imB = F2G1re * Math . sin ( Math . PI * x /100) + F2G1im * Math . cos ( Math . PI * x /100) ;

50 51

52

F2G2reB = F2G2re * Math . cos ( Math . PI *( x + y ) /100) - F2G2im * Math . sin ( Math . PI *( x + y ) /100) ; F2G2imB = F2G2re * Math . sin ( Math . PI *( x + y ) /100) + F2G2im * Math . cos ( Math . PI *( x + y ) /100) ;

53 54 55

re_amplituda [ x ][ y ] = F1G1re + F1G2reB + F2G1reB + F2G2reB ; im_amplituda [ x ][ y ] = F1G1im + F1G2imB + F2G1imB + F2G2imB ;

56 57 58

re_amplituda [ x +100][ y ] = F1G1re + F1G2reB - F2G1reB - F2G2reB ; im_amplituda [ x +100][ y ] = F1G1im + F1G2imB - F2G1imB - F2G2imB ;

59 60 61

re_amplituda [ x ][ y +100] = F1G1re - F1G2reB + F2G1reB - F2G2reB ; im_amplituda [ x ][ y +100] = F1G1im - F1G2imB + F2G1imB - F2G2imB ;

62 63

64

re_amplituda [ x +100][ y +100] = F1G1re - F1G2reB - F2G1reB + F2G2reB ; im_amplituda [ x +100][ y +100] = F1G1im - F1G2imB - F2G1imB + F2G2imB ;

65 66

67

68

69

intenzita [ x ][ y ] = ( Math . pow ( re_amplituda [ x ][ y ] ,2) + Math . pow ( im_amplituda [ x ][ y ] ,2) ) /256; intenzita [ x +100][ y ] = ( Math . pow ( re_amplituda [ x +100][ y ] ,2) + Math . pow ( im_amplituda [ x +100][ y ] ,2) ) /256; intenzita [ x ][ y +100] = ( Math . pow ( re_amplituda [ x ][ y +100] ,2) + Math . pow ( im_amplituda [ x ][ y +100] ,2) ) /256; intenzita [ x +100][ y +100] = ( Math . pow ( re_amplituda [ x +100][ y +100] ,2) + Math . pow ( im_amplituda [ x +100][ y +100] ,2) ) /256;

70 71

72 73 74 75 76 77 78

intenzita_prum = intenzita_prum + intenzita [ x ][ y ] + intenzita [ x +100][ y ] + intenzita [ x ][ y +100] + intenzita [ x +100][ y +100]; i f ( x == y ) { g . setColor ( Color . black ) ; g . fillRect (275 ,75 ,325 ,125) ; g . setColor ( Color . white ) ; g . drawString ( String . valueOf ( x ) + " % " ,295 ,95) ; }

79 80 81

} } Na řádcích 1 a 3 se objevují cykly, které zajistí, že hodnoty proměnných x a y proběhnou od 0 do 99. Pro každou dvojici hodnot x a y je v souladu se vztahem (17) proveden výpočet F1G1, F1G2, F2G1 a F2G2, což je zajištěno cykly na řádcích 13 a 15, ve kterých běží indexy j a k od hodnoty 0 do hodnoty 24. Všimněme si, že až na rozdílné indexy ve funkci f (j, k) je obsah sum ve vztazích (17) totožný. Vždy se totiž objevuje stejný fázový člen eπi(xj+yk)/2M . Z tohoto důvodu je výhodné nejprve podmínkou na řádku 17 ověřit, zda bude v daném proběhnutí cyklů potřeba fázový člen počítat a pokud ano, je na řádcích 18 a 19 vypočtena jeho reálná resp. imaginární složka a hodnoty jsou uloženy do proměnných 9


11. června 2006

FxGxre resp. FxGxim. Na řádcích 22 – 40 vidíme, jak je fázový člen přičten k proměnným F1G1, F1G2, F2G1 nebo F2G2, jsou-li otevřeny na apertuře otvory o příslušných souřadnicích. Cykly na řádcích 13 a 15 představují pouze sumy ze vztahů (17). Vidíme, že je však ještě třeba započítat korekční fázi, která se v těchto vztazích vyskytuje před znamením sumy. Násobení korekční fází v algoritmu zajišťují řádky 45 – 52, kde je třeba provést násobení komplexního čísla komplexním číslem podle předpisu ω1 = (a1 + ib1 ) ∈ C, ω2 = (a2 + ib2 ) ∈ C =⇒ ω1 ω2 = (a1 a2 − b1 b1 ) + i(a1 b2 + a2 b1 ).

(19)

Následně jsou na řádcích 54 – 64 vypočteny reálné a imaginární složky amplitud podle vztahů (18). Z těchto amplitud poté na řádcích 66 – 69 algoritmus vypočte intenzity bodů stínítka. Řádek 71 slouží k výpočtu průměrné intenzity užité dále při škálování. Na řádcích 72 – 78 je programován krátký kód, který zajišťuje zobrazení stavu výpočtu podobně, jak je tomu v případě základního algoritmu. Optimalizovaný algoritmus představuje sice významné urychlení výpočtů, ale přináší s sebou i některá omezení. Klíčová pro odvození algoritmu je periodicita funkce eiφ . Pokud bychom umožnili uživateli měnit vlnovou délku, která se v případě základního algoritmu objevuje právě v exponenciální funkci, zavedli bychom do našich výpočtů neodstranitelnou aperiodicitu, při které celé výše zmíněné odvození optimalizace selže. Stejné problémy se ztrátou výhodné periodicity by se objevily v případě, kdy by uživatel mohl do výpočtu zasáhnout prostřednictvím lupy tak, jak má tuto možnost v základním algoritmu. Tato omezení jsou důvodem, proč má uživatel k dispozici volbu mezi prvním a druhým algoritmem. Nemá-li zvláštní požadavky na lupu nebo na vlnovou délku, může využít rychlejšího optimalizovaného algoritmu. V případě tohoto algoritmu je vynechání úvah o vlnové délce znázorněno použitím bílé barvy při vykreslení difrakčního obrazce. Studium závislosti difrakce na vlnové délce a použití lupy je k dispozici u pomalejšího základního algoritmu. 2.2.4

Vykreslení difrakčního obrazce

Po vypočtení intenzit ve všech bodech může začít vykreslení difrakčního obrazce (zdrojový kód 3). Zdrojový kód 3: Vykreslení difrakčního obrazce 1 2

( i n t x = 0; x < 200; x = x + velikostBodu ) { f o r ( i n t y = 0; y < 200; y = y + velikostBodu ) {

for

3 4 5 6

7

intenzitaBodu = intenzita [ x ][ y ]/ intenzita_prum * zdroj /2; i f ( intenzitaBodu > 1) { intenzitaBodu = 1;} g . setColor ( n e w Color (( i n t ) Math . round ( intenzitaBodu * RGB [0]) ,( i n t ) Math . round ( intenzitaBodu * RGB [1]) ,( i n t ) Math . round ( intenzitaBodu * RGB [2]) ) ) ; g . fillRect ( x +201 , y , velikostBodu , velikostBodu ) ;

8 9 10

} } Vykreslení je již nezávislé na volbě použitého algoritmu a je zajištěno cykly na řádcích 1 a 2, které proběhnou díky indexům x a y celé stínítko. Na řádku 4 je nejprve pro každý bod přeškálována hodnota intenzity pomocí hodnoty průměrné intenzity intenzita_prum a pomocí proměnné zdroj, jejíž hodnotu nastavuje přímo uživatel a která slouží k drobným korekcím jasu. Po tomto výpočtu je na řádku 6 provedeno nastavení barvy, která bude použita. Barvu ovlivňuje hodnota intenzity i hodnota pole RGB[] určená vlnovou délkou. V případě, kdy uživatel vybral výpočet prostřednictvím optimalizovaného algoritmu, je barva nastavena na RGB[i] = 255, pro i = 0, 1, 2. Posledním krokem je vykreslení plošky na stínítku funkcí fill.Rect() na řádku 7.

10


2.2.5

11. června 2006

Porovnání časové náročnosti představených algoritmů

Nyní uzavřeme část věnovanou algoritmickému řešení apletu porovnáním časové náročnosti výpočtů pomocí základního a optimalizovaného algoritmu. Za tímto účelem byl proveden počítačový experiment, při kterém byl z otvorů na apertuře sestaven čtverec o straně délky n (viz obrázek 3). Pro různé

n

n Apertura A

Obrázek 3: Schéma apertury pro účely měření rychlosti algoritmů. hodnoty n byly měřeny časy potřebné k výpočtu pomocí základního algoritmu (tz ) a pomocí optimalizovaného algoritmu (to ) při maximálním rozlišení stínítka. Naměřené hodnoty jsou sestaveny v tabulce 1. V této tabulce je také vypočítán koeficient zrychlení výpočtu k definovaný předpisem k=

n tz [ms] to [ms] k

1 1 488 524 2,84

5 2 469 649 3,80

10 5 838 827 7,06

15 11 483 1 236 9,29

20 19 373 1 637 11,83

tz . to

25 29 639 2 344 12,64

(20)

30 42 497 2 933 14,49

35 56 924 3 969 14,34

40 74 106 4 760 15,57

45 94 186 6 142 15,33

50 115 329 7 136 16,16

Tabulka 1: Časy potřebné k výpočtu n2 prvků s využitím obou algoritmů. Průběh závislosti tz (n) a to (n) je znázorněn na obrázku 4. Naměřené hodnoty je třeba chápat pouze jako orientační, protože měření probíhalo na osobním počítači4 za běžného chodu, kdy je kromě apletu spuštěno mnoho různých procesů běžících na pozadí. Tyto procesy mohly zkreslit výsledek, avšak z obrázku 4 je jasně patrný trend, kdy s rostoucím počtem prvků n dosahuje optimalizovaný algoritmus přibližně 16× kratšího výpočetního času. To je možné vysvětlit tím, že v případě optimalizovaného algoritmu je potřeba, aby každý ze čtyř použitých indexů (x, y, j, k) proběhl v porovnání se základním algoritmem pouze polovinu hodnot. Tím, že se jedná o vnořené cykly, obdržíme pro koeficient urychlení k hodnotu přibližně 24 = 16. Zdrojový kód tohoto simulačního programu obsahuje necelých 850 řádků.

2.3

Návod k použití apletu

Uživatelské rozhraní simulačního apletu je znázorněno na obrázku 5. Okno lze rozdělit na části A, B a C. V části A se nachází plocha, do níž uživatel myší zakresluje otvory v apertuře. Plocha je rozdělena na 50×50 otvorů a každý otvor může být buď otevřený, nebo uzavřený. Program znázorňuje otevřené otvory bílou barvou a otvory uzavřené barvou černou. Otevření otvoru provede uživatel umístěním 4

Parametry použitého počítače: procesor AMD Thunderbird 1GHz, operační systém Linux Mandriva 2006, úroveň běhu 5, verze Javy 1.5.0 06.

11


11. června 2006

12 11

to

10

tz

9 t × 104 [ms]

8 7 6 5 4 3 2 1 0

0

10

20

n

30

40

50

Obrázek 4: Závislost časů tz a to na hodnotě n.

myši nad příslušný bod apertury a stisknutím levého tlačítka. Obdobně lze stisknutím pravého tlačítka otvor uzavřít. Podrží-li uživatel levé tlačítko myši stisknuté, jsou otevřeny všechny otvory přes které myší přejede. Stejným způsobem lze při podržení pravého tlačítka otvory uzavírat. Část B poskytuje uživateli ovládací prvky, kterými lze nastavovat parametry simulace. K dispozici je zde „Lupaÿ, která slouží k detailnějšímu pohledu na prostřední část difrakčního obrazce. Lupa nabízí prostřednictvím posuvníku 4 různé úrovně zvětšení: 1×, 2×, 3× a 4×. Další ovládací prvek označený jako „Vlnová délkaÿ nastavuje použitou vlnovou délku. Posuvníkem si uživatel volí hodnoty v rozsahu 390–760 nm. Ovládací prvek umožňuje sledovat, jak se bude měnit nejen barva světla, ale především difrakční obrazec v závislosti na použité vlnové délce. Aktuálně nastavená vlnová délka se zobrazuje pod difrakčním obrazcem. Po spuštění apletu je jako výchozí nastavena hodnota 555 nm. Posuvník označený „Intenzita zdrojeÿ jemně dolaďuje √ jas difrakčního obrazce. Měřítko je logaritmické a jeden skok představuje násobení hodnotou 3 10. Program automaticky počítá vhodnou hodnotu jasu a tu ztotožní s prostřední polohou posuvníku. Uživatel získává možnost jemně zvýšit i snížit jas podle aktuálních podmínek a potřeb (např. nastavení monitoru). Další v pořadí následuje posuvník „Rozlišeníÿ, kterým si uživatel volí rozlišení difrakčního obrazce. Levá krajní pozice posuvníku odpovídá nejnižšímu rozlišení a pravá krajní pozice odpovídá rozlišení nejvyššímu. Přitom čím vyšší je rozlišení, tím delší dobu trvá difrakční obrazec vykreslit. Výchozí nastavení rozlišení po startu programu je druhé nejvyšší. V literatuře je často sledována difrakce na několika typických otvorech. Jedná se například o difrakci na kruhovém otvoru nebo na otvoru obdélníkového tvaru. Aby měl uživatel možnost rychle simulovat tyto typické případy, je v programu sada významných otvorů předem připravena. K výběru z této sady slouží výběrový seznam nazvaný „Štěrbinaÿ. Nabídka otvorů obsahuje obdélníkový otvor, kruhový otvor, dva obdélníkové otvory, dva kruhové otvory a otvor tvaru kružnice. Dále je k dispozici série bodových otvorů pro studium interference v provedení postupně 2, 3, 5, 7, 11 a 13 bodů rozmístěných po obvodu kružnice a následují ještě otvory tvaru hvězdice, obdélníku v rámu a rámečku. Další volba s názvem „Algoritmusÿ umožňuje vybrat mezi dvěma algoritmy implementovanými v apletu – základní a optimalizovaný. Základní algoritmus je pomalejším algoritmem, ale umožňuje

12


11. června 2006

Obrázek 5: Vzhled apletu pro simulaci difrakčních jevů

měnit vlnovou délku a použít lupu. Oproti tomu optimalizovaný algoritmus je, co se výpočtu týče, rychlejší, ale volba vlnové délky a lupy není při výběru tohoto algoritmu dostupná. V části B jsou ještě dvě tlačítka. První s nápisem „Vyčisti všeÿ zavře všechny otevřené otvory na apertuře a druhé s nápisem „Překresliÿ generuje difrakční obrazec na základě navolených parametrů. Poslední část C představuje stínítko, kde může uživatel difrakční obrazec pozorovat. Během výpočtu se na stínítku zobrazuje stav v procentech a po dokončení se zde vykreslí difrakční obrazec. Závěrem si představíme na obrázku 6 několik ukázek difrakčních obrazců, které je možné pomocí předkládaného apletu vytvořit. Simulací získané obrazce se dobře shodují s experimentálními daty [15, 16].

a)

b)

c)

d)

Obrázek 6: Ukázky několika difrakčních obrazců vytvořených prezentovaným apletem: a) difrakce na kruhovém otvoru, b) difrakce na obdélníkovém otvoru, c) Youngův pokus - interference dvou bodových zdrojů, d) interference 11 bodových zdrojů rozložených po obvodu kružnice.

13


11. června 2006

Použitá literatura [1] Diffraction: http://fr.wikipedia.org/wiki/Diffraction [2] BORN, M.; WOLF, E.: Principles of Optics 6th edition, Pergamen Press, Oxford 1980, p. 370–386. [3] Théorie de la diffraction: http://fr.wikipedia.org/wiki/Théorie de la diffraction [4] Huygens Principle: http://www.mathpages.com/home/kmath242/kmath242.htm [5] A brief history of Optics: http://members.aol.com/WSRNet/D1/hist.htm [6] Encyklopedie Coto.je: http://www.cotoje.cz [7] SALEH, B. E. A.; TEICH, M. C.: Fundamentals of photonics, A Wiley-Interscience publication John Wiley sons Inc., NewYork 1991, p. 127–130. [8] HENRY, M.: Optique ondulatoire, Serveur Techniques de l’ingénieur, Dossier A 191, p. 10–15 http://www.techniques-ingenieur.fr/dossier/optique ondulatoire interferences diffraction polarisation/A191

[9] Fraunhofer Diffraction and 2D Fourier Transform: http://oldsite.vislab.usyd.edu.au/CP3/Four6/node1.html [10] WEISSTEIN, E. W.: Discrete Fourier Transform, From MathWorld – A Wolfram Web Resource. http://mathworld.wolfram.com/DiscreteFourierTransform.html [11] TVERSKOJ, P. N.: Optické, elektrické a akustické jevy v atmosféře., Naše vojsko, Praha 1955 in http://www.pef.zcu.cz/pef/kof/cz/st/dp/hosnedl/html/Slunecni.htm [12] COOLEY, J. W.; TUKEY, J. W.; An Algorithm for the Machine Calculation of Complex Fourier Series, Math. Comp. 19 (1965), p. 297–301 [13] WEISSTEIN, E. W.: Fast Fourier Transform, From MathWorld – A Wolfram Web Resource. http://mathworld.wolfram.com/FastFourierTransform.html [14] MARSHALL, D.: The Fast Fourier Transform Algorithm http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL COPIES/MARSHALL/node20.html [15] LIPSON, S. G.; LIPSON H.; TANNHAUSER D. S.: Optical Physics 3rd edition, Cambridge Univ. Press, Cambridge 1995, p. 185, 186 [16] TOMASINO, A. et al.: Physique Term. S, Nathan, Paříž 1995, p. 300–304

14


3 3.1

11. června 2006

Počítačová simulace šíření světelného paprsku Teoretický základ

V předchozí části jsme viděli, jak prostřednictvím difrakce demonstrujeme vlnovou povahu světla. Často je možné zanedbat vlnový charakter světla a omezit se na jednodušší představu světelného paprsku. V této souvislosti hovoříme o geometrické optice, kdy předpokládáme, že energie přenášená světlem se šíří po přesně definovaných drahách. Experimenty ukazují, že tvar těchto paprsků bude záviset na vlastnostech prostředí, kterým se světlo šíří. Veličina spojená s vlastnostmi prostředí ovlivňující šíření světla se nazývá index lomu a definuje se pomocí vztahu n=

c , v

(1)

ve kterém v označuje fázovou rychlost šíření světla v daném místě prostředí a c rychlost šíření světla ve vakuu [1, 2, 3, 4]. V homogenním prostředí (n = konst) se světlo šíří po přímce. Při odrazu světla na odrazné ploše (viz obrázek 1a) platí zákon odrazu ve tvaru α1 = α10 ,

(2)

kde α1 označuje úhel dopadu a α10 úhel odrazu. V případě, kdy světlo dopadne na rozhraní dvou prostředí s hodnotami indexů lomu n1 v prvním prostředí a n2 ve druhém prostředí, bude pro úhel dopadu α1 a úhel lomu α2 (viz obrázek 1b) platit Snellův zákon n1 sin α1 = n2 sin α2

(3)

pojmenovaný podle holandského fyzika Willebrorda van Roijen Snella (1580 nebo 1591 – 30. 10. 1626). V případě, kdy dochází k šíření světla z prostředí s vyšším indexem lomu do prostředí s nižším indexem A

B

A

A

B

α1 n1 n2

α1′ α1

α1′ α1 α2

a) Odraz paprsku

B

n1 n2

b) Lom paprsku

´ y odraz c) Upln´

Obrázek 1: Ilustrace k zákonům odrazu a lomu lomu, tedy n1 > n2 , může dojít pro určitou hodnotu úhlu dopadu k lomu paprsku rovnoběžně s rovinou rozhraní. Úhlu dopadu, při kterém k takto popsané situaci dojde, říkáme mezný úhel αm a ze Snellova zákona pro něj platí podmínka n2 sin αm = . (4) n1 Pro úhly dopadu větší než je mezný úhel (α1 > αm ), již nebude docházet k lomu světla do druhého prostředí, ale paprsek se odrazí zpět do prvního prostředí (viz obrázek 1c). Říkáme, že došlo k úplnému odrazu. Jak bude dále ukázáno, simulační program testuje při svém běhu na každém rozhraní, zda dochází k lomu, nebo k úplnému odrazu. Odraz světla i lom světla mají jednu důležitou společnou vlastnost. Jak již mezi svými axiomy postuluje Isaac Newton (4. 1. 1643 – 31. 3. 1727), rovina dopadajícího paprsku je totožná s rovinou paprsku odraženého i lomeného. Při hlubším studiu problému lomu se ukazuje, že kromě lomeného paprsku lze také pozorovat paprsek odražený, jehož směr se řídí 15


11. června 2006

zákonem odrazu. Energie, kterou odražený paprsek přenáší, je však za běžných okolností mnohonásobně menší než energie přenášená lomeným paprskem. Z tohoto důvodu simulační aplet zanedbává odraz v případě, kdy dochází k lomu a nevykresluje odražený paprsek [1, 2, 3, 4, 5]. Simulace si klade za cíl znázornit chod paprsku světla v libovolném uzavřeném mnohoúhelníku. Mnohoúhelník, který si uživatel zvolí, definuje hranici mezi vnějším a vnitřním prostředím. Každé toto prostředí může být charakterizováno jiným indexem lomu. Paprsek je vyslán ze zdroje, jehož polohu i směr může uživatel během simulace upravovat. Paprsek se bude šířit přímočaře, dokud nenarazí na rozhraní mezi vnitřním a vnějším prostředím. V reálném materiálu může navíc docházet k disperzi, tedy k rozkladu bílého světla na spektrální složky. Simulace dává možnost studovat disperzi ve vnitřním prostředí. Uživatel si volí ze tří režimů: vypnutá disperze, reálná disperze a zvýrazněná disperze. V případě, kdy je disperze vypnuta, bude program vykreslovat pouze jeden paprsek bílé barvy, který se bude řídit uživatelem zadaným indexem lomu vnitřního prostředí. Je-li nastaven režim reálné disperze, provede se vykreslení tří paprsků (červený, zelený a modrý). Červený paprsek bude pracovat s hodnotou indexu lomu, kterou pro vnitřní prostředí zadá uživatel. Pro zelený paprsek je použita hodnota o 0,01 vyšší. Pro trasování modrého paprsku je použita hodnota indexu lomu o 0,02 vyšší, než je hodnota nastavená uživatelem. V režimu zvýrazněné disperze je situace obdobná, ale rozdíly mezi použitými indexu lomu jsou 0,05 [2].

3.2

Algoritmické řešení apletu

Jádrem celé simulace je cyklus představený ve zdrojovém kódu 4. Tento cyklus, začínající řádkem 1, provádí vykreslení paprsku tak, jak by se skutečně světlo šířilo. Cyklus běží neustále dokola a v každém svém opakování provede znázornění chodu paprsku mezi dvěma rozhraními vnějšího a vnitřního prostředí. Ukončen je v případě, když už v cestě paprsku neleží žádné další rozhraní, na které by mohl dopadnout, nebo v případě, že je již provedeno 100 běhů cyklu. Toto opatření chrání uživatele před případnou nekonečnou smyčkou. Zdrojový kód 4: Cyklus zajišťující znázornění chodu paprsku 1

while

( cyklus < 100 && (( kresli_zdroj == f a l s e && kresli_rotace == ) || rezimCHB . getState () == t r u e ) ) { cyklus ++; i f (! najdi_rozhrani () ) b r e a k ; / / s p u s t í h l e d á n í r o z h r a n í , n e n a j d e

false 2 3 4

- li

žádné ,

ukončí

běh

cyklu

5 6

g . setColor ( barva_cary ) ;

7 8

if

9

if

( velky_cyklus == 1 && nedisperzni_cyklus [1] >= cyklus ) { g . setColor ( Color . cyan ) ;} ( velky_cyklus == 1 && nedisperzni_cyklus [0] >= cyklus ) { g . setColor ( Color . white ) ;}

10 11

g . drawLine ( paprX , paprY , prusX , prusY ) ; / / paprsku

s

nalezeným

průsečíkem

s

spojí

aktualní

pozici

rozhraním

12 13

snell () ; / / provede

14

if

spustí

funkci ,

výpočet

která

směrnice

souladu po

se

Snellovým

střetu

s

g . setColor ( Color . red ) ; 16

zákonem

rozhraním

( nedisperzni_lom == t r u e && velky_cyklus > 1) { n edi sp er zn i_ cy kl us [3 - velky_cyklus ]++;}

15 16

v

paprsku


11. června 2006

g . fillRect ( paprX -4 , paprY -4 , 8 , 8) ; / / }

17 18

znázorní

bod

dopadu

Při každém svém běhu volá nejprve cyklus pomocnou funkci boolean najdi_rozhrani(). Tato funkce určí, na které rozhraní mezi vnějším a vnitřním prostředím paprsek dopadne. Zároveň určí bod dopadu paprsku na toto rozhraní. Při tomto vyhledávání prochází funkce boolean najdi_rozhrani() všechna rozhraní, hledá průsečík rozhraní a paprsku a poté vybere to s nejbližším průsečíkem. V případě, kdy je rozhraní nalezeno, vrátí funkce hodnotu true a cyklus běží dál. V opačném případě je vrácena hodnota false a cyklus je příkazem na řádku 4 ukončen. Řádky 6–9 zajišťují správný výběr barvy, kterou bude paprsek znázorněn. V případě vypnuté disperze je použita bílá barva. V případě reálné nebo zvýrazněné disperze se postupně provede vykreslení paprsků červenou, zelenou a modrou barvou. Může se ovšem stát, že dojde k úplnému odrazu, kdy se disperze neprojeví. V takovém případě je i nadále použita bílá barva. Další možnou situací je lom červené barvy a úplný odraz zelené a modré. Paprsek odpovídající zelené a modré barvě dohromady je vykreslován azurovou barvou. Po volbě barvy je znázorněn příkazem na řádku 11 chod paprsku mezi předchozí polohou a bodem na vybraném rozhraní, do kterého dopadne. Dále je nezbytné určit chování paprsku na tomto rozhraní. Výpočet nového směru šíření paprsku zajišťuje pomocná funkce snell() volaná na řádku 13. Zmíněná funkce vyhodnotí, zda dojde k lomu, nebo k úplnému odrazu a podle Snellova zákona (případně podle zákona odrazu) určí novou směrnici paprsku. Poslední činností, kterou cyklus provádí, je znázornění bodu dopadu. Na řádku 16 je nastavena červená barva a na řádku 17 je vykreslen červený čtvereček odpovídající bodu, ve kterém paprsek dopadl na rozhraní. Zdrojový kód tohoto simulačního programu obsahuje přibližně 750 řádků.

3.3


Vzhled apletu je představen na obrázku 2. V horní části je vykreslována samotná simulace a v dolní části je připravena sada ovládacích prvků, kterými uživatel nastavuje parametry prováděné simulace. První řádek ovládacích prvků obsahuje textová políčka pro zadání indexu lomu vnějšího a vnitřního

Obrázek 2: Vzhled apletu pro simulaci šíření paprsku optickým prostředím

17


11. června 2006

prostředí. Index lomu může nabývat hodnoty minimálně 1 a desetinná místa jsou oddělena tečkou.

a)

b)

c)

d)

Obrázek 3: Ukázky několika simulací chodu paprsku: a) rozklad světla na hranolu, b) trasování paprsku pentagonálním odražečem, c) ukázka funkce koutového odražeče, d) trasování paprsku tvarem obrácené „Uÿ.

Na druhém řádku v ovládacích prvcích nastavuje uživatel polohu a směr svícení zdroje. Zdroj paprsku je pevně vázán na levý okraj apletu. Uživatel však může volit jeho vertikální polohu pomocí položky „Výška zdrojeÿ. Směr, kterým bude paprsek ze zdroje vycházet, nastavuje uživatel další volbou „Sklon zdrojeÿ. Údaj je ve stupních. Hodnota 0◦ odpovídá paprsku vyslanému horizontálně a hodnota 90◦ paprsku vyslanému ze zdroje vertikálně směrem dolů. Kromě ovládacích prvků z druhého řádku lze k nastavení zdroje použít i myš. Zdroj je v apletu znázorněn jako zelené kolečko. Jeho tažením pomocí myši může uživatel měnit polohu zdroje. Sklon zdroje je znázorněn polohou zeleného čtverečku. Pomocí myši může uživatel čtverečkem hýbat a měnit tak sklon zdroje. Na třetím řádku ovládacích prvků je prováděno nastavení tvaru mnohoúhelníku, kterým bude paprsek trasován. Položka „Tvarÿ obsahuje nabídku z několika často studovaných mnohoúhelníků. Jsou to hranol, optické vlákno (zjednodušeno jako podlouhlý obdélník), koutový odražeč, čtverec, šestiúhelník a dvanáctiúhelník. Tvar optické vlákno lze také použít jako planparalelní destičku. Chce-li uživatel studovat jiný mnohoúhelník, vybere volbu „uživatelskéÿ. Poté do textového políčka označeného jako „Uživatelské rozhraníÿ zapíše souřadnice vrcholů jím požadovaného mnohoúhelníku. Syntaxe pro tento zápis je následující: každý vrchol má x-ovou a y-ovou souřadnici. Počátkem souřadnic je levý horní roh apletu. X-ová souřadnice roste směrem doprava a y-ová směrem dolů. Mezi x-ovou a y-ovou souřadnici vloží uživatel dvojtečku. Jednotlivé vrcholy pak odděluje uživatel středníkem a středník je také třeba zapsat nakonec za souřadnice posledního vrcholu. Na ukázku uveďme, jak může uživatel vytvořit například rovnoběžník: 100:100;400:100;350:200;50:200;. Poslední řádek ovládacích prvků obsahuje volbu „Úroveň disperzeÿ, kde si uživatel volí ze tří

18


11. června 2006

možností: bez disperze, reálná disperze a zvýrazněná disperze. Význam jednotlivých voleb je popsán v teoretické části. Další volbou je zaškrtávací políčko „Režim kresleníÿ. Je-li vybráno, bude aplet provádět vykreslení paprsku současně i během změny polohy a sklonu zdroje myší. Protože může někdy být takovéto vykreslování v reálném čase náročné na výpočetní výkon, je uživateli umožněno jej vypnout a vykreslení poté proběhne, až uživatel dokončí nastavení zdroje myší (pustí tlačítko). Poslední položkou je tlačítko „Překresliÿ provádějící překreslení cesty paprsku, je-li na něj kliknuto myší. Na závěr představme několik ukázek trasování paprsku tak, jak je simulační aplet znázorňuje (viz obrázek 3).

Použitá literatura [1] BORN, M.; WOLF, E.: Principles of Optics 6th edition, Pergamen Press, Oxford 1980, p. 109–132. [2] SALEH, B. E. A.; TEICH, M. C.: Fundamentals of Photonics, A Wiley-Interscience publication John Wiley sons Inc., NewYork 1991, p. 1–12. [3] GOODMAN, D. S.: Geometrical Optics (in Handbook of Optics 2nd edition), McGraw-Hill Professional, New York 1994, p. 7–29. [4] Encyklopedie Coto.je: http://www.cotoje.cz [5] NEWTON, I.: Opticks: or, a Treatise of the Reflexions, Refractions, Inflexions and Colours of Light. . ., Sam Smith and Benj. Waleord, Londýn 1704, p. 4, 13, 53.

19


4

11. června 2006

Počítačová simulace Lissajousových křivek

4.1

Teoretický základ

Uvažujme bod, který současně vykonává harmonické kmity ve dvou na sebe kolmých směrech. Zajímavé bude studovat jeho trajektorii za předpokladu, že frekvence těchto kolmých kmitů jsou obecně různé a je mezi nimi libovolný fázový posun. Tímto problémem se v průběhu 19. století zabývali zejména dva vědci. Prvním z nich byl již v roce 1815 Nathaniel Bowditch (26. 3. 1773 – 16. 3. 1838). Americký matematik a fyzik, známý též jako autor námořních map, byl první, kdo úlohu řešil. Využíval při tom speciální kyvadlo tvaru Y5 schopné zaznamenávat trajektorii kmitajícího bodu. Na podobném principu pracující přístroje, tzv. harmonografy, se v průběhu 19. století hojně využívaly při fyzikálních experimentech, ale také jako zdroje esteticky zajímavých obrazců. Jejich vývojem a výrobou se komerčně zabývalo i několik firem. Harmonografy využívaly k zápisu křivek nejen různých psacích potřeb, ale také například jemný písek, který se v průběhu experimentu odsypával z kyvadla na podložku. Tato metoda umožňovala z tloušťky pískové stopy odhadnout navíc i rychlost pohybu. Nezávisle na Bowditchovi se studiem dvou kolmých harmonických kmitů zabýval francouzský fyzik Jules-Antoinne Lissajous (4. 3. 1822 – 24. 6. 1880). Ten svou práci publikoval sice až v roce 1857, ale jeho hlavním přínosem bylo využití světla při zobrazení složených kmitů. Experiment, při němž nechal svazek světla postupně odrazit na dvou zrcadlech a posléze dorazit na stínítko, lze považovat za první realizaci osciloskopu. Pokud byla zrcadla připevněna ke dvěma ladičkám, mohl být touto technikou studován i zvuk. Veřejnost měla možnost Lissajousovy experimenty zhlédnout roku 1867, kdy je představil na Světové výstavě v Paříži [1, 2, 3, 4]. V následující části se zaměříme na analytické vyjádření složených kmitů. Ztotožníme-li kmity v jednom směru se směrem osy x a kmity ve druhém směru se směrem osy y, lze polohu kmitajícího bodu zapsat pomocí rovnic x = A1 sin(ω1 t + φ1 )

(1)

y = A2 sin(ω2 t + φ2 ).

(2)

V předchozích vztazích označujeme písmenem A1 amplitudu kmitů ve směru osy x a A2 amplitudu kmitů ve směru osy y. Dále ω1 označuje úhlovou frekvenci kmitů v ose x a ω2 úhlovou frekvenci pro kmity v ose y. Protože mohou být kmity obecně posunuty o libovolnou fázi, je také nutné ji do úvah zahrnout. Pro kmity v ose x je počáteční fáze označena φ1 a pro kmity v ose y φ2 . Proměnná t v rovnicích označuje čas [2, 5]. Z matematického pohledu lze na rovnice (1) a (2) nahlížet jako na parametrické vyjádření trajektorie kmitajícího bodu, kde parametrem je čas. Speciální třídu křivek budou tvořit takové křivky, kdy platí, že ω1 ∈ Q, (3) ω2 tedy poměr frekvencí obou kolmých kmitů je racionální číslo. V tomto případě bude trajektorií bodu uzavřená křivka, kterou nazveme Lissajousovou křivkou6 . Příklady několika vybraných Lissajousových křivek jsou na obrázku 1. Trajektorie bodu, vytvářející Lissajousovu křivku, se uzavře po uplynutí periody složených kmitů. Ta se rovná nejmenšímu společnému násobku period obou skládajících se kmitů. Je vhodné zdůraznit, že tvar Lissajousovy křivky závisí pouze na vzájemném poměru frekvencí skládajících se kmitů. Dodejme, že v případě, kdy nebude splněna podmínka (3), projde postupně kmitající bod všemi body ležícími v obdélníku h−A1 ; A1 i × h−A2 ; A2 i [2]. Zajímavou praktickou aplikací Lissajousových křivek je zpětné odečtení poměru frekvencí generujících kmitů. Jak je znázorněno na obrázku 2, Lissajousovu křivku protneme vertikální přímkou p 5

Toto kyvadlo bývá v literatuře zmiňováno jako Blackburnovo kyvadlo. Prvním doloženým konstruktérem takovýchto kyvadel byl však James Dean. 6 Některé zdroje hovoří vzhledem k Bowditchovu prvenství o tzv. Bowditchových křivkách [6]

20


11. června 2006

Obrázek 1: Příklady Lissajousových křivek pro poměry frekvencí ω1 : ω2 a) 2:1 b) 3:2 c) 5:3. Tyto příklady byly pořízeny simulačním apletem prezentovaným v této práci.

(kolmou na osu x) a horizontální přímkou q (kolmou na osu y). Přímky volíme tak, aby počet jejich průsečíků s Lissajousovou křivkou byl maximální. Ukazuje se, že poměr počtu průsečíků Lissajousovy křivky s přímkou p a q je roven poměru mezi frekvencemi ω1 : ω2 .

Obrázek 2: Grafické znázornění odečtu poměru frekvencí z Lissajousovy křivky.

Uvedené metody lze využít pro určení neznámé frekvence kmitů. Stačí neznámé kmity složit kolmo s kmity o známé frekvenci a výše zmíněná metoda nám umožní určit hledanou frekvenci [2]. Lissajousovy křivky mají v dnešní době široké spektrum uplatnění. Bylo například ukázáno, jak lze pomocí Lissajousových křivek zpřesnit stanovení koncentrace kyslíku v krvi [7]. Praktické aplikace Lissajousových křivek jsou hlavní motivací k jejich studiu a k simulaci například prostřednictvím Java apletů. Lissajousovy křivky lze zobecnit i na trojrozměrný případ. Této problematice jsou například věnována sdělení [8, 9].

4.2

Algoritmické řešení apletu

Programovací jazyk Java není schopen algebraických úprav a řešení rovnic. Grafická primitiva7 tohoto jazyka zahrnují pouze základní geometrické útvary. Simulaci trajektorie je třeba provést vykreslováním bod po bodu, jako bychom skutečně sledovali kmitající těleso. Při vykreslování je použito následujícího cyklu, který při každém proběhnutí vykreslí jeden bod křivky. Zdrojový kód 5: Algoritmus pro vykreslení Lissajousových křivek 1 2

while

( t < stop_t +1) { pomer = ( d o u b l e ) t / stop_t ;

7 Grafickými primitivy se rozumí soubor implementovaných grafických prvků, které je možné příslušnými příkazy vykreslovat [10].

21


11. června 2006

X = ( i n t ) Math . round ( meritko *100* ampX * Math . sin (( pomer * Math . PI *2) * freqX + fazerad ) +120) - posunX ; Y = ( i n t ) Math . round ( meritko *100* ampY * Math . sin (( pomer * Math . PI *2) * freqY ) *( -1) +120) - posunY ;

3

4

5

i f ( X <= 240 && Y <= 240) { g . drawRect (X , Y , 1 , 1) ;} t ++; }

6 7 8

Proměnná t označující čas postupně mění své hodnoty v intervalu od nuly až do maximální hodnoty stop_t, kterou uživatel může měnit ovládacím prvkem označeným jako „Citlivostÿ. Při takto realizované technice vykreslovaní křivek platí, že čím více bodů je na vykreslení použito, tím větší přesnosti lze dosáhnout. Proto je také volba, určující kolikrát cyklus proběhne a tedy kolik bodů se vykreslí, označena jako „Citlivostÿ. Při každém opakování cyklu je vypočtena hodnota nové proměnné pomer, která je definována jako podíl t a stop_t, a následně je také hodnota t zvětšena o jedničku (t++). Dále je třeba podle vztahů (1) a (2) spočítat souřadnice X, Y kmitajícího bodu. Tyto souřadnice jsou v Javě vyjádřeny pomocí proměnných celočíselného typu (integer), avšak výsledkem operace sinus, kterou v Javě zajišťuje funkce Math.sin, je číslo reálné (typ double). Je tedy nezbytné výsledek funkce sinus zaokrouhlit a přetypovat jej tak na celočíselnou hodnotu. Při tomto kroku dochází nevyhnutelně ke zkreslení výpočtu, protože ztrácíme veškerou přesnost obsaženou v desetinných místech. V obou osách se může aplet při zaokrouhlování dopustit maximální chyby 0,5 px. Prostor pro vykreslení křivky má rozměry 200×200 px. Z celkového rozsahu os tak může být zaokrouhlením způsobena chyba maximálně 0,25 %. Ukazuje se, že tuto nepřesnost lze zanedbat. Po provedení výpočtu souřadnice je programem testována poloha bodu, který má být vykreslen. Cílem je ověřit, zda bod leží v oblasti apletu určené pro kreslení křivky a ne jinde (třeba v oblasti ovládacích prvků). Splňuje-li bod tuto podmínku, je vykreslen jako malý čtvereček o rozměrech 1×1 px. K tomu je využita funkce drawRect. Abychom se ujistili, že naprogramovaný cyklus skutečně projde celou časovou periodu složeného kmitu, je v argumentu funkce sinus použit součin uživatelem zvolené frekvence, proměnné pomer a konstanty 2π. Je-li přitom uživateli umožněno za frekvenci volit minimálně hodnotu jedna, máme jistotu, že cyklus skutečně alespoň jedenkrát celou periodu projde. To, že nemá uživatel možnost volit frekvence nižší než jedna, nikterak nesnižuje obecnost simulace, protože, jak bylo zmíněno v teoretickém rozboru, výsledný tvar Lissajousovy křivky záleží pouze na vzájemném poměru frekvencí skládajících se kmitů. Jako příklad uveďme, že volbu frekvencí 1 a 1/3 lze úspěšně nahradit volbou frekvencí 3 a 1, aniž by došlo ke změně tvaru Lissajousovy křivky. Ostatní proměnné, které se ve výše uvedeném zdrojovém kódu vyskytují, představují různé matematické korekce související pouze s obrazovým výstupem. Uživatel bývá zvyklý na počátek souřadného systému uprostřed zobrazovací plochy a na osu y, jejíž hodnoty rostou směrem nahoru. Vzhledem k tomu, že v Javě je počátek souřadného systému umístěn do levého horního rohu a osa y roste směrem dolů, je vhodné tento rozdíl korigovat pro zvýšení uživatelova pohodlí. Celý zdrojový kód tohoto apletu zahrnuje téměř 300 řádků. Na většině z nich je programováno a spravováno uživatelské rozhraní, které s apletem umožní snadno pracovat. Aplet využívá obvyklé ovládací prvky, jako jsou textová políčka a posuvníky.

4.3


Vzhled apletu je znázorněn na obrázku 3. Ovládání je zcela intuitivní. Uživatel má k dispozici sadu interaktivních prvků (textová pole a posuvníky), pomocí kterých může nastavit parametry simulace. Každý parametr simulace lze zadat z klávesnice prostřednictvím textového políčka nebo jej lze nastavit myší prostřednictvím příslušného posuvníku. Základními parametry simulace jsou amplitudy

22


11. června 2006

Obrázek 3: Vzhled apletu pro simulaci Lissajousových křivek

kmitů v osách x a y. Amplitudy nastavuje uživatel prvními dvěma posuvníky nebo textovými poli a mohou nabývat hodnot od 1 do 10. Zde má také uživatel možnost ověřit, že výsledný tvar Lissajousovy křivky nezávisí na jednotlivých hodnotách amplitud, ale až na jejich poměru. Další dvojice ovládacích prvků určuje frekvence kmitů v osách x a y. Uživatel má možnost volit frekvence v rozsahu 1 až 20. Tyto volby umožní, mimo jiné, pozorovat, že pro tvar Lissajousovy křivky je rozhodující pouze poměr mezi těmito frekvencemi. Následující prvek je označený jako „Fázový posunÿ a dovoluje uživateli přidat do kmitů v ose x libovolný fázový posun v rozsahu -90◦ až +90◦ . Tak může uživatel simulovat závislost tvaru křivky na fázovém posunu. Poslední prvek prvního sloupce ovládá již zmíněnou citlivost, a to tak, že udává počet bodů, pomocí kterých bude křivka vykreslena. Zde však znamená více bodů také delší výpočetní čas. Druhý sloupec obsahuje volby, které nesouvisejí přímo s tvarem Lissajousovy křivky, ale s jejím vykreslením na monitoru. První položka, nazvaná „Lupaÿ, umožňuje zvětšovat nebo zmenšovat velikost Lissajousova obrazce. Posuvník, kterým se lupa nastavuje, využívá logaritmické měřítko. Další dvě volby („Posun Xÿ a „Posun Yÿ) slouží k posunutí uživatelova zorného pole horizontálním a vertikálním směrem.

23


11. června 2006

Použitá literatura [1] MAOR, E.: Trigonometric Delights, Princeton University Press, Princeton 1998, p. 145-149 [2] BAJER, J.: Mechanika 3, VUP, Olomouc 2006 (v tisku), p. 256-262. [3] http://www.mathcurve.com/courbes2d/lissajous/lissajous.shtml [4] WHITAKER, R. J.: Harmonographs. I. Pendulum design. Am. J. Phys. 69 (2001), p. 162-173. [5] GOLDSTEIN, H. et al.: Classical Mechanics (3rd edition). Addison Wesley, San Francisco 2002, p. 437-440. [6] http://poncelet.math.nthu.edu.tw/chuan/zoo/toc.html [7] ADDISON, P. S. and WATSON, J. N.: A novel time-frequency-based 3D Lissajous figure method and its application to the determination of oxygen saturation from the photoplethysmogram. Meas. Sci. Technol. 15 (2004), p. 15-18. [8] BOGLE, M. G. V.; HEARST, J. E.; JONES, V. F. R. and STOILOV, L.: Lissajous knots. Journal of Knot Theory and Its Ramifications 3(2) (1994), p. 121-140. [9] http://www.mathcurve.com/courbes3d/lissajous3d/lissajous3d.shtml [10] HEROUT, P.: JAVA - grafické uživatelské prostředí a čeština, Kopp, České Budějovice 2001.

24


5

11. června 2006

Závěr

V této práci jsem představil tři simulační aplety, které jsem vytvořil. Při jejich tvorbě se ukázalo, že vždy existuje prostor pro rozšíření stávajících simulací o další zajímavé možnosti. Také se ovšem ukázalo, že vývoj simulačních programů je časově poměrně náročný a to zejména ve fázi finálního „laděníÿ programů. Při programování se projevily některé přednosti (snadná obsluha uživatelských akcí), ale také nevýhody (absence předdefinovaných objektů jako jsou např. vektory) programovacího jazyka Java. Věřím ovšem že přednosti převažovaly. Do budoucna se programování apletů plánuji věnovat alespoň okrajově. Již nyní (duben 2006) je ve fázi vývoje další aplet, který studuje statickou rovnováhu uživatelem postaveného bloku cihel.

25

Počítačové simulace vybraných fyzikálních jevů

Recommend Documents