MASARYK UNIVERSITY F9602 Interakce elektromagnetického pole se živou hmotou
GSVIT hands-on tutorial
Petr Steindl
Faculty of Science Department Condensed matter physics
1. září 2017
Obsah
1 Teorie elektromagnetického pole 2
1.1 Maxwellovy rovnice............................................. 2
1.2 Optické vlastnosti látek.......................................... 3
2 FDTD solver GSVIT 5
2.1 Metoda FDTD............................................... 5
2.1.1   Volba velikosti kroku........................................ 6
2.2 GSVIT.................................................... 7
2.3 Jak seznam parametrů vytvořit...................................... 7
3 COMPUTATIONAL DOMAIN 9
4 BASIC PARAMETERS 10
5 BOUNDARY CONDITIONS 12
5.1   type..................................................... 12
6 SOURCE 14
7 MEDIA 17
7.1 Vytvoření souboru pro zkoumaný předmět................................ 17
7.1.1 code_shape ............................................. 18
7.1.2 shape_volume_def.......................................... 18
7.1.3 Načtení ze souboru......................................... 19
7.1.4 material............................................... 20
7.2 MATMODE_CHECK ........................................... 21
8 OUTPUTS 22
9 NEAR FIELD to FAR FIELD 25
10 Příklad: Odraz na rovinném rozhraní 27
10.1 Zadání problému.............................................. 27
10.2 Parametrový soubor............................................ 27
10.3 Výsledky simulace............................................. 30
Literatura 33
Projevy interakce světla s materiálovými předměty pozorujeme v každodenním životě a k jejich popisu byla vybudována komplexní teorie interakce elektromagnetického záření v materiálovém prostředí Maxwellem [1]. Pro řešení interakce světla při dopadu na rovinné rozhraní či pro interakci s kulovým předmětem existují analytická řešení, pro složitější geometrie je třeba použít numerické řešení.
V posledních letech se oblíbenou metodou na řešení jevů spojených s šířením elektromagnetického záření v materiálovém prostředí se složitou geometrií stala metoda konečných diferencí v časové oblasti (Finite Difference Time Domain, FDTD).
V první kapitole připomeneme základní teoretický popis pomocí Maxwellových rovnic. V kapitole 2 nastíníme fungování FDTD metody a představíme si FDTD solver GSVIT. Kapitoly 3-9 slouží jako příručka ke GSVITu a přibližují jeho používání. V poslední kapitole demonstrujeme řešení příkladu odrazu elektro-magnetického vlnění na rovném rozhraní kov-vzduch.
1
Kapitola 1
Teorie elektromagnetického pole
1.1   Maxwellovy rovnice
Interakci světla s materiálovým objektem lze popsat pomocí Maxwellových rovnic:
V D = pu (1.1) V-B   =   0, (1.2)
VxE   =   --, (1.3) dD
VxH   =   Jf + ^- (1-4)
Tyto makroskopické rovnice vzájemně propojují elektrickou D a magnetickou B indukci s intenzitou elektrického E a magnetického H pole skrze hustotu volného náboje p f a hustotu volného proudu Jf.
Elektrické a magnetické pole při šíření v materiálovém prostředí vyvolává elektrickou (polarizaci P) a magnetickou (magnetizaci M) odezvu, která modifikuje intenzitu celkového elektrického a magnetického pole dle následujících vztahů
D   =   e0E + P, (1.5)
H   =   —B-M, (1.6)
Mo
kde eo a po jsou elektrická permitivita a magnetická permeabilita vakua.
Makroskopické (globálni) Maxwellovy rovnice (1.2)-(1.7) jsou doplněny třemi rovnicemi mikroskopickými (lokálními)
Jf = <rE, (1.7) M = /ioArH, (1.8) P   =   e0xE, (1.9)
kde a je tenzor elektrické vodivosti, pr je tenzor relativní permeability a % je tenzor elektrické susceptibility, který udává, jak snadno se materiál vlivem vnějšího elektrického pole polarizuje.
2
Teorie elektromagnetického pole
3
Dosazením rovnice (1.10) do rovnice (1.6) přepíšeme závislost elektrické indukce D na elektrické intenzitě E pomocí standardnější veličiny relativní permitivity éT = 1 + X-
D = e0érE (1.10)
B = /ioArH., (1.11)
Veličiny a, fiľ a éT jsou jedinečné pro každý materiál, hovoříme o nich jako o materiálových parametrech, a jejich znalost nám postačuje k vyřešení Maxwellových rovnic v materiálovém prostředí.
Aby se elektromagnetické pole šířilo materiály spojitě, je třeba, aby na rozhraní materiálů byly splněny následující lokální podmínky, které spojitost zaručují. Uvádíme rovnice příslušející rozhraní dvou materiálů, kde materiály odlišujeme indexováním:
(Da-DO-n   =   pf1, (1.12) (B2-Bi)-n   =   0, (1.13) nx(E2-E!)   =   0, (1.14) nx(H2-H!)   =   Jfm, (1.15)
(1.16)
kde pf1 symbolizuje volný náboj na rozhraní materiálů a J™ značí volný proud mezi materiály; n je jednotkový vektor kolmý k rozhraní materiálů.
1.2   Optické vlastnosti látek
Optické vlastnosti látek ovlivňují šíření elektromagnetických vln k jejichž popisu se obvykle používá komplexní index lomu N a komplexní dielektrická funkce er. Tyto dvě optické veličiny jsou svázány vztahy
N = n + ik =       = V^i + ie2, (1-17)
kde jednotlivé komponenty lze navzájem přepočítat pomocí matematických úprav. Uvádíme vyjádření komponent dielektrické funkce v závislosti pouze na složkách komplexního indexu lomu lze psát
2      i 2
£i = n — k , é2   = Ink.
(1.18)
(1.19)
Jelikož veličiny N a er popisují tytéž optické vlastnosti, je pouze na zvyku, které se pro popis kterých fenoménů používají. Například pro vyjádření reflexe nebo transmise elektromagnetického záření bývá preferován popis pomocí komplexního indexu lomu N. Tatáž veličina se upřednostňuje k popisu šíření elektromagnetické vlny prostředím, kde n vystupuje ve vztahu pro fázovou rychlost vlny a k popisuje útlum vlny během šíření. Naopak dielektrická funkce se volí při popisu rozptylu na částicích, jejichž velikost je malá ve srovnání s vlnovou délkou záření, nebo při popisu mikroskopických mechanismů zodpovědných za optické jevy.
FDTD solver GSVIT
4
Na obrázku 1.1 vidíme spektrální průběh reálné a imaginární části dielektrické funkce pro typické kovy používané v plazmonice ve viditelném spektru: zlato a stříbro.
Obrázek 1.1: Závislost reálné (plná čára) a imaginární části (přerušovaná) dielektrické funkce pro zlato (červeně) a stříbro (modře) [2]. Šipky naznačují příslušnost křivek k jednotlivým osám.
Kapitola 2
FDTD solver GSVIT
V této kapitole jednoduše přiblížíme fungování metody FDTD a zaměříme se na základní strukturu solveru GSVIT, kterou v dalších kapitolách budeme podrobněji komentovat.
Metoda FDTD je založena na diskretizaci Maxwellových rovnic (1.4), (1.5) a (1.12) v prostoru i čase. Na ID příkladě budeme demonstrovat fungování této metody [3].
Uvažujme lineárně polarizovanou elektromagnetickou vlnu šířící se ve směru osy z, jejíž vektor elektrické intenzity E kmitá ve směru x a vektor magnetické intenzity H ve směru y. Nejprve je třeba prostor i čas diskretizovat, k tomuto účelu zavedeme indexy m, n, pomocí nichž rozdělíme prostor (ve směru osy z) a čas tak, že m-tý prostorový krok a n-tý časový krok budou s polohou z a časem t souviset vztahy
K diskretizaci elektromagnetického pole nahradíme parciální derivace v Maxwellových rovnicích (1.4), (1.5) a (1.12) centrálními diferencemi (předvádíme pro elektrickou intenzitu E, pro magnetickou intenzitu H by vztahy vypadaly obdobně),
2.1   Metoda FDTD
z = mAz,
t = nAt.
(2.1)
(2.2)
dE
dz dE
AE
Az
(2.3)
Az
AE
Az
(2.4)
dt
At
Po této náhradě vypadaj Maxwellovy rovnice (1.4) a (1.5) následovně:
Enx+\m)-E^{m)
(2.5)
Aí
(m+§)-ffr*(™+|) Aí
ExH{m + l)-E2H{m)
HqAz
e0Az
(2.6)
5
FDTD solver GSVIT
6
Rovnice (2.5) a (2.6) se následně střídavě řeší v jednotlivých prostorových a časových krocích. Střídání řešení rovnic v prostorových a časových krocích vede k tomu, že časový vývoj jednoho pole v následujícím kroku je počítán z prostorového vývoje druhého pole v kroku předcházejícím.
Zaměříme-li se na vývoj pole v určitém místě v prostoru, lze hodnotu v každém časovém kroku určit ze znalosti polí v předchozím časovém kroku, viz. obrázek 2.1. Toto bude zřejmé, pokud vyjádříme pole pomocí diskretizovaných Maxwellových rovnic (2.5) a (2.6) v posledním časovém kroku
1/2 3/2 5/2 7/2
Ex
Obrázek 2.1: Demonstrace výpočtu elektromagnetického pole Yeeho algoritmem, kdy se střídají v jednotlivých časových krocích výpočty elektrického a magnetického pole. Obrázek převzat z [4].
2.1.1   Volba velikosti kroku
Velikost kroku je pro FDTD zásadním parametrem, ovlivňuje totiž nejen kvalitu získaných výsledků, ale taktéž náročnost a délku výpočtu. Obecně platí, že každý rozměr buňky by měl být menší, než desetina nejmenší v simulaci použité vlnové délky Amin,
(2.9)
FDTD solver GSVIT
Dále je třeba dbát na to, aby na simulované struktury připadlo dostatečné množství buněk. V neposlední řadě pamatujme na to, že diskretizací jsme prostor rozdělili na kvádry, takže pokud z nich chceme vysvavět zakřivený povrch, je nutné taktéž uzpůsobit (zjemnit) krok.
Časový krok se řídí tzv. Courantovou podmínkou
Aí<^, (2.10) 2c
kde c značí rychlost světla ve vakuu.
2.2 GSVIT
GSVIT[5] je jeden z mnoha solverů k řešení Maxwellových rovnic v časové doméně. Jeho výhodou oproti komerčním solverům (např. Lumerical[6]) je jeho přístupnost, neboť jde o open source solver. Tato výhoda je ovšem znevážena horším grafickým rozhraním pro nastavení řešené úlohy.
XSVIT, grafické rozhraní GSVITu, ve skutečnosti generuje pouze textový soubor .par, který se následně spouští pomocí binárního programu, což přináší možnost si parametrové vstupní soubory snadno skriptovat.
2.3   Jak seznam parametrů vytvořit
Seznam parametrů doporučujeme vytvářet v libovolném textovém editoru (případně v grafickém prostředí XSVIT). Po naplnění dokumentu kýženými parametry, přidáme při ukládání souboru příponu .par, díky čemuž bude k obsahu přistupovat GSVIT jako k seznamu parametrů.
Aby GSVIT věděl, který parametr si má kam zařadit, je důležité dodržovat určitou syntax (viz dále), která se opírá o několik hesel, psaných kapitálkami.
Zde uvádíme možnou ukázku parametrového souboru main.par:
PDDL
200   200   200   le-06   le-06 le-06
CDMP 800
VERBOSE 1
THREADS -1
GPU 0
UGPU 0
SDURCE_TSF
10   10   10   190   190   190  0  0  0   1     5.5e-07 1.0
BDUNDARY_ALL
FDTD solver GSVIT
8
liao
DUT_IMAGE
Ey   10   -1   -1   100 Ey_XY DUT_FILE
mfmball_depGWY_mfm_10pixelle - 6 . gwy
MEDIUM_VECTDR ball.txt
Spuštění simulace
Spuštění simulace je možné pomocí příkazové řádky zadáním následujícího příkazu, v ukázce naznačujeme spuštění souboru main.par, jehož naplnění předvedeme v následujících kapitolách,
steindl:   r- $   . /gsvit   main . par
případně z grafického rozhraní XSVIT.
Kapitola 3
COMPUTATIONAL DOMAIN
POOL
PDDL
200 250 200  le-06 le-06 le-06
Heslo POOL odkazuje na definici výpočetního objemu.
První tři argumenty (označeny červeně) rozdělí výpočetní kvádr na daný počet segmentů v daném směru; (v ukázce postupně 200 segmentů ve směru x, 250 ve směru y a 200 ve směru z).
Další tři argumenty (označeny modře) nastavují velikost jednoho segmentu v každém směru; (v ukázce jde o segment - krychličku o délce hrany le-06 metru).
9
Kapitola 4
BASIC PARAMETERS
COMP
CDMP 800
Heslo COMP charakterizuje výpočetní čas, jako argument přijímá počet časových kroků (v ukázce jsme nastavili 800 časových kroků).
DT_MULT
DT_MULT 0.5
Heslo DT_MULT nastavuje multiplikativní faktor pro časový krok. Přidání tohoto parametru je především nutné pro pokročilé modely kovů.
Výše jsme zvolili hodnotu multiplikativního faktoru 0.5, čímž jsme dvakrát zkrátili velikost časového kroku. VERBOSE
VERBDSE 1
Heslo VERBOSE je spojeno s množstvím vypsaných informací od GSVITu v průběhu výpočtu. VERBOSE = 0 nevypisuje nic.
VERBOSE = 1 vypisuje pouze číslo kroku, který probíhá.
VERBOSE > 1 vypisuje informace o stavu dokončení jednotlivých fází (především stav započtení okrajových podmínek pro jednotlivé složky elmag. záření).
10
BASIC PARAMETERS
11
THREADS
THREADS -1
THREADS souvisí s počteni použitých jader k výpočtu. Lze nastavit hodnotu -1, která odpovídá automatické detekci počtu jader a následnému využití všech automaticky dohledaných jader.
GPU
GPU
0
GPU souvisí s počtem použitých grafických karet, při nastavení GPU =< @0@ > probíhá výpočet na procesoru počítače.
Kapitola 5
BOUNDARY CONDITIONS
Při výpočtu FDTD hrají zásadní roli okrajové podmínky. V GSVITu se dají nastavit okrajové podmínky pro celý početní objem stejného typu BOUNDARY.ALL, pro každou stěnu výpočetního kvádru podmínky jiné BOUNDARY, případně podmínky periodické MBOUNDARY.
Z počátku pojednáme o možných typech okrajových podmínek.
5.1 type
Typ none odpovídá hranici bez odrazu, tedy představuje rozhraní mezi výpočetním objemem a ideálním elektrickým vodičem.
Typ liao odpovídá druhému řádu absorpce, poskytuje základní absorpci, ovšem ne tak kvalitní jako okrajové podmínky cpml.
Typ cpml (Convolutional Perfectly Matched Layer) očekává parametry depth, power, crmax, "maj, Kmax a pomocí nich dopočítává absorbci. Tato metoda je navržena tak, aby rozhraní neodráželo žádné vlny zpět do výpočetního objemu.
Parametry crmax, amax, Kmax slouží ke generování souřadnic protahování a polynomální tlumení stupně power na vrstvě tloušťky depth.
Při nastavení crmax = — 1, bude parametr crmax optimalizován. Aby byl výpočet stabilní, musí být parametry crmax, amax kladné a Kmax > 1.
BOUNDARY
V této sekci se zaměříme na zápis podmínek pro jednotlivé hranice výpočetního objemu. Postupně jsou ukázány zápisy pro různé typy okrajových podmínek.
12
BOUNDARY CONDITIONS
13
BDUNDARY_XO liao
BDUNDARY_XN liao
BDUNDARY_YO
cpml    10 3 -1 0.03 4
BDUNDARY_YN
cpml  10 3 -1 0.03 4
BDUNDARY_Z0 none
BDUNDARY_ZN none
Výpočetní objem bude tedy ohraničen ve směru osy x okrajovými podmínkami liao, ve směru y cpml (s parametry depth=10, power=3, crmax=-l, amax=0.03, Kmax=4), ve směru z vlna volně opouští na obou koncích výpočetní objem.
BOUNDARY_ALL
V případě, že chceme pro všechny hranice nastavit stejný typ okrajových podmínek, použijeme příkaz BOUN-DARY_ALL. Ke spuštění potřebujeme jen upřesnit typ okrajových podmínek (typy okrajových podmínek jsou uvedeny v sekci type). Uvádíme příklad použití liao podmínek (okrajové podmínky none).
BDUNDARY_ALL liao
MBOUNDARY
GSVIT umožňuje nastavit i umělé periodické uvnitř výpočetního objemu. Ukážeme zápis pro umělé podmínky vložené na pozici x = 10 a x = 90, námi nenastavené hodnoty jsou automaticky nastaveny na okraj výpočetního objemu.
MBDUNDARY_X0 periodic 10
MBDUNDARY_XN periodic 90
Zde je třeba poznamenat, že nastavení typu MBOUNDARY.ALL, jak bychom mohli intuitivně očekávat, které funguje pro nastavení všech okrajových podmínek na stejný typ, není podporováno.
Kapitola 6
SOURCE
Total/Scattered field source TSF_SOURCE
SDURCE_TSF
area_from area.to *angles mode
Rovinná vlna (TSF) má 4 základní argumenty - odkud kam se v objemu šíří, pod jakým úhlem dopadá na objem a vlastnosti záření.
Tedy na místě area_from (x_start, y_start, z_start) a area_to (x_end, y_end, z_end) se očekávají souřadnice k vymezení objemu, ve kterém se bude záření šířit.
*angles (6, (p, ip) první dva úhly určují orientaci šíření vlny výpočetním objemem, viz obrázek 6.1. V případě nastavení nulových úhlů dopadá záření kolmo na stěnu výpočetního kvádru a není nutné tento argument zadávat. Třetí argument udává polarizaci propagující rovinné vlny. Uhly je třeba zadávat v radiánech!
Abychom určili charakter záření, lze volit ze tří módů:
mode = 0 umožňuje načíst průběh elektrického pole ze souboru, je nutné, aby data v souboru byla uspořádána do dvou sloupců, kde první sloupec odpovídá jednotlivým časovým krokům (zde je třeba dbát, aby v souboru bylo obsaženo alespoň stejně hodnot, kolik jsme zadali časových kroků v COMP). Tato varianta si vyžádá zadat do zdrojového souboru.
SDURCE_TSF
area^from area.to *angles 0 název_souboru
V mode = 1 si GSVIT vytvoří sinusový signál sám, postačí mu dodat vlnovou délku [v metrech] a amplitudu [ve Voltech] záření. Do zdrojového souboru tedy zapíšeme.
SDURCE_TSF
area_from area.to *angles 1 vlnová_délka amplituda
mode = 2 je obdobou mode = 1 s gaussovským tlumením. Signál si opět GSVIT vytvoří sám, musíme mu jen dodat vlnovou délku záření [v metrech], šířku gaussovské obálky [v počtu časových kroků] a amplitudu [ve Voltech] záření. Do zdrojového souboru se zapisuje
14
SOURCE
15
direction	polarisation	theta [deg]	phi [deg]	psi [deg]	zi /<	k E \	
x axis	y	90	0	0			
x axis	z	90	0	90			r k
y axis	x	90	90	0			
y axis	z	90	90	90			—► y
z axis	x	0	0	90	*ŕ x		
z axis	y	0	0	0			
Obrázek 6.1: Obrázek zachycující smysl úhlů a způsob nastavení polarizace podél příslušné osy. Obrázek převzat
z dokumentace GSVIT [5].
SDURCE_TSF
area_from area.to *angles 2 vlnová_délka šířka_gauss_obálky amplituda
TSF.SKIP
TSF_SKIP umožňuje nastavit chování záření na rozhraní s objemem, kterým jsme záření omezili. Tento příkaz rozhoduje, zda se na určité hraně uplatní okrajové podmínky, nebo nikoliv. V následující ukázce vynecháváme okrajové podmínky ve stěně yz připadající x = x_start (iO) a ve stěně xy připadající z = z_end (in).
TSF_SKIP iO
TSF_SKIP ku
SOURCE.POINT
SDURCE_PDINT position mode angles
Pro šíření rovinné vlny z bodového zdroje z místa position, je nutné kromě souřadnic (xq, yo, zo) dodat také úhly angles (6, cp) [úhly odpovídají obr. 6.1], do kterých se záření šíří a charakterizovat záření pomocí volby modu.
Mody rozlišujeme, podobně jako u TSF_SOURCE.
mode = 0 umožňuje načíst časový průběh elektrického a magnetického pole ze souboru, je nutné, aby data v souboru byla uspořádána do sloupců v pořadí (čas, Ex, Ey, Ez, Hx, Hy, Hz). I zde dbejte, aby v souboru bylo obsaženo alespoň stejně hodnot, kolik jsme zadali časových kroků v COMP). Tato varianta si vyžádá zadat do zdrojového souboru.
SDURCE_PDINT
position 0 název_souboru angles
SOURCE
16
V mode = 1 si GSVIT vytvoří časovou závislost se sinusovým průběhem sám (hodnoty uloží do aktuální složky pod názvem tmpsource), postačí mu dodat vlnovou délku [v metrech] a amplitudu [ve Voltech] záření. Do zdrojového souboru tedy zapíšeme
SDURCE_PDINT
position 1 vlnová.délka amplituda angles
mode = 2 je obdobou mode = 1 s gaussovským tlumením. Signál si opět GSVIT vytvoří sám, musíme mu jen dodat vlnovou délku záření [v metrech], šířku gaussovské obálky [v počtu časových kroků] a amplitudu záření [ve Voltech]. Do zdrojového souboru se zapisuje
SDURCE_PDINT
position 2 vlnová.délka šířka_gauss_obálky amplituda angles
Kapitola 7
MEDIA
FDTD simulací sledujeme obvykle interakci elmag. záření s vloženým předmětem. Tento předmět je charakterizován tvarem a materiálem. V této části se zaměříme, jak GSVITu vnutit námi požadovaný tvar a materiál předmětu, který budeme sledovat. V počátku si ukážeme, jak definovat materiál pixel po pixelu (MEDIUM_LINEAR), následně se budeme věnovat definování materiálu pomocí řádky parametrů (MEDIUM_VECTOR).
7.1   Vytvoření souboru pro zkoumaný předmět
GSVIT umožňuje vkládat do výpočetního objemu zkoumaný předmět ve dvou režimech. MEDIUM_LINEAR, který je vhodný pro výpočet komplikovaných nebo plynule se měnících struktur. Pro jednodušší geometrické objekty se doporučuje použít režim reprezentovaný řádkou parametrů MEDIUM_VECTOR.
Dále budeme předpokládat, že soubory medium_linear.txt a medium.txt jsou umístěny ve stejném adresáři jako hlavní soubor main.par, pokud tomu tak není, bude nutné název materiálového souboru nahradit v následujícím celou cestou k tomuto souboru.
MEDIUM_LINEAR
MEDIUM_LINEAR přijímá jako argument název souboru (medium_linear.txt), který je vyplněn binárními daty reprezentující materiál pixel po pixelu.
MEDIUM_LINEAR medium_linear .txt
Strutura souboru musí být následující (symboly reprezentují sloupečky):
position     e     a     n <jm
position (3 sloupečky reprezentující souřadnice voxelů v x, y a z) - segmentace v x, y, z musí odpovídat počtu voxelů, na které jsme v části 3 rozdělily výpočetní objem. V těchto bodech musí být vyčísleny materiálové parametry relativní permitivita e, relativní permeabilita a, elektrická /i a magnetická vodivost om- Souřadnice musí být reprezentovány integeremi v 32-bitovém zápisu, materiálové parametry fioaty v 32-bitovém zápisu.
17
MEDIA
18
MEDIUM.VECTOR
Doporučujeme oddělit hlavní tělo parametrů od souboru týkajícího se předmětu vloženého do elmag. pole. To provedeme založením nového souboru, např. medium, txt. Do hlavního souboru parametrů main. par zaneseme informaci, že GSVIT má hledat parametry zkoumaného předmětu v souboru medium.txt následovně (zde předpokládáme, že soubor medium.txt je umístěn ve stejném adresáři jako hlavní soubor main.par, pokud tomu tak není, je zde třeba vyplnit celou cestu k souboru medium.txt)
MEDIUM_VECTDR medium.txt
Dále veškeré informace v této kapitole budeme zapisovat do souboru medium.txt. Informace zapisované do souboru mají následující charakter.
code_shape shape_volume_def material
7.1.1 code_shape
Pomocí code_shape nadefinujeme, jaký tvar chceme do výpočetního objemu vložit.
Tabulka 7.1: Kody náležící k různým tvarům vkládaným do GSVITu. * Gwyddion výškové pole.
tvar	koule	kvádr	válec	kužel	uříznutý kužel	čtyřstěn	soubor	Gwyddion*
code_shape	4	8	7	107	108	12	21	22
7.1.2 shape_volume_def
V předchozím jsme vybrali tvar, který vložíme do výpočetního objemu, zde mu nadefinujeme ukotvení v prostoru a rozměrově ho vymezíme.
Koule
Koule je v prostoru charakterizována umístěním středu (x, y, z) a jejím poloměrem (rádius).
4 x y z rádius materiál
Kvádr
Kvádr charakterizumeje pomocí souřadnic dvou bodů, které jsou vůči sobě úhlopřičně umístěny. Jeden bod nechť charakterizují souřadnice (xO, yO, zO), druhý bod (xN, yN, zN). Píšeme potom:
8 xO yO zO xN yN zN material
MEDIA
19
Válec
Válec udejme pomocí pozice středu jedné podstavy (xO, yO, zO), středu druhé podstavy (xN, yN, zN) a poloměru podstav rádius.
7 xO yO zO xN yN zN rádius materiál
Kužel
Kužel zadejme pomocí pozice vrcholu (xO, yO, zO), středu podstavy (xN, yN, zN) a jejím poloměrem rádius.
107 xO yO zO xN yN zN rádius materiál
Uříznutý kužel
Uříznutý kužel je do výpočetního prostoru umístěn pomocí definování středů (xO, yO, zO), (xN, yN, zN) a poloměrů radiusO, radiusN obou podstav.
108 xO yO zO xN yN zN radiusO radiusN materiál
Čtyřstěn
Čtyřstěn definujeme pomocí souřadnic všech jeho vrcholů.
12 xO yO zO xl yl zl x2 y2 z2 x3 y3 z3 materiál
7.1.3   Načtení ze souboru
Soubory obsahující informace k načtení tvaru musí být ve formátu .node a .ele (v ukázce f ilebase. node, f ilebase . ele). K vytvoření využijeme příslušných programů. Abychom mohli nastavit správnou velikost a umístění v prostoru, posunutí v prostoru zajistíme pomocí xshift yshift zshift, velikost naškálujeme pomocí xmult ymult zmult. Soubory f ilebase .node, f ilebase. ele mohou obsahovat indexy materiaLindex, kterými se definuje podmnožina bodů, která bude sloužit jako materiál. Tohoto parametru se dá využít v případě, že soubory obsahují několik indexů a každému se přiřadí jiný materiál. Píšeme:
21 filebase atribute_num material_index xshift yshift zshift xmult ymult zmult materiál
Zde je příklad kodu odpovídající použití materiaLindex.
21 cube.l 01000 111 materiall 21  cube.l 00000111 material2
# my   little cube 6     4 1
# elements   (tetrahedrons),   links  to   corresponding  node list
#
1 6 7 15 1
MEDIA
20
2 1 8 3 4 0
3 3 7 1 2 1
4 1 7 6 2 0
5 1 7 3 8 1
6 7 8 1 5 0
Symbolem l značíme příslušnost k materiálu material 1, podtržítko zde je jen pro zvýraznění. 7.1.4 material
Použití materiálů ukážeme pro jednoduchost na kuličce o poloměru 10 voxelů vložené na pozici 50, 50, 50. vakuum
V případě, že budeme vkládat do výpočetního objemu objekt z vakua, použijeme
4 50 50 50 10 0
ideální elektrický vodič (pec)
4 50 50 50 10 10
lineární prostředí
U lineárního prostředí je třeba zadat relativní permitivitu e, relativní permeabilitu /i a elektrickou a a magnetickou vodivost (tm-
4 50 50 50 10  1 e a fi <rM
CP model
Drudeho model s kritickými body (CP model, 5) s předpisem permitivity [7, 8]
ecpM=eoo-   2,P.     +2-AMň-^~ + o   .     xt     ■ i7"1)
U!z + vyuj \ S lp — U! — ll p       S lp + U! + ll p J
vyžaduje zadání statické permitivity e^, parametry související s Drudeho modelem (plasmová frekvence cjp, frekvence tlumení 7) a dvě sady popisující kritické body (amplituda A, fáze cp, rezonanční frekvence uj a frekvence tlumení V).
5 50 50 50 10 4 £co uP 7 AI 01 ul TI A2 02 u2 T2
Pro zlato G SVIT přebírá parametry z [7]
MEDIA
21
4 50 50 50 10 5 1.11683 0 13.1839el5 -1.18299 3.88123el5 0.452005el5	0	,0207332el5	3	.04155 -1.09115	4	20737el5	2	35409el5	0	273221
a pro stříbro										
4 50 50 50 10 5 0.89583 0 13.8737el5 -1.48944 6.15009el5 0.659262el5	0	,0207332el5	1	.3735 -0.504659	7.	59914el5	4	28431el5	0	304478
načtení z databáze
GSVIT využívá informací z SOPRA (www. sspectra. com/sopra. html) databáze a z webových stránek https: //refractiveindex. info . Volbou 99 určíme, že chce vybírat materiál z databáze GSVITu a následným argumentem vybereme konkrétní sloučeninu. Seznam dostupných sloučenin v databázi nejlépe zjistíte, když rozkliknete adresář, kam jste GSVIT nainstalovali a dále cestou /gsvit/gsvit-1.8. l/data/spectra dojdete do databáze. V ukázce vybíráme z databáze A10N.
4 50 50 50 10 99 A1DN
Do databáze lze přidávat další materiály, v takovém případě doporučujeme postupovat podle http://gsvit.net/ spectra.php, kde je problematika popsána.
7.2 MATMODE_CHECK
MATMDDE_CHECK 1
Při volbě 1 zkontroluje nastavení materiálu a pokud je to možné, alokuje menší část paměti. Především se jedná o použití jednodušších rovnic v případě nemagnetického materiálu nebo při výpočtech ve vakuu, čímž lze zkrátit výpočetní čas. Defaultně je tento parametr nastaven na 0.
Kapitola 8
OUTPUTS
V této části popíšeme sondy integrované do GSVITu. Poznamenejme, že výběr správných sond a jejich vhodné nastavení je zásadní pro námi simulovaný jev. Dodejme, že vypisování výsledků stojí nějaký čas, proto se snažíme vypisování omezit a soustředit se pouze na podstatné veličiny.
OUT_FILE
Příkaz OUT_FILE slouží k definování výstupního souboru, do kterého má GSVIT zapisovat 2D data. Pro soubor se volí formát Gwyddionu, uveďme příklad se zápisem do souboru gwyout.gwy.
DUT_FILE gwyout.gwy
OUT.POINT
Chceme-li sledovat vybranou složku elmag. pole v určitém místě, využijeme funkce OUT_POINT, které specifikujeme, o kterou složku elmagjype a v kterém bodě position se zajímáme, s jakým časovým krokem time_step a do jakého souboru filename se má zapisovat.
DUT_PDINT
elmag_type nskip position filename
GSVIT umožňuje sledovat v libovolném bodě výpočetního objemu jednotlivé složky elmag. pole. Jednotlivé značky elmag. pole najdeme v tab. 8.1.
Tabulka 8.1: Tabulka symbolů odpovídajícím jednotlivým složkám elmag. pole.
složka	X	y	z
elektrická	Ex	Ey	Ez
magnetická	Hx	Hy	Hz
všechny	All		
Mimo chování jednotlivých složek elmag. vlny umožňuje GSVIT v tomto výstupu zaznamenávat i protékající proud, relativní permitivitu a permeabilitu a vodivost. Odpovídající značení pro zdrojový kod je v tab. 8.2.
22
OUTPUTS
23
Tabulka 8.2: Další možné parametry sledování v případě OUT_POINT režimu výstupu.
permitivita	Epsilon
permeabilita	Mu
vodivost	Sigma
Sigast	Sigast
Uvádíme příklad zápisu všech složek v bodě zájmu na pozici (50, 50, 50) se zápisem do souboru all_output.txt vždy po 10 hodnotách a zápis pouze elektrické x-ové složky na stejné pozici se zápisem do ex_output.txt po 5 hodnotách. V případě zápisu pomocí více OUT_POINT (ať už z důvodu více bodů zájmu nebo sledování průběhu více složek, je nutné volit různé výstupní soubory, jinak dochází k přepisování dat)!
DUT_PDINT
All 10 50 50 50 all_output.txt Ex 5 50 50 50 ex.output. txt
OUTJMAGE
DUT_IMAGE
elmag_type timestep position description
OUTJMAGE je funkce zapisující do zvoleného gwy souboru, který byl zvolen v rámci OUT_FILE, 2D data odpovídající řezu výpočetním objemem podél jedné souřadnicové osy. Tedy v případě nastavení řezu například na x = 10, musíme zbylé souřadnice nastavit na hodnotu -1. Argument elmag_type se volí dle popisu v části OUT_POINT.
Ukážeme příklad pro nastavení řezu v místě x = 10, kdy dochází k zápisu každých 5 časových kroků. Zapisovat Ex komponentu pole budeme do souboru gwyout.gwy, kde po otevření Gwyddionem budou tato data popsána jako
Ex_out.
DUT_FILE gwyout . gwy
DUT_IMAGE
Ex 5  10 -1 -1 Ex.out
OUT.PLANE
je obdobou OUTJMAGE. Opět zaznamenává 2D data, ale ukládaje do formátu raw, odpovídající řezu výpočetním objemem (obdobně jako výše). Nadefinujeme o kterou složku máme zájem pomocí elmag_type (diskutováno výše). Dalšími argumenty vybereme časový krok timestep ukládání, začátek a konec detekce timestart timestop. Argument ASCI=1 zapisuje výsledky v ASCI kodu.
DÁLE NUTNO OVĚŘIT FUNKČNOST.
DUT_PLANE
elmag_type timestep ASCI timestart timestop position filename
OUTPUTS
24
OUT.SUM
Tato varianta zaznamenává součet elektrické intenzity pole (buď jednotlivé složky Ex, Ey, Ez a nebo všech dohromady All) nebo absorpci Abs. Sumování probíhá ve vymezeném kvádru s danými materiálovými konstantami e,
er, [i a <jm-
DUT_SUM
elmag_type timestep position e a fi um filename
Pozice kvádru je určena pomocí souřadnic dvou úhlopříčně umístěných bodů. Na místo position vkládáme následující
xO yO zO xN yN zN
OUT_FORCE
Pro zachycení optické síly působící na definovaný objem (xO yO zO xN yN zN) slouží OUT_FORCE. Hodnota přiřazená každému bodu definovaného objemu odpovídá časovému průměru přes periodu zdroje. Optická síla je pak dopočtena pomocí Maxwellova stress tenzoru.
DUT_FDRCE
timestep xO yO zO xN yN zN filename
Kapitola 9
NEAR FIELD to FAR FIELD
Transformace blízkého pole na vzdálené pole je metoda používaná k výpočtu hodnoty pole mimo výpočetní objem. Pro použití není ovšem nezbytné, aby bod počítaný jako v dalekém poli ležel mimo výpočetní objem, proto můžeme porovnávat hodnoty získané pro daleké pole s hodnotou získanou z Yeeho algoritmu (to je úkolem příkladu na http://gsvit.net/xexl6_4.php).
NFFF
Příkaz NFFF definuje objem, ze kterého se bude počítat transformace NEAR FIELD to FAR FIELD. Funkce vyžaduje zadat souřadnice dvou bodů (xO yO zO xN yN zN) umístěných úhlopříčně a začátek sum_from a konec sum_to integrace.
NFFF
xO yO zO xN yN zN sum_from sum_to
NFFF_SKIP
Tento příkaz určuje oblast (hranici), která bude z NFFF výpočtu vyloučena. Pomocí parametrů vybereme o kterou hranici jde a následně z ní můžeme 'vykrojit' obdélník (vybereme ho pomocí 4 souřadnic ze 6 iO jO kO iN jN kN, kde první dvě hodnoty značí horní levní roh a další dvě hodnoty dolní pravý roh obdélníku), který se nezapočíta do NFFF výpočtu.
NFFF_SKIP
boundary iO JO iN jN
Argument boundary volíme dle tab. 9.1.
Tabulka 9.1: Značení pro hranice.
hranice	xO	yO	zO	xN	yN	zN
značení	iO	JO	kO	iN	JN	kN
25
NEAR FIELD to FAR FIELD
26
NFFF_RAMAHI_POINT
Umožňuje zvolit bod, do kterého se přesunou hodnoty vypočtené pomocí NFFF výpočtu. Je třeba definovat pozici bodu (x y z) a název souboru filename pro výpis hodnot elektrického pole. Tento bod může být umístěn i mimo výpočetní objem.
NFFF_RAMAHI_PDINT x y z filename
Kapitola 10
Příklad: Odraz na rovinném rozhraní
V této kapitole se zaměříme na konkrétní simulaci v GSVITu. Vybrali jsme řešení problému odrazu elektromagnetického záření na rovinném rozhraní. Budeme sledovat odrazivost a lokální zesílení pole na tomto rozhraní. Výsledky z GSVIT simulace porovnáme s analytickým řešením.
10.1    Zadání problému
Vyšetřete spektrální průběh odrazivosti na rovinném rozhraní vzduch-kov (zlato nebo stříbro) a studujte lokální zesílení elektrického pole přímo na tomto rozhraní.
Lokální zesílení LFIEF se definuje jako podíl celkové intenzity na rozhraní a dopadajícího záření
LFIEF=|E^+^'2, (10.1)
kde (E'a) značí elektrickou složku dopadající (odražené) vlny, přičemž pro kolmý dopad se dá vztah vyjádřit pouze jako podíl dielektrických konstant (index M odpovídá kovu, A vzduchu)
LFIEF = 4
(10.2)
10.2   Parametrový soubor
V této části se zaměříme na řešení navrženého problému pomocí GSVITu. Ke spuštění vytvoříme soubor parametrů planarsurface.par.
PDDL
60 60 120 le-9 le-9 le-9
CDMP
9001
DT_MULT
0.8
27
Příklad: Odraz na rovinném rozhraní
28
VERBDSE 0
THREADS
-1
GPU
0
SDURCE_TSF
10 10 15 50 50 220 0 0 0  1 340e-6 1
TSF_SKIP JO
TSF_SKIP jn
TSF_SKIP iO
TSF_SKIP in
TSF_SKIP kn
BDUNDARY_ALL liao
MBDUNDARY_XO periodic 15 MBDUNDARY_XN periodic 45 MBDUNDARY_YO periodic 15 MBDUNDARY_YN periodic 45
DUT_FILE
p lan ar _ surf ace . gwy DUT_PDINT
Ey 100 30 30 10 outputR.txt DUT_PDINT
Ey 100 30 30 50 output_surface.txt
MEDIUM_VECTDR deska.txt
Parametrový soubor 10.1: Ukázka souboru parametrů planarsurface.par pro GSVIT - detekce odrazivosti a
zesílení na povrchu rozhraní vzduch-kov.
Nastavení výpočetního objemu
V první řadě musíme vytvořit výpočetní objem, to znamená nastavit velikost objemu a jeho členění, což se dělá v části POOL a počet časových kroků, které se nastavují v části COMP.
V našem příkladu jsme vytvořili objem o 60 x 60 x 120 voxelech, kde jsme voxely volili jako krychličky o hraně
1 nm (le-9).
V části COMP jsme zvolili 9001 časových kroků. Zkušenosti uživatelů simulujících úlohy s kovovými materiály přiměly vývojáře přidat od verze GSVIT 1.8.2 parametr DT_MULT. Jedná se vlastně o multiplikátor regulující délku časového kroku. Pro kovové materiály je třeba krok zkrátit, proto volba 0.8.
Příklad: Odraz na rovinném rozhraní
29
Vypisování a hardware
Parametr VERBOSE umožňuje nastavit, kolik informací má GSVIT vypisovat v průběhu simulací. Parametr je kaskádován od 'Vypiš všechno' (3) až po 'Nevypisuj nic' (0).
Volba THREADS '-1' nastaví maximální počet dostupných jader (GSVIT si jejich počet zjistí sám). GPU umožňuje pouštět GSVIT na grafických kartách. Zvolením '0' spustíme výpočet na procerosu.
Zdroje záření
Aby se definovaným objemem mohla šířit elektromagnetická vlna, musíme ji do objemu vložit. GSVIT disponuje několika typy elektromagnetického zdroje, více informací najdete v kapitole 6 nebo dokumentaci [5]. Pro naši úlohu jsme použili postupnou rovinou vlnu, kterou charakterizuje oddíl SOURCE_TSF. Pomocí 10 10 15 50 50 220 jsme prostorově vymezili šíření elektromagnetické vlny. Další tři argumenty 0 0 0 značí úhly (v radiánech), pomocí nichž se definuje směr a polarizace elektromagnetické vlny, dle obrázku 6.1.
Poslední tři číselné argumenty umožňují nastavit mód, vlnovou délku a aplitudu vlny. Módy 0/1/2 odpovídají postupně možnosti načtení profilu elmag. záření ze souboru/ sinusový signál/sinusový signál s gaussovským tlumením. Pro úlohu jsme nechali vygenerovat GSVIT sinusovou vlnu (1) o vlnové délce wl = 340e — 6 (v metrech) a amplitudou 1 (ve Voltech).
Podotkněme, že GSVIT umí simulovat pouze jednu vlnovou délku, chce-li uživatel získat spektrální závislost, musí simulaci pustit opakovaně, vždy pro jednu vlnovou délku.
Okrajové podmínky
Jelikož metodou FDTD řešíme diferenciální rovnice, je řešení zásadně ovlivněno nastavením okrajových podmínek. Celý výpočetní objem jsme obalili absorpčními okrajovými podmínkami liao. Dále jsme díky symetrii úlohy zavedli do objemu periodické okrajové podmínky MBOUNDARY, kterými jsme ve směrech x a y vymezili objem (mezi voxely 15 a 45 v obou směrech).
TSF.SKIP vyloučí zadanou stěnu TSF objemu z aplikace TSF.SOURCE. Výstupy
Abychom získali informace o rozložení a intenzitě elektromagnetického pole, musíme do objemu zavést vhodné sondy. Pro účely naší úlohy byly zvoleny dva výstupy OUT_POINT, do kterých se zapisuje zvolená komponenta elektromagnetického pole (Ey) v závislosti na čase vyčítaná z určitého místa v prostoru (30 30 10). Data se každých 100 časových kroků zapíší do souboru outputR.txt.
Výstup outputR.txt zapisuje odraženou vlnu v dalekém poli, výstup output_surface.txt zaznamenává vlnu v blízkém poli (na povrchu rozhraní). Aby bylo možné výstupy vypisovat, potřebuje mít GSVIT zadaný OUTJľlLE.
Příklad: Odraz na rovinném rozhraní
30
Vložené materiály
Abychom mohli studovat interakci elektromagnetické vlny s materiálem, je třeba ho do výpočetního objemu vložit, k tomu slouží oddíl MEDIUM_VECTOR. Standardně se jako argument používá textový soubor, ve kterém jsou uložena všechna specifika, co se materiálů týká.
8  10 10 50 50 50 210 model materiál
Parametrový soubor 10.2: Ukázka souboru deska.txt pro definování materiálů použitých v simulaci.
GSVIT má předdefinované základní geometrické útvary, které do něj lze vložit, těmto útvarům odpovídá první argument v materiálovém souboru (8 - kvádr). Tento geometrický útvar se dále umisťuje do prostoru sadou argumentů (u kvádru jde o pozice vrcholů ležících na tělesové úhlopříčce, 10 10 50 50 50 210). Dále je nutné přiřadit geometrickému útvaru správné fyzikální konstanty, aby ho elektromagnetická vlna vnímala jako vyrobený ze specifického materiálu. K tomu je třeba vybrat model, dle kterého si GSVIT přepočítá materiálové parametry. Každý model si vynucuje nějaké vlastní parametry, kterých je zapotřebí k simulování konkrétního materiálu. Pro naše simulace jsme použili Drudeho model s kritickými body (CP model, model=5), který vyžaduje zadání materiálových parametrů na místě materiál GSVIT přebírá z [7], pro stříbro (následující sekvenci je třeba zadat do deska.txt na místo materiál):
0.89583 0 13.8737el5 0.0207332el5 1.3735 -0.504659 7.59914el5 4.28431el5 0.304478 -1.48944 6.15009el5 0.659262el5
Parametrový soubor 10.3: Materiálové parametry pro stříbro [7].
a pro zlato:
1.11683 0 13.1839el5 0.0207332el5 3.04155 -1.09115 4.20737el5 2.35409el5 0.273221 -1.18299 3.88123el5 0.452005el5
Parametrový soubor 10.4: Materiálové parametry pro zlato [7].
Spuštění simulace
Spuštění simulace je možné pomocí příkazové řádky zadáním následujícího příkazu,
steindl:   r- )in$   ./gsvit planar_surface.par
případně prostřednictvím grafického rozhraní XSVIT.
10.3   Výsledky simulace
Nyní srovnáme výsledky simulace pro kolmý dopad elektromagnetické vlny na rovné rozhraní vzduch-kov s analytickou formulí (10.2).
GSVIT nám po spuštění simulace zapíše do souboru dva sloupečky, které odpovídají časovému kroku a (v našem případě, viz nastavení OUT_POINT) elektrické intenzitě Ey.
Ze získaných dat dopočteme maximální intenzitu v daném místě, čímž získáme reflektivitu i zesílení na kovovém povrchu (LFIEF = \Esnr + ^ur|2/|£sur|2).
Příklad: Odraz na rovinném rozhraní
31
Obrázek 10.1: Komponenta Ey po 9000 časových krocích v simulovaném objemu pro úlohu rovného rozhraní vzduch-stříbro. Použité elektromagnetické záření má vlnovou délku 340 nm. V obrázku je dále naznačena geometrie a uspořádání úlohy. Číselné údaje v obrázku jsou ve voxelech nebo nanometrech.
Obrázek 10.1 zachycuje rozložení komponenty Ey po 9000 krocích simulace. Největší (absolutní) hodnotu vidíme u rozhraní. Stejný obrázek zachycuje i grafický návrh úlohy.
V obrázku 10.2 sledujeme srovnání výsledků z FDTD simulace s analytickým výpočtem (rovnice (10.2)) a literaturou [9]. Simulace s analytickým výsledkem se velmi dobře shodují, to je způsobeno jednak jednoduchostí úlohy a především použitím stejných materiálových parametrů [2]. Při pohledu na křivky odpovídající Schlůckerově knize [9] vidíme posunutí struktur ve spektru stříbra. Toto posunutí vůči ostatním výsledkům souvisí s tím, že Schlůcker použil jiné optické konstanty pro stříbro.
V obrázku 10.2b) je vyneseno lokální zesílení LFIEF, které je obecně pro rovné rozhraní vzduch-kov menší než 1 ve srovnání se situací, že bychom zde rozhraní neměli. Ideální bezeztrátový kov by vytvořil na povrchu pole, které by vyrušilo dopadající vlnu (Ei + ET = 0), muselo by dojít tedy k absolutnímu odrazu. Obecně tedy můžeme říct, že malé LFIEF souvisí s velkou odrazivostí. Tento trend je patrný i na obrázku 10.2, kde LFIEF > 1 pouze pro stříbro okolo 330 nm, kde zároveň stříbro přestává dobře odrážet světlo.
Příklad: Odraz na rovinném rozhraní
32
1.0 -
0.8 -
1 1	1 1	i	i i
,//	y____v-//:' ///		-
			Au Ag
i i			
- - Analytical
— GSVIT ...... Schlucker
200 300
400     500     600 700 Wavelength (nm)
800     900 1000
Obrázek 10.2: Stříbro modře, zlato červeně: a) Porovnání reflektivity z FDTD simulace získané pomocí GSVITu (plná čára) s tabulovanou reflektivitou [2] (přerušovaně) a literaturou [9] (tečkované), b) Srovnání lokálního zesílení na povrchu LFIEF rovného rozhraní vzduch-kov: plná čára odpovídá FDTD výpočtu, přerušovaná analytickému výpočtu dle rovnice (10.2) s parametry dle [2], tečkovaná čára značí výsledek ze Schlůckerovy učebnice[9].
Literatura
[1] J. Clerk Maxwell. A Dynamical Theory of the Electromagnetic Field. Phil. Trans. R. Soc. Lond., 115:459-512, 1865.
[2] P. B. Johnson and R. W. Christy. Optical constants of the noble metals. Phys. Rev. B, 6:4370-4379, Dec 1972. doi: 10.1103/PhysRevB.6.4370. URL https://link.aps.Org/doi/10.1103/PhysRevB.6.4370.
[3] D. W. Ward and K. A. Nelson. Finite-difference time-domain (FDTD) simulati- ons of electromagnetic wave propagation using a spreadsheet. Comput. Appl. Eng. Educ, 13:213-221, 2005.
[4] Obrázek - princip Yeeho algoritmu. zhnotes.wordpress.com/2013/08/09/ one-dimensional-fdtd-simulation-of-laser-pulse-in-vacuum, 2017.
[5] GSVIT 1.8.2. http://gsvit.net/index.php, 2017.
[6] Lumerical. https ://www. lumerical. com/, 2017.
[7] Chun Kyungwon, Kim Huioon, Kim Hyounggyu, and Chung Youngjoo. PLRC and ADE implementations of Drude-critical point dispersive model for the FDTD method. Progress In Electromagnetics Research, 135: 373-390, 2013.
[8] A. Vial and T. Laroche. Comparison of gold and silver dispersion laws suitable for FDTD simulations. Applied Physics B, 93:139-143, June 2008.
[9] Sebastin Schliicker. Surface Enhanced Raman Spectroscopy. WILEY-VCH Verlag Co. KGaA, 2011.
33