Анализ дифракции света на элементах цилиндрической микрооптики объединенным методом конечных элементов Галеркина и граничных элементов

АНАЛИЗ ДИФРАКЦИИ СВЕТА НА ЭЛЕМЕНТАХ ЦИЛИНДРИЧЕСКОЙ МИКРООПТИКИ ОБЪЕДИНЕННЫМ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ ГАЛЕРКИНА И ГРАНИЧНЫХ ЭЛЕМЕНТОВ

Д.В. Нестеренко, В.В. Котляр Институт систем обработки изображений Российской академии наук Самарский государственный аэрокосмический университет

Аннотация

Рассматривается задача дифракции плоской электромагнитной волны на двумерном (цилиндрическом) объекте с поглощением, размеры которого сравнимы с длиной волны. Для приближенного решения этой задачи разработан объединенный метод конечных элементов Галеркина и граничных элементов. Проводится исследование зависимости относительной погрешности объединенного метода на примере моделирования дифракции на диэлектрическом и проводящем цилиндрах.

Введение

Для задач моделирования дифракции света на оптических элементах в однородном пространстве с размерами порядка длины волны свет должен рассматриваться как электромагнитное излучение. Численные методы моделирования дифракции могут классифицироваться как дифференциальные [1], разностные [2, 3], интегральные [4-8, 30], вариационные [9-12], дискретных источников [13], рассмотрение хода лучей [14].

Для определения электромагнитных полей в точке пространства интегральные методы объединяют вклады в эту точку от поля источников по объему или на поверхности. Популярность интегральных методов основывается на их способности решать неограниченные полевые задачи, т.к. условие излучения Зоммерфельда, безусловно, удовлетворяется в формулировке задачи. Более того, интегральные методы требуют знания поля только на поверхности дифракционного элемента, а не полного поля в пространстве, что минимизирует число неизвестных. В работах [15, 16] представлен объединенный метод на основе метода граничных элементов. В [17] разработан метод расчета дифракции на плоскопараллельной пластине с неоднородностью на основе интегральных уравнений, связанный с численным решением тензора Грина. Недостаток таких методов в том, что они приводят к полностью заполненным матрицам и, следовательно, требуют большего объема компьютерной памяти и длительное время вычисления. Также границу применения методов необходимо брать по физической границе объекта. Если неоднородность имеет нетривиальную форму, то это приводит к увеличению числа неизвестных.

Разностное решение дифференциальных уравнений Максвелла рассматривалось в работах [18-21, 3]. В работе [22] описано разностное решение волнового уравнения. Недостатками такого подхода являются невозможность использования условий излучения, ограничения на шаги сетки. Для моделирования стационарных задач прохождения излучения разностными схемами используется конечное количество длин волн падающего импульса, что искажает спектр волны. Использование в качестве граничных условий для неограниченных задач дифракции

условий поглощающей границы [23, 10] позволяет приближенно решать уравнения Максвелла разностными схемами, точность решения зависит от количества слоев на искусственной границе и от степени ее замкнутости.

В отличие от методов разностного решения системы уравнений Максвелла интегральные и вариационные методы не требуют конструирования сложных поглощающих граничных условий [2, 10].

Вариационные методы в задачах с ограниченной областью задачи определяют решения уравнения Гельмгольца путем минимизации функционального соотношения. В работе [10] уравнение Гельмгольца решалось методом конечных элементов Галеркина с использованием граничных условий сложного вида, зависящих от неизвестного параметра, что потребовало применение границы, определенной формы. Кроме того, данный метод также не включает в себя условия излучения Зоммерфельда.

В работе [24] представлен гибридный метод на основе метода конечных элементов, сформулированного через метод Ритца, и метода граничных элементов. В данном гибридном методе метод конечных элементов применяется для решения уравнения Гельмгольца во внутренней части неоднородного диэлектрического элемента микрооптики и применяется интегральная методика, метод граничных элементов, к области, внешней к элементу, где должны удовлетворяться условия излучения. Оба метода соединяются на границе между внешней и внутренней частями, с удовлетворением условий непрерывности поля. Использование метода конечных элементов для определения поля внутри объекта приводит к матрице трехдиагонального вида, что требует меньше компьютерной памяти и времени вычисления, чем методы объемных интегралов [25]. Результатом использования метода граничных элементов для определения поля на границе является более точное решение, чем применение метода конечных элементов с условиями поглощающей границы. Но применение метода Ритца к решению уравнения Гельмгольца некорректно, т.к. он накладывает требование положительности оператора решаемого уравнения. Вывод о знакоопределенности оператора уравнения Гельмгольца сделать нельзя.

В предыдущих работах авторов [26-29] разработаны методы анализа и синтеза элементов двухмерной (цилиндрической) диэлектрической микрооптики, основанные на методе, объединяющем метод конечных элементов Галеркина (МКЭГ) и метод граничных элементов (МГЭ). В данной работе проводится исследование разновидности объединенного метода МКГЭ-МГЭ для случая дифракции ТЕ- и ТМ-поляризованных волн на проводящих и диэлектрических цилиндрических элементах микрооптики.

Описание метода расчета В рассматриваемой задаче источник, находящийся в пространстве, освещает цилиндрическую структуру. В отсутствии структуры этот источник создал бы падающее поле. При наличии структуры он создает другое поле, называемое полным полем. Рассеянное поле определяют как разность между полным полем и падающим полем. Цель задачи -определить полное или рассеянное поле, характеризующие структуру.

Любое двумерное поле может быть разложено на Ег-поляризованное и Яг-поляризованное поля. В области дифракции поле описывается системой дифференциальных уравнений, различных для случаев ТЕ- и ТМ-поляризаций. Для ТЕ-поляризации

(Е (х, у) = (0,0, Ег (х, у))) комплексная амплитуда

и(х,у) обозначает полное электрическое поле Ег(х,у), которое направлено вдоль оси г (вдоль образующей цилиндрического оптического элемента), координаты (х,у) лежат в плоскости нормального сечения.

Для ТМ-поляризации (н(х,у)=(0,0,Яг(х,у))) комплексная амплитуда и(х,у) обозначает полное магнитное поле Яг(х,у).

Полное поле иа(х,у) в области О должно удовлетворять уравнению:

-1-Уи„ (х, у) +

Р

V-

+К Ч (х> у) иа (х> у) = /а

где /а = А10 Л

ТЕ-поляризации и

Р (х, у)=цг , ч (х, у)=ег для

/а=-

1 У,

д_

+ ду

' 1 У."

Р (х, у)=ег, Ч (х, у)=Цг

для ТМ-поляризации. Константы и гг представляют собой отношение магнитной и диэлектрической проницаемостей среды к аналогичным показателям свободного пространства, т.е., цг=ц/ц0 и ег=е/е0, к0 представляет собой волновое число волны в свободном пространстве

, / \ 1/2 ю 2п

к0 =ю(^0е0 ) = - = 7", С Ап

10 = -\/н-0 / е0 - импеданс свободного пространства, У - плотность электрического тока источника.

В рассматриваемой задаче область расчета бесконечна. Однако, как известно, МКЭ применим только к конечной или ограниченной области. Таким образом, чтобы решить уравнение (1), бесконечная область внешняя к рассеивателю, должна быть ограничена введением искусственной границы Г. Соответственно, для единственного решения задачи на данной искусственной границе должны быть введены граничные условия. Такие условия должны сделать границу как можно более прозрачной для рассеянного поля или, другими словами, должны минимизировать нефизические отражения от границы. Один из классов граничных условий, разработанных для этих целей, может быть получен из граничных интегральных уравнений, применяемых к внешней области. Эти граничные условия являются глобальными по своей природе, т.е. они соотносят поле в одном граничном узле с полем на всей границе. Данные граничные условия предотвращают отражение на границе для всех углов падения волн и приводят к точному решению.

Таким образом, требуется определить полное поле и(х,у) в областях а (внутренняя) и Т (внешняя) удовлетворяющее указанным условиям.

Метод Галеркина решения уравнения (1) основан на решении соотношений вида:

Я(-"Р Аиа^- Чк'иаУ- /ау) аа = 0

где у - произвольная функция из области определения уравнения (1).

Используя первую формулу Грина:

Л РАд аа = | р ^ПП а/ = JJvрvg аа,

а г ап а

для функций Р и д, где а - область плоскости х, у; Г - ее граница, обходимая против часовой стрелки;

ад

ап

производная в направлении внешней нормали

к кривой Г, получим:

^аГ!^ -р VuQ (ху )VY- чк \ (х у аг = ГГ /аУ аа.

а V Р

у duа

-Г

г Р а

Систему базисных функций для а обозначим

( а / , N

{юк/ (х, у)} и систему базисных функций для Г

обозначим {ю^ (х,у)}М1, где Ых, Ыу - число узлов

сеточного покрытия прямоугольной области а по оси х и у, соответственно, М - число узлов сеточного покрытия границы Г.

Подставляя в соотношение (4) вместо произвольной функции у систему базисных функций для метода Галеркина можно записать систему линейных уравнений:

а

Ли + Бу = С/, где и = (и!, ..., иЫхЫу ) коэффициентов {и

вектор, составленный из

разложения:

иП(х у )= Е ик,, (х у).

к ,1=0

Вект°р / = ( f1... /МШу )Т - вектор^

ный из коэффициентов разложения:

/«(у) = Е л,<,(х,У).

составлен-

Хотя равенства (6) и (7) действительны для всех точек (х, у) в области «, необходима отдельная обработка величин поля и его частных производных на границе Г от значений во внутренней области. Разложение, аналогичное (6) и (7), для поля и его частных производных на границе имеет вид:

иГ (х у ) = Е (х У^

т=1

М

уГ

т=1

М

Л Г(х, у ) = Е ЛЮ (х, у),

где (х, у) е Г, у=(у1, ..., уМ) - вектор, составленный из коэффициентов разложения ук = дик/ дп.

Элементы матрицы А вычисляются по формулам:

аА,

Ыу(к )+1, Ыу (/)+у

=л

1

Р

дю?,

дх

дх

дюк,

ду ду

- к02 д

к,/ = [1,мх], I,у = [1,ку],

где у - область разбиения области включающая узлы сети к иу.

Элементы матрицы В вычисляются по формулам:

К* = - $ ютю; ^

Гт ,

т,;

= [1,М ]

где Гт, ; - линейная область границы Г, включающая узлы границы т и ;.

Элементы матрицы С вычисляются по формулам:

СМу

«к,,

к, I = [1, Ых ], ,, у = [1, N ],

где у - область разбиения области включающая узлы сети к и у.

В качестве кусочно-линейного базиса была определена функция вида:

I. (х, у) =

1 "к 1 -1 + 1 + 1+ 1 -

хк - х у, - у

к

хк - х

к

у, - у

к

хк - х

к

хк - х

к

у, - у

к

у, - у к '

к

если х, у е «

если х, у е £2Ы2

если х, у е « к

если х, у е «к

если х, у е « к

если х, у е

где к - длина сегмента сети покрытия.

Элементы (ак1) матрицы А, элементы (Ьт;) матрицы В и элементы (е'к1) матрицы С вычисляются

из уравнений (11), (12) и (13), соответственно. Тогда систему уравнений (10) можно представить в виде:

и«

[ Л«,«] [ Л Г,«] 0

[ Л «,Г] [ Л Г,Г] [Б ]

[С О,«] [СГ,«] [С«,Г] [СГ,Г]

Система уравнений (15) не имеет единственного решения, так как она состоит из N равенств с Ы+М неизвестными: Ы=ЫхЫу - общее число узловых величин поля ик,,(х, у) в области «и М производных по нормали на граничных узлах ук,,(х, у).

Определим поля в свободном пространстве ¥ (вне области ^ с ее границей Г). Так как это пространство однородное, то поле может быть сформулировано в терминах граничных интегралов с соответствующей функцией Грина. Полное поле и¥(х, у) в области ¥ должно удовлетворять уравнению:

V-

- Vuт © Р

+ кО2 дир (£,) = , ^еР,

где = укогоЛ , Р (х у) = , д (х у) = ег для ТЕ-поляризации и

/р=-

д I 1

-Л

+±[ ^ л4

ду г

дх ^е

р (х, у)=ег, д (у) = ^г

для ТМ-поляризации. Л ¥ - плотность электрического тока источника в свободном пространстве.

Для построения граничного интегрального уравнения для поля и его нормальной производной вве-

т=1

а

дем функцию Грина и , удовлетворяющую условиям излучения Зоммерфельда и являющуюся фундаментальным решением уравнения Гельмгольца:

V2 и'

Фундаментальное решение для уравнения Гельмгольца в двумерном однородном пространстве известно и равно

и ' =

где г ^[(л) -^(О]2 +[[(п) - I)] 2 ,

Н0(1) (кг) = 30 (кг + гУ0 (кг)) - функция Ханкеля первого рода и нулевого порядка, где 30 - функция Бесселя первого рода нулевого порядка, У0 - функция Неймана нулевого порядка.

Для построения граничного интегрального уравнения для рассеянного поля и его нормальной производной воспользуемся теоремой Грина в следующем виде:

I [V V (I, п) + к2 и* (I, п)] и(п) ат =

г

+1 и(п)4* (I, п)аг+1 (п)и* (I, п) ат.

г т

Функции в обоих интегралах в правой части уравнения (19) - это д(п) = ди(п)/дп' - нормальные производные амплитуды поля.

Подставляя равенство (17) в уравнение (19) и осуществляя предельный переход точки наблюдения 4 из внутренней в граничную, найдем

с (I) и (I) = -| (п) и* (I, п)ат + +| д (п) и* (I, п)аг -1 и (п) д* (I, п) аг.

Это уравнение обеспечивает функциональную связь между функциями и и ее нормальной производной д на границе Г. Функция с в уравнении (20) равна:

с(0 = 1 --Пф ,

2п

где ф - внутренней угол кусочно-линейной границы в точке I. Первое слагаемое в правой части является полем, создаваемым источником / в свободном пространстве, и может быть обозначено как падающее поле и £ .

Таким образом, уравнение (20) записывается в следующем виде:

с (I) и (I) = | д (п) и* (I, п)аг-

г (22) и (п)д*(I,п)аг + и"©.

г

Подставляя вместо функции комплексной амплитуды в уравнении (22) на границе ее аппроксимацию базисными кусочно-линейными функциями

при (8) и (9), получим соотношение (23), которое может быть представлено в матричном виде (24) с элементами матриц [Б] и [в] вида (25) и (26).

Интегралы в равенствах (25) и (26) могут быть оценены численно. Объединяя системы уравнений (15) и (24), получим замкнутую систему линейных алгебраических уравнений для решения задачи дифракции плоской волны на цилиндрическом микрообъекте (27).

СА =Ё к У* |г( + [Р*+1 - Р* ]УК»

*=1 [ [ 0

1

»(Рт , Р* + [Р*+1 - Р* ]У)аУ +|юг (р*-1 +

0

+ ([Р* - Р*-1 ]'у) и* (Рт , Р*-1 + [Р* - Р*-1 ]Т)ат] --и* 10 г(р* + [Р*+1 - Р* ]у)х

(Рт . Р* + [Р*+1 - Р* ]У) а дп'

+10 г (-1 + [р *- р *-1]у)х

,,ди' (Рт > Р*-1 + [Р * - Р *-1]У ) дп'

т = [1, М ]. [Б] йг+[в] уг = .

ау

dm. = - к

+ [Р*+1 - Р*] у)>

Хди' (Р т > Р * + [Р*+1 - Р*] У) дп'

+£*-1 + [Р* -Р*-1]у)х

,,ди'

дп'

ёу

- Ст 5т*>

= к

.К

Хи* (Рт = Р* + [Р Р* ] У)ёу + 1

+/г (Р*-1 + [Р* -Р*-1]у)»

0

Хи* (Рт = Р*-1 + [Р* - Р*-1] У)ёУ] т, * = [1, М ].

[Аа.п] [Аг,О] 0 "

[ Ап.г] [[ Аг,г] [В]

0 [О] [С]

[Со,о] [Сг>п] 0 [Со,г] [[Сг,г] 0 0 0 [Е]

где подматрица Аап размерностью (М - М)*(М - М) включает в себя коэффициенты соотношений поля во внутренних узлах сети разбиения. Подматрицы Ааг и АГ,П размерностью (М - М) *М и М*(М - М), соответственно, включают в себя коэффициенты связи поля в граничных и внутренних узлах. Подматрица АГ, Г размером М*М включает в себя коэффициенты связи поля в граничных узлах. Подматрица В размером М*М включает в себя коэффициенты соотношений между производными поля на границе и полем во внутренних узлах сети разбиения. Подматрица Б размером М*М включает в себя коэффициенты связи поля свободного пространства в граничных узлах. Подматрица О размером М*М включает в себя коэффициенты соотношений между производными поля на границе и полем свободного пространства во внутренних узлах сети разбиения. Подматрица Сп,п размерностью (М - М)*(М - М) включает в себя коэффициенты соотношений источников поля во внутренних узлах сети разбиения. Подматрицы СП,Г и СГ,П размерностью (М - М) *М и М*(М - М), соответственно, включают в себя коэффициенты соотношений источников поля в граничных и внутренних узлах. Подматрица СГ, Г размером М*М включает в себя коэффициенты соотношений источников поля в граничных узлах. Подматрица Е - единичная матрица размером М*М Векторы ип и иГ - вектора напряженности поля во внутренних и граничных узлах сети, уГ - векторы нормальных производных поля в граничных узлах сети. Векторы и 1Г - векторы источников поля во внутренних и граничных узлах сети, и 'г" - вектор напряженности внешнего падающего поля в граничных узлах сети. Таким образом, размерность системы уравнений (27) - (М + М)х(М + М).

После определения значений поля и его производных на границе Г поле в любой точке области ¥ определяется соотношением (22), где с(|)=1.

Результаты численного моделирования Рассмотрим дифракцию плоской волны (ТЕ- и ТМ-поляризаций) на диэлектрическом и проводящем однородных цилиндрах с круглым сечением для экспериментального исследования сходимости объединенного метода. Сходимость алгоритма зависит от длины сегмента к, длины волны источника X, и гг среды. Поскольку магнитная и диэлектрическая проницаемости среды являются переменными в задаче, рассмотрим зависимость решения от параметра X/к , определяющего количество сегментов сети на одной длине волны. Пусть плоская волна падает на цилиндр, длина волны Х0 =0,5 мкм. Радиус цилиндра равен 0,25 мкм. Относительная диэлектрическая проницаемость проводящего цилиндра из алюминия е= - 41,4+/11,9. Относительная диэлектрическая проницаемость диэлектрического цилиндра е=2,25. Окружающее цилиндр однородное пространство имеет параметры е=ц=1. Расчет производился объединенным методом МКЭГ-МГЭ. В качестве области О выбиралась квадратная область, контур Г - периметр области О. Источники внутри

О отсутствуют. Область О была покрыта квадратной сетью 105^105 узлов. Процесс расчета занял 10 минут на персональном компьютере Реп1шт4.

На рис. 1 представлены результаты численного моделирования дифракции ТЕ- и ТМ-волны на диэлектрическом цилиндре. На рис. 2 - результаты моделирования дифракции на проводящем цилиндре.

Рис. 1. Распределение интенсивности поля дифракции на диэлектрическом цилиндре

Рис. 2. Распределение интенсивности поля дифракции на проводящем цилиндре

Для оценки дифракционных процессов используем диаграмму направленности рассеяния, зависящей от угловой координаты ф, определенную в бесконечно удаленных точках (р^-да) как

ст(ф) = Нт 2пр-

р^да I

В частности, диаграмма направленности определялась в прямом направлении (ф=0), в обратном направлении (ф=п) и в поперечном направлении (ф=п/2) при р = 10000Х. Ниже приведены зависимости ст/Х от параметра X/к . Результаты представлены

в дБ (стдБ=10 ^10ст). На рис. 3 представлены результаты оценки дифракции ТЕ- и ТМ-волны на диэлектрическом цилиндре.

На рис. 4 представлены зависимости значений диаграммы направленности дифракции ТЕ- и ТМ-волны на проводящем цилиндре.

На рис. 5 приведена зависимость от параметра X/ к относительного отклонения значений диаграммы направленности рассеяния ТЕ- и ТМ-волны на диэлектрическом цилиндре от значений при X/к = 100 для углов ф равных 0, п/2, п. На рис. 6 приведена зависимость от параметра X/к относительного отклонения значений диаграммы направленности рассеяния ТЕ- и ТМ-волны на проводящем цилиндре.

и

2

и

а)

о/к,дБ 30

20 10 0 -10 -20

í\ — ф=о° Ф =90° Ф=180°

—

\

ГЦ 1

—

0

20

40

60

80 k/h

б)

о/к,дБ 30

20 10 0 -10 -20

\ — ф=0° Ф=90° Ф=180"

1 \ —

\

Л 1 i

1

1 !

0

20

40

60

80 k/h

а)

о/к,дБ 30 20 10 0 -10

Рис. 3. Зависимость диаграммы направленности рассеяния на диэлектрическом цилиндре от параметра Я/И для ТЕ-поляризации

о/к,дБ

í\ -ф=0° -Ф =90° --- ф=180°

i \

н

0

20

40

60

80 k/h

б)

30 20 10 0 -10 -20

-Ф=0°

-Ф=90°

---ф=180"

0

20

40

60

80 k/h

Рис. 4. Зависимость диаграммы направленности рассеяния на проводящем цилиндре от параметра Я/И для ТЕ-поляризации

Ла/а,%

а)

0,30 0,25 0,20 0,15 0,10 0,05

0

-Ф=0° -Ф=90° ---ф=т°

I

I 1

j \

\

1 ____ ----

Ла/а,% 0,30

0,25

0,20

0,15

0,10

0,05

1 -Ф=0° -Ф=90° --- ф=180°

1

1

1 1 U 1 у\ |

1\ i

20

20

40

60

а)

Ла/а,% 0,30

0,25

0,20

0,15

0,10

0,05

0

40 60 80 k/h б) 0

Рис. 5. Зависимость относительного отклонения значений диаграммы направленности рассеяния на диэлектрическом цилиндре от параметра Á/h для ТЕ-поляризации

Ла/а,% 0,30

0,25

0,20

0,15

0,10

0,05

80 k/h

1 1 -Ф=0° -ф=90° --- ф=180°

1 1 -L

1 i 1

л V

г \

1 \ Ы

-Ф=0° -Ф=90° --- ф=180°

[\ 1

1

1 1

1 1 1 (

\

20 40 60 80 Ш б) 0 20 40 60

Рис. 6. Зависимость относительного отклонения значений диаграммы направленности рассеяния на проводящем цилиндре от параметра Я/И для ТЕ-поляризации

Серия экспериментов с диэлектрическим цилиндром показала, что относительное отклонение значений диаграммы направленности рассеяния становится менее 5 % при ИИ>40 и менее 1% при Х/И>80 для обеих поляризаций. Эксперименты с проводящим цилиндром по-

казали, что относительное отклонение значений диаграммы направленности рассеяния становится менее 5% при УН>30 для ТЕ-поляризации и при УН>50 для ТМ-поляризации и менее 1% при Х/И>50 для ТЕ-поляризации и при Я/И > 80 для ТМ-поляризации.

Таким образом, состояние поляризации практически не влияет на результаты моделирования диэлектрических структур предложенным методом, но оно должно быть учтено при выборе длины сегмента сети покрытия для расчета проводящих структур с соответствующей относительной погрешностью.

Заключение

В работе получены следующие результаты. Разработан объединенный метод конечных элементов Галеркина и граничных элементов, с помощью которого можно моделировать дифракцию плоской электромагнитной волны на неоднородных двумерных (цилиндрических) диэлектрических и проводящих объектах микрооптики, размер которых сравним с длиной волны света.

Проведено исследование зависимости относительной погрешности результатов моделирования от отношения длины волны падающего света к длине сегмента сети покрытия.

Благодарности

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» (грант CRDF RUX0-014-SA-06), а также грантов РФФИ 05-0850298, 07-07-97601 и 07-07-97600.

Литература

1. Montiel F., Neviere M. Differential theory of gratings: extension to deep gratings of arbitrary profile and permittivity through the R-matrix propagation algorithm // J. Opt. Soc. Am., 1994. - Vol. 11. - P. 3241-3250.

2. Taflove A. Computational electromagnetics: the finite-difference time domain method // Artech House, Boston, 1995.

3. Головашкин Д.Л., Дегтярев А.А., Сойфер В.А. Моделирование волноводного распространения света оптического излучения в рамках электромагнитной теории // Компьютерная оптика. - Самара, 1997. - В. 17.

4. Brebbia C.A. The boundary Element Method for Engineers // Pentech Press, London; Halstead Press, New York, 1978 (Second edition, 1980).

5. Choi M. K. Numerical calculation of light scattering from a layered sphere by the boundary-element method // J. Opt. Soc. Am., 2001. Vol. 18. No. 3.P. 577-583.

6. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. - М.: Мир, 1987. - 524 с.

7. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. - М.: Мир, 1987. - С. 104-107.

8. Бреббия К., Уокер С. Применение метода граничных элементов в технике. - М.: Мир, 1982. - 248 с.

9. Davies J. B. Finite element analysis of waveguides and cavities - a review // IEEE Trans. Magn., 1993. - Vol. 29. - P. 1578-1583.

10. Lichtenberg B., Gallagher N. Numerical modeling of dif-fractive devices using the finite element method // Opt. Eng., 1994. - Vol. 33. - No. 11. - P. 3518.

11. Wei X., Wachters A.J.H., Urbach H.P. Finite-element model for three-dimensional optical scattering problems // J. Opt. Soc. Am. A., 2007. - V. 24. - No. 3. - P. 866-881.

12. Koshiba M., Saitoh K. Finite-element analysis of birefringence and dispersion properties in actual and idealized holey-fiber structures // App. Opt., 2003. - Vol. 42. - No. 31.

13. Blaike R. J., McNab S. J. Evanescent interferometric lithography // App. opt., 2001. - Vol. 40. - No. 4. -P. 1692-1698.

14. Voznesensky N. Simulation model for light propagation through nanometer-sized structures // Optical Memory and Neural Networks, 2000. - Vol. 9. - No. 3. - P. 175-183.

15. Prather D. W., Mirotznik M. S., Mait J. N. Boundary integral methods applied to the analysis of diffractive optical elements // J. Opt. Soc. Am., 1997. - Vol. 14. -P. 34-43.

16. Tanaka M., Tanaka K. Computer simululation for two-dimensional near-field optics with use of a metal-coated dielectric probe // J. Opt. Soc. Am., 2001. - Vol. 18. -No. 4. - P. 919-925.

17. Paulus M., Martin O. J. F. Light propagation and scattering in stratified media: a Green's tensor approach // J. Opt. Soc. Am., 2001. - Vol. 18. - No. 4. - P. 854-861.

18. Dou W. B., Yung E. K. N. Diffraction of an electromagnetic beam by an aperture in a conducting screen // J. Opt. Soc. Am., 2001. - Vol. 18. - No. 4. - P. 801-806.

19. Lee J.-F., Palandech R. and Mittra R. Modeling three-dimensional discontinuities in waveguides using nonorthogonal FDTD algorithm // IEEE Trans. Microwave Theory Tech., 1992. - Vol. 40. - P. 346-352.

20. Prather D. W., Shi S. Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements // J. Opt. Soc. Am., 1999. - Vol. 16. - No. 5. - P. 1131-1142.

21. Shi S., Tao X., Yang L., Prather D. W. Analysis of diffrac-tive optical elements using a nonuniform finite-difference time-domain method // Opt. Eng., 2001. - Vol. 40. -No. 4. - P. 503-510.

22. Gruzdev V., Gruzdeva A. Finite-difference timedomain modeling of laser beam propagation and scattering in dielectric materials // Proceedings of SPIE, 2001. - Vol. 4436. - P. 27- 38.

23. Berenger G.P. A perfectly matched layer for the absorption of electromagnetic waves // J. of Comp. Phys., 1994. -Vol. 114. - P. 185-200.

24. Mirotznik M., Prather D., Mait J. A hybrid finite element-boundary element method for the analysis of diffractive elements // Journal of Modern Optics, 1996. - Vol. 43. -No. 7. - P. 1309-1321.

25. Ильинский А.С., Кравцов В.В., Свешников А.Г. Математические модели электродинамики // - М.: Высшая школа, 1991. - 223 c.

26. Kotlyar V.V., Nesterenko D.V. A finite element method in the problem of light diffraction by micro-optics // Optical Memory and Neural Networks, 2000. - Vol. 9. - No. 3. -P. 209-219.

27. Kotlyar V.V., Nesterenko D.V. Analysis of light diffraction by binary micro-optics using a combination of boundary element method and finite element method // Proceedings of SPIE, 2001. - Vol. 4242. -P. 125-132.

28. Котляр В.В., Нестеренко Д.В. Дифракция электромагнитной волны на круговом диэлектрическом цилиндре: расчет по аналитическим формулам и методом конечных элементов - граничных элементов // Физика волновых процессов и радиотехнические системы. - Самара: ПГАТИ, 2000. - Т. 3. - В. 3-4. - С. 25-28.

29. Нестеренко Д.В., Котляр В.В. Объединенный метод конечных элементов Галеркина и граничных элементов для анализа дифракции ТМ-поляризованной плоской волны на цилиндрических оптических элементах // Компьютерная оптика. - Самара: ИСОИ РАН, 2002. - В. 24. - С. 17-25.

30. Colton D., Kress R. Integral equation methods in scattering theory // John Wiley&Sons, New York, 1983.