Распространение радиально-ограниченных вихревых пучков в ближней зоне. Часть i. алгоритмы расчёта

РАСПРОСТРАНЕНИЕ РАДИАЛЬНО-ОГРАНИЧЕННЫХ ВИХРЕВЫХ ПУЧКОВ В БЛИЖНЕЙ ЗОНЕ.

ЧАСТЬ I. АЛГОРИТМЫ РАСЧЁТА

Хонина С.Н.1'2, Устинов А.В.1, Ковалёв А.А.1'2, Волотовский С.Г. 1 1 Учреждение Российской академии наук Институт систем обработки изображений РАН, 2 Самарский государственный аэрокосмический университет имени академика С.П. Королёва

Аннотация

На примере дифракции плоской волны на круглой апертуре в ближней зоне (порядка нескольких длин волн) проведено сравнение алгоритмов расчёта с использованием векторного интегрального преобразования Рэлея-Зоммерфельда (РЗ) и разложения по плоским волнам (ПВ) по точности и скорости вычислений.

В скалярном случае, что соответствует вычислению поперечных компонент электрического поля, результаты различаются только в очень близкой к апертуре области.

В векторном случае при расчёте продольной компоненты в методе ПВ возникает особенность в области спектральных частот, радиус которых близок к единице. Предложены различные варианты обхода этой особенности. На расстоянии нескольких длин волн результаты двух рассматриваемых алгоритмов совпадают и отличаются от полученных с помощью конечно-разностного временного метода (FDTD) только масштабно (среднеквадратичное отклонение с учетом масштаба составляет менее 2%). Таким образом, рассмотренные в данной работе алгоритмы позволяют получать за существенно меньшее время структурно верную (но несколько завышенную по амплитуде) картину дифракции в ближней зоне. Такая амплитудная «завышенность» может быть связана с тем, что рассмотренные методы РЗ и ПВ не подразумевают наличия у-компоненты в дифракционной картине изначально х-поляризованного поля. Во второй части статьи рассматривается модификация метода ПВ, позволяющая учесть наличие всех векторных компонент.

Ключевые слова: дифракционный интеграл Рэлея-Зоммерфельда, разложение по плоским волнам, дифракция на круглой апертуре, вихревой пучок.

Введение

В связи с уменьшением размеров оптических устройств большое внимание в последнее время уделяется описанию непараксиального распространения световых полей [1-9] и разработке алгоритмов моделирования такого распространения [10-21].

Дифракция на круглой апертуре плоской волны является тестовым примером для разрабатываемых алгоритмов как случай, имеющий аналитическое решение на оптической оси [22, 23], и исследовалась с различной степенью точности во многих работах [11, 24-29]. Таким образом, данный пример может использоваться для начальной оценки точности рассматриваемых в данной работе алгоритмов. В данной работе дифракция на апертуре моделируется как распространение в свободном пространстве ограниченных по радиусу в начальной плоскости полей с постоянной амплитудой и вихревой фазой. Частным случаем таких полей является ограниченная плоская волна.

Непараксиальная скалярная модель, основанная на теории Рэлея-Зоммерфельда (РЗ) [30], позволяет получать согласующиеся с экспериментами результаты на очень близких расстояниях от апертуры [3134]. В частности, в работе [33] на примере дифракции на цилиндрических тонких фазовых объектах (что практически соответствует скалярной задаче дифракции на щели) было показано, что в зависимости от фазового набега различные типы дифракционных интегралов - Кирхгофа, Рэлея-Зоммерфельда I или II типа - позволяют получить корректные результаты на расстоянии всего полдлины волны от

объекта. На расстоянии более двух длин волн все типы интегралов дают одинаковый результат.

Однако при уменьшении поперечных размеров светового поля (или его характерных деталей) до размеров порядка длины волны также необходимо учитывать векторный характер светового поля. В этом случае применяют векторный вариант интегралов РЗ [35, 2, 3, 7, 12, 27-29] или метод разложения по плоским волнам [36, 37, 4, 20, 29]. Хотя эти подходы являются аналогичными, они подразумевают использование различных алгоритмов расчётов, достоинства и недостатки которых в ближней зоне дифракции было бы полезно выяснить.

Т.к. аналитическое решение тестовой задачи (дифракция плоской волны на круглой апертуре) известно только на оси, в данной работе оценка работы рассматриваемых алгоритмов во внеосевой области выполняется по сравнению с результатами конечно-разностного временного метода FDTD, реализованного в программном продукте фирмы К-8ой.

При численной реализации интегралов РЗ сложно воспользоваться какой-либо симметрией световых полей, поэтому для упрощения расчётов или получения аналитических выражений приходится прибегать к различным аппроксимациям [12, 24, 37, 38]. Однако, как показано в данной работе, использование даже непараксиальных аппроксимаций в ближней зоне является некорректным, поэтому для расчётов был использован алгоритм, разработанный в [21] и обобщённый на векторный случай в [39].

При использовании метода разложения по плоским волнам возможно применение алгоритма быст-

рого преобразования Фурье (БПФ) [10, 17], что значительно сокращает время расчёта. Однако применение алгоритма БПФ имеет свои недостатки, связанные с фиксированной дискретностью сигналов на входе и выходе, а также возможностью вычислять только поперечные распределения. При необходимости получения продольных (вдоль оптической оси) распределений использовать алгоритм, основанный на БПФ, нецелесообразно.

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

Также разработка быстрых алгоритмов прямого вычисления дифракционного интеграла РЗ до сих пор является актуальной темой для многих исследователей [11, 13, 14, 16, 17, 19-21]. В частности, алгоритм, предложенный в [21], основан на кусочно-постоянном представлении входного поля, что позволяет часть операций выполнить аналитически. Однако апробации этого алгоритма на тестовых примерах в ближней зоне дифракции не проводилось.

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

1. Непараксиальная скалярная модель: сравнение алгоритмов

Интегральное преобразование Рэлея-Зоммер-фельда первого типа [30, 40] в декартовых координатах имеет следующий вид:

z e'u ( 1 ^

E(u, v, z) = -—j| Е^ y) — [ ik - - ) dxdy

e

l2

где Е0(х,у) - входное поле,

I = ^1 (и - х)2 + (у - у)2 + г2 , Х0 - область, в которой задано входное поле, к = 2р /1 - волновое число,

1 - длина волны.

Скалярный непараксиальный оператор распространения с использованием разложения по плоским волнам записывается следующим образом [41]:

Е(u,v,z) = j jF(X,h)exp(ikzVl4W )x x exp \jk (Xu + hv)] dX dh, F(X, h) = ir jj Eo (x, y) exp [-ik(Xx + hy)] dx dy,

где ^(Х,^) - спектр разложения входного поля по

плоским волнам, : 01 < у/Х2 + Л2 < °2 - область учитываемых пространственных частот. При 01 = 0, о2 = 1 рассматриваются только распространяющиеся волны, а при о1 = 1, 02 > 1 - только затухающие волны.

Выражения (1.1) и (1.2) можно свести друг к другу, однако в последнем случае возможно применение различных быстрых алгоритмов, в том числе БПФ [10, 17], что значительно сокращает время расчёта, несмотря на удвоенное по сравнению с (1.1) интегрирование.

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

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

В данном разделе на тестовом примере дифракции плоской волны на круглой апертуре проводится сравнение этого алгоритма с прямым вычислением дифракционного интеграла (1.1) на основе алгоритма, разработанного в [21].

При дифракции плоской волны на круглой апертуре в скалярном случае известно аналитическое решение для значений интенсивности вдоль оптической оси [11]: 2

I (0,0, z) = 1 +

2 2 r + z2

2 z

V

r2 + z2

-cos{k(^r02 + z2 - z)},

где г0 - радиус ограничивающей апертуры.

1.1. Быстрый алгоритм, учитывающий радиальную симметрию

Для плоской освещающей волны, ограниченной круглой диафрагмой радиуса г0, входное поле в (1.1) и (1.2) будет иметь следующий вид:

Eo( x, y) = Eo(r, j) = Eo(r) =

1, r £ r0, 0, r > r„.

где (г, ф) - полярные координаты во входной плоскости.

В более общем случае, когда входное поле является «вихревым» и представимо в виде:

Е0

выражение (1.2) можно упростить:

S

S

о

Е(р, 9, z) = -im+1k2 exp(im9) х

^0 '0

хj jE0(r Jm (ksr)r dr

хехр(¡кгу] 1 -о2)(кор)оdo,

где (р, 8) - полярные координаты в выходной плоскости, о - радиальная координата в частотной плоскости, о0 - радиус учитываемых пространственных частот.

При численной реализации по теореме Найквиста 00 определяется дискретизацией входного поля Дг :

S <

2Ar

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

Рассмотрим интеграл (1.2) в полярных координатах только в области затухающих волн:

Е(р,9,z) = jexp(-Wо2 -1 )х

i

Г 2p

хI j F(о,f)exp[ikpocos(9-f)]df

о do.

Проанализируем функции в (1.8), зависящие от полярного угла. Экспоненциальный множитель по модулю равен единице. Спектральная функция ^(о, ф) (при фиксированном ф) убывает не медленнее, чем 1/ о , иначе не будет выполняться равенство Парсеваля. Таким образом, интегрирование по углу даст функцию, которая не возрастает с ростом о, и для дальнейшего анализа её можно заменить на константу:

I =

j exp (-kWо2 -1 )о do = ——

1 (кг)2

Абсолютная погрешность при замене верхнего предела на конечное значение ог :

А =

+-

1

kz (kz)2 относительная погрешность:

exp

А

e = I = '

Величина е монотонно убывает с ростом ог (можно доказать, взяв производную). Сделав замену / = кго2 -1, получим функцию, не зависящую

от определённых значений 1 и г . Для нахождения допустимой границы отсечения нужно задаться определённой погрешностью е и решить уравнение (1.11). В частности, для е = 0,04 получаем / = 5 , т.е. выбор в качестве верхней границы частот:

о^ = (kz) +1

обеспечивает погрешность расчёта (I.8) не выше 5%.

1.2. Численное моделирование при дифракции плоской волны на круглой апертуре

Для сравнения численных результатов с аналитическим выражением (I.3) были выполнены расчёты осевой интенсивности плоской волны (I.4), прошедшей через круглое отверстие радиусом r0=101, на различных расстояниях с помощью интегрального преобразование РЗ (I.1) и разложения по плоским волнам (I.6) при m = 0.

Для получения значений на оптической оси или в продольном сечении применение алгоритма БПФ неэффективно, т.к. для этого потребуется вычислить столько поперечных сечений, сколько желательно получить значений вдоль оси распространения пучка.

При вычислении (I.1) был использован алгоритм быстрого расчёта [21] с дискретизацией на входе Ar = 0,11, а для выражения (I.6) (m = 0) выполнялось численное интегрирование с дискретизацией на входе Ar = 0,011 и при о0 = 3 .

На рис. 1 приведены сравнительные результаты вычисления осевой интенсивности в ближней зоне дифракции, из которых видно, что результаты обоих методов практически полностью совпадают для расстояний больших длины волны. При этом временные затраты для вычисления ста точек на графике (рис. 1 а) на основе выражения (I.6), учитывающего радиальную симметрию моделируемой оптической схемы, на порядок меньше, чем для расчётов по формуле (I.1): в первом случае понадобилось 1,6 с, а во втором - 30 с (расчёты проводились на персональном компьютере Intel Pentium4 CPU 2,40GHz, RAM 512 МЬ).

Однако на расстояниях меньших длины волны метод разложения по плоским волнам несколько отличается от аналитического выражения, особенно при z < 0,11. По формуле (I.12) для получения корректных результатов при z = 0,11 необходимо учитывать пространственные частоты вплоть до oz = 8 , а при z = 0,011 д. б. oz > 79, что резко увеличивает временные затраты, т. к. в этом случае необходимо также увеличивать дискретизацию входного поля.

Для получения результатов, приведённых на рис. 1 в, верхний предел учитываемых спектральных частот был принят о0 = 20 , что в совокупности

х

0

0

/

0

с уменьшением шага дискретизации привело к увеличению расчётного времени до 46 с и превысило время прямого расчёта (1.1) по алгоритму [21]. Заметим, что для получения корректных результатов в выражении (1.1) при г < 0,11 понадобилась на порядок меньшая степень дискретизации, чем при использовании формулы (1.6). Таким образом, в данном случае алгоритм, разработанный в [21], оказывается предпочтительным.

На рис. 1г показано внеосевое продольное распределение интенсивности при дифракции плоской волны на круглой апертуре в области г е [0,11, 101], х е [-101, 101], рассчитанное по формуле (1.6) за 14 с. На расстояниях, больших одной десятой длины волны, алгоритм разложения по плоским волнам, учитывающий радиальную симметрию задачи, полностью совпадает с аналитическим решением и значительно опережает по скорости алгоритм [21].

Также предлагаемый быстрый алгоритм (см. раздел 1.1) удобнее БПФ. Удобство заключается в том, что он позволяет с произвольной дискретизацией (не связанной с бинарной системой) получать распределения на любых поверхностях за малое время. Недостатком является требование наличия радиальной симметрии или вихревой угловой зависимости входного пучка.

Более детально внеосевое распределение интенсивности в поперечных направлению распространения плоскостях на различных расстояниях от апертуры показано в табл. 1. Среднеквадратичное отклонение в расчётах, получаемых двумя реализованными методами, составляет менее 5%.

Расчёт поперечной области размером 200x200 отсчётов для интеграла РЗ (I.1) занимает 73 с. Для сравнения при использовании алгоритма, описанного в [42], расчёт полей размером 100x100 отсчётов требует 179 с, а в работе [10] расчёт оценён в 2000 с (для восьмипроцессорного компьютера SGI Origin 2000, CPU 300-MHz, RAM 2048 Mb), при этом прямое численное интегрирование методом прямоугольников без какой-либо оптимизации займёт более недели [42].

При использовании алгоритма через разложение ПВ (I.6), учитывающего радиальную симметрию, время расчёта составляет 22 с. Применение алгоритма БПФ позволяет вычислять распределения в поперечных плоскостях за несколько секунд, но в продольной плоскости это займёт значительное время, пропорциональное числу отсчётов вдоль оптической оси.

Таким образом, каждый из методов и их реализаций обладает достоинствами и недостатками. Относительно медленный алгоритм прямого расчёта интеграла РЗ играет важную роль на очень близких расстояниях от апертуры, а также как универсальный по отношению к типу падающей волны и форме апертуры тестовый инструмент. Метод разложения ПВ, реализованный через БПФ, является самым быстрым, хотя и требующим много памяти [10], и также универсален. Радиальная реализация (I.6) характеризуется высокой скоростью расчёта, очень экономична в требованиях к свободной памяти и удобна для расчётов распределений на любых поверхностях, но круг решаемых задач ограничен.

а)

б)

0,04

0,08 1

Рис. 1. Сравнение распределения осевой интенсивности, рассчитанной с помощью (1.1) (пунктирная линия) и (1.6) (точечная линия), с аналитическим выражением (1.3) (черный цвет сплошная линия): (а) на отрезке г е [0,11, 101] при и0 = 3, (б) на отрезке г е [0,011,11] при и0 = 5 и (в) на отрезке г е [0,011,0,11] при и0 = 5 (точечная линия) и а0 = 20 (штрихпунктирная линия); (г) фокусирующее поведение круглой апертуры

Таблица 1. Распределение интенсивности в поперечных направлению распространения плоскостях на различных расстояниях от апертуры

Амплитуда

Фаза

Сечение интенсивности для РЗ (сплошная линия) и ПВ (точечная линия) методов

N

8 10 12 X

с<

5

N

8 10 12 "к

с<

6

N

10 12 X

2. Непараксиальная векторная модель

Интегральная теорема Рэлея-Зоммерфельда первого типа (1.1) в векторной форме записывается следующим образом [35, 2]:

г е'ы ( 1 ^

Ех (u, v, г) = — Ц Е0 х(х У) ( ¡к -11 dх d>',

г е'ы ( 1 ^

Еу

Ег (и,V, г) = -2- Л [Е0х(х,у)(и - х) +

2— Е

е'ы ( 1 ^

+Е0у(х у)(v - у)]—г [ ¡к - dхdУ,

где ЕСх(х,у) и Е0у(х,у) - комплексные амплитуды х- и у-компонент входного электрического поля, г-ком-понента предполагается нулевой за счёт выбора системы координат. При линейной х-поляризации у-компонента в этой модели всегда отсутствует.

Во всех этих интегралах присутствует выражение для сферической волны ехр(/'к/)//. Как правило, это не позволяет аналитически вычислить интегралы (1.13). В [43] есть разложение сферической волны в ряд по функциям Бесселя (выражение 8.533), но замена подынтегрального выражения на бесконечный ряд, включающий спецфункции, не облегчает численный расчёт, поэтому в данном разделе для прямого вычисления (1.13) используется обобщение алгоритма [21], описанное в [39].

2.1. Непараксиальная аппроксимация дифракционных интегралов РЗ Рассмотрим непараксиальную аппроксимацию интегралов (1.13), часто используемую при вычислениях [44, 45, 7]. Такая аппроксимация является более точной, чем хорошо известная параксиальная аппроксимация.

При подстановке в выражениях (1.13) в подынтегральные экспоненты

г2 - 2

I» Я + -

2Я

где К = Vи2 + V2 + г2 , г = ^X + у2 , а в остальные

1

к3

получаем:

„2

места I » К , а также пренебрегая слагаемым с 3

Ех (р, 9, г) = - |ехР | Г—

1К2

х

х и Е0х (г, ф)ехр

-'к

гр соБ(ф-9) К

dф> г dr,

-2 А

'г ехр ('кК) г0 | г I

Еу(р, аг)=-!ехр Iя" -шIх

х 1| Е0у (г, ф)ехр

к гр соБ(ф-9) ' К

dф> г dr,

' ехр ("кК) г0 | г Ег (р, 9, г) = ^2 | ехр I 'к —

1К

2 А

г

2К 1Х

ХI [ Е0х (г, ф) (р СОБ 9-г СОБ ф)-

0

+Е0у (г, ф) (р бш 9 - г бш ф)] х

х ехр

-'к

гр СОБ(ф-9) К

dфf г dr.

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

2'Я

I

XС9 ехр

ехр ['а соб(ф - 9)] dф =

= 2РХ С/" ехр('т,9)/т, (а),

9

то выражения

Ех (р, 9, г) = --

К2

ехр ('т9)х

'т+1 гк ехр ('кК) . . Еу (р, 9, г) =--^-1 ехр ('т9) х

К2

|ехр I" 2к А Е0.- «г)^ |1 г

Ег

2К

К 2 1 4 ' з I 2К1 0 0 у л т { К

ехр ('т9) |ехр { гк ^ Е0х (г) |е'9 ^ { "ТР ^ - ^ ^ {"Г )| г ^ -

,.2

'т+1кехр

К

К

2К

ехр ('т9) |ехр {'к ^| Е0у (г) \в?Зт+Х |"р 1

К I + ^т-1 { "Р » г 2* .

Из выражений (1.17) видно, что при учёте векторного характера электромагнитного поля для |т| = 1 возможно ненулевое значение интенсивности

на оптической оси, причём за счёт г-компоненты. Данный факт также отмечался ранее [46, 7].

2.2. Быстрый алгоритм расчёта распространения вихревых полей на основе разложения по плоским волнам Полученные в предыдущем разделе формулы (1.17) просты в реализации и обеспечивают высокую скорость вычислений, но они являются приближенными в соответствии с (1.14), и при уменьшении расстояния от входной плоскости, т. е. в ближней зоне дифракции, будут давать неверный результат.

Для получения более точных выражений рассмотрим векторный вариант оператора распространения, использующего разложение по плоским волнам [36, 37]:

| Ех (х, у, г) Е( х, у, г) = Еу (х, у, г) =

{Ег( л; ^ г) |

гг™* ч |Р(X,п)А = I1 М(Х,п) {р (X, п) 1ехр

х ехр['к (XX+пу)] dX dn,

1 -

где

| Р (X, п) А =

{ Ру (X, п)1

<11

Е0х

Е0 у

2„ у

ех

:р [-'к (XX + пу)] dх dy

I

- спектры тангенциальных компонент входного электрического поля, определённого в области Х0, и матрица преобразования:

9

0

0

0

0

0

0

X

M E (X, h) =

л/i -

В полярных координатах выражения (1.18)-(I.20) принимают следующий вид:

e(p, 9, z) =

1 о2 2p

= — j jР(о,f)exp[ikopcos(9-f)]x 1 J J

х exp

ikzyj 1 - о2

о do df,

Р(о, f) =

1 0 0 1 о cos f о sin f

ТГо2 VTo

' Fx

' Fx

х exp [-ikro cos (j-f)] r dr dj. Учитывая вихревой вид компонент входного поля:

( Fx (о, f)N

Fy

= 2pim exp(imf) х

4

X

Jm (kro) r dr,

0 l 0 y v ' J

выражения (I.21)-(I.23) также можно упростить:

2p I -2m

E(p, 9, z) = J iim exp(im9)х

о2 p _-,

х j exp ikz^J 1 -о2 Qm(kop,9)оdo

где

Qm (t, 9) =

Jm (t) 0

oCm

'P

lpy (o)J J0 IE0y (r) J

0

Jm

vr-o2

I

Г Px

l py (o).

Jm (kro)rdr

Ст 8) = 2[-'8^т+1 () - ^т-1 ()] ,

(Г, 8) = 1 [-''8(г) + е-'8Зт_х (/)] , t = кор .

2.3. Сравнение алгоритмов расчёта при дифракции плоской волны на круглой апертуре в случае линейной х-поляризации

При дифракции линейно-поляризованной вдоль оси х плоской волны на круглой апертуре выражения (1.17) сводятся к следующему виду (т=0):

¡гк ехр (¡кЯ) Ех (р, 8, г) =--Я^ х

хН (¡к ] 7 0(Я)г dг,

Еу (р, 8, г) = 0,

к ехр

Ez (p, 9, z) =

R

jexp Iik 2R ]х

х cos

j ГxJ+rJ.(f J

r dr.

В осевую интенсивность будет вносить вклад только х-компонента (на оптической оси в (1.14) Я = г):

ik exp(ikz) r°t I r i Ex (0,0, z) =--^—'- j exp | ik— | r dr =

= exp

( f r 2Л I ' f

exp

l l

2 z

-1

Аналогичный, но более точный результат получается при использовании выражений (1.24)-(1.26), которые в рассматриваемом случае сводятся к следующему виду:

Г Р(о) Л Г0 ( 1 Л { Г0 71(кОГ0) Л

х = Л0 1^(кго)гdг =

lP (о) J 01 0

E(p, 9, z) = j exp [ikW 1 -o / J

k о

0

J 0(k op) 0

Ю

^cos

Jj(k or0)do.

J

При вычислении продольной составляющей ЕТ в интегральном выражении (1.30) возникнет особенность при о = 1, которую можно устранить интегрированием по частям:

x

h

о 0

0

о

о

2

х

о

(

\

х

2

2

Ег (р, 8, г) =

2-Г0 1

Г 2 = ^ОТБ 8 (

7 Л

Э

Эо

ехр

81 ехр (¡кгл/Г

о1

-¡о

'Л-

31 (к ор) 31 (ког0) d о =

г

= ^соэ 8

С учётом

31 (к ор) 31 (ког0) d о =

-о г °2 _ э

ехр (¡кг>/1 -о2) 31 (к ор) 31 (к ог0) - — со8 81 ехр (1кг41 -о2)—[ 31 (к ор) 31 (ког0)] d о.

[31 (ах)31 (Ьх)] = а30 (ах)31 (Ьх) +

+Ь31 (ах) ,70 (Ьх)-

231

окончательно получаем:

Е2 (р,8,г) = -^со88х

ехр

-—со88|ехр(гкг>/1 -о2)х

г о

х[р^

к о

dо.

Получившаяся формула более громоздка, чем (1.30), но она не содержит особенности в подынтегральном выражении при о = 1, а неопределённость при о = 0 легко устраняется, так как 31(г) ~ г/2 при малых значениях г.

Выражение (1.31) получено для конкретного вида падающей волны, а именно плоской. Далее рассмотрим способ обхода особенности для любого входного поля.

Интеграл в (1.30) для продольной компоненты является несобственным, но сходящимся. Теоретически всегда возможно преобразовать сходящийся несобственный интеграл в собственный, но в нашем случае это делать неудобно, так как сетка дискретизации уже задана заранее. По этой же причине затруднительно выделить узкую окрестность данной окружности и в ней вычислить интеграл приближённо-аналитически. Однако есть возможность обойти потенциальное деление на нуль, не отказываясь от квадратурной формулы прямоугольников и без заметной потери точности.

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

1

(с учётом замены переменной). Вычисляя

интеграл

dх

^ Г

./Т J

dх

(х полагается ма-

лым) по формуле Ньютона-Лейбница, получаем значение, равное \/2х , а при использовании квадратурной формулы центральных прямоугольников этот интеграл будет равен

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

в формуле (1.30) умножается на

л/2 .

Как видно из (1.30), вклад в осевую интенсивность будет давать только х-компонента:

Ех (0,0, г) =

2-г0

о2 р ---

<|ехр ¡кг41 -о2 31(ког0)dо,

которая не вычисляется аналитически, но в разделе 1.1 было показано её совпадение с аналитическим решением (1.3) при о1 = 0 и достаточно большом о2.

Сравнение выражений осевой интенсивности, получаемых при использовании аппроксимации (1.28):

1х (0,0, г) = 2 - 2соб

Г г2 Л к-0-2г

и аналитического решения (1.3), позволяет сделать вывод, что аппроксимация (1.17) будет верна только при г02 << г2. В этом случае

'+£ - -

Г г2 1 +

, 2г

2г

и выражение (1.3) будет совпадать с (1.33).

Численное сравнение (1.28) и (1.32) (практически полное совпадение последнего с аналитическим решением показано в разделе 1.1) приведено на рис. 2а.

На рис. 2 показаны сравнительные результаты расчёта дифракции линейно-поляризованной плоской волны на круглой апертуре радиусом 101 в ближней зоне на основе выражений (1.17) и (1.25)-(1.26). Как видно из результатов моделирования, аппроксимацию (1.17) некорректно применять на расстояниях сравнимых с размером апертуры.

2

2

о

г

о

х

о

х

о=о

х

1

о

2

г

0

г

г

- г =

2.4. Дифракция плоской волны на круглой микроапертуре: сравнение с ГОТО

В разделе 1 было показано, что при радиусе апертуры 101 в скалярном случае рассмотренные алгоритмы РЗ и ПВ практически полностью совпадают. При уменьшении размеров апертуры до сравнимых с длиной волны необходимо учитывать векторный характер светового поля.

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

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

Параметры расчёта при использовании алгоритма РБТБ, реализованного в программном продукте Я-8ой: длина волны 1 = 532 нм, радиус, ограничи-

вающий плоскую волну на входе r0 = 21, размер расчётной области х е [-51, 51 ], y е [-51, 51 ], z е [0, 51 ]. Толщина PML - 1, шаг дискретизации по пространству: 1/15, шаг дискретизации по времени: 1 /(30c) (c - скорость света в вакууме), положение плоскости регистрации: z = 0,31 (»160 нм) и z = 41 (»2 мкм). Длительность моделирования на компьютере Intel(R) Celeron(R) CPU 3,06 Ггц, ОЗУ 512 Мб примерно 29 минут.

В табл. 2 приведены сравнительные результаты расчёта дифракции плоской волны при линейной x-поляризации с помощью метода FDTD, дифракционных интегралов РЗ (1.13) и разложения по плоским волнам (I.25)-(I.26).

Для алгоритма ПВ в соответствии с (I.12) на расстоянии z = 0,31 от апертуры учитывались спектральные частоты в радиусе s0 = 3 и применялся обход особенности на основе формулы Ньютона-Лейбница, а на расстоянии z = 41 от апертуры полагалось s0 = 1. В табл. 2 приведены сравнительные

графики и среднеквадратичные погрешности 8 сечений полученной интенсивности для каждого метода: точечной линией показан график, полученный с помощью дифракционных интегралов, а сплошной -полученный с помощью метода FDTD.

3 2 1

(в) 0

/ /

v / / / к\ \\ >

/ t

2 3 4 5 6 7

9 10 11 12 1

\

\ А

( >

(д) 0

3 4 5 6 7 8 9 10 11 12 1

Рис. 2. Сравнительные результаты расчёта дифракции линейно-поляризованной плоской волны на круглой апертуре: осевое распределение суммарной интенсивности, рассчитанное по

и суммарной

Таблица 2. Сравнительные результаты расчёта дифракции плоской волны при линейной х-поляризации с помощью метода РОТО и дифракционных интегралов на расстоянии г = 0,31 и г = 41

л н о о W о о ч

с

о

Сечение амплитуды х-компоненты, |Ех|

Сечение амплитуды

z-компоненты,

E

Сечение суммарной интенсивности, |Е|

ГО

о" II

N

со Ин

-2-1 0 1 2 X

5=3,4%, 5о=2,9%

J 1

-2-10 1 2 X 5=52,3%, 50=26,5%

-2-1 0 1 2 X

5=11,4%, 50=7,5%

<N

-2-1 0 1 2 X

5=3,4%, 50=2,9%

0,3 0,2 0,1 0

| -i ¡t 4 í

П ь

Я/

-2-10 1 2 X 5=33,1%, 50=17,3%

-2-1 0 1 2 X

5=6,9%, 50=4,9%

с<

4

N

со Ин

-2-1 0 1 2 X

5=4,4%, 50=2,5%

-2-10 1 2 X 5=13,6%, 50=7,4%

-2-1 0 1 2 X

5=9,1%, 50= 1,8%

<N

-2-1 0 1 2 X

5=4,4%, 50=2,5%

-2-1 0 1 2 X

5=13,1%, 50=6,5%

-2-1 0 1 2 X

5=8,9%, 50=1,8%

* расчёт на компьютере Intel Pentium4 CPU 2,40GHz, RAM 512 МЬ

Как следует из приведённых в табл. 2 результатов, основную погрешность в расчётах вносит именно продольная компонента. Результаты для х-компо-ненты (что соответствует скалярному случаю) совпадают для обоих алгоритмов и незначительно отличаются от предсказанных ГБТБ.

Погрешность в вычислении г-компоненты растёт с уменьшением расстояния до апертуры. При этом алгоритм ПВ (1.25)-(1.26) демонстрирует лучшие результаты, чем алгоритм РЗ [39], на расстояниях, меньших длины волны.

При удалении от апертуры всего на несколько длин волн результаты этих двух алгоритмов фактически совпадают (табл. 2 для г = 41) и отличаются от

полученных с помощью ГБТБ только масштабно, т. е. рассматриваемые в данной работе алгоритмы позволяют получать правильное распределение, но несколько завышенное по амплитуде. Среднеквадратичное отклонение приведённых к одному масштабу распределений 50 составляет в обоих случаях менее 2% для суммарной интенсивности.

Такая амплитудная «завышенность» может быть связана с тем, что рассмотренные методы РЗ и ПВ не подразумевают наличия у-компоненты в дифракционной картине изначально х-поляризованного поля. Метод же ГБТБ показывает, что эта компонента присутствует, хотя её доля энергии в суммарном электрическом поле невелика.

Во второй части статьи рассматривается модификация метода ПВ, позволяющая учесть наличие всех векторных компонент.

Заключение

На примере дифракции плоской волны на круглой апертуре в рамках непараксиальной скалярной модели показана эффективность разработанного в [21] алгоритма быстрого расчёта дифракционного интеграла РЗ в очень близкой к апертуре области (на расстоянии z < 0,11). Универсальный по отношению к типу падающей волны и форме апертуры относительно быстрый (ускорение по сравнению с [42] почти в 10 раз) прямой расчёт дифракционного интеграла полезен как тестовый инструмент в различных задачах.

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

Хотя самым быстрым и в то же время достаточно универсальным является метод, реализованный через БПФ, он требует много памяти. Кроме того, к недостаткам метода относятся фиксированная дискретность сигналов на входе и выходе, а также возможность вычислять только поперечные распределения. При необходимости получения продольных (вдоль оптической оси) распределений использовать алгоритм, основанный на БПФ, нецелесообразно.

При расчёте дифракции на микроапертурах также необходимо учитывать векторный характер светового поля. В работе было показано, что при численной реализации векторных интегралов РЗ использование даже непараксиальных аппроксимаций в ближней зоне является некорректным, поэтому для расчётов был использован алгоритм, описанный в [39].

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

Анализ результатов применения указанных выше методов при дифракции плоской волны на микроапертуре показал, что основную погрешность в расчётах вносит именно продольная компонента. Погрешность в вычислении z-компоненты растёт с уменьшением расстояния до апертуры.

Однако всё же на расстояниях меньших длины волны (но больше одной десятой длины волны) алгоритм ПВ (1.25)-(1.26) демонстрирует лучшие результаты как по точности, так и по скорости, чем алгоритм РЗ [39], хотя в этой области ресурсозатрат-ность алгоритма возрастает обратно пропорцио-

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

При удалении от апертуры всего на несколько длин волн результаты этих двух алгоритмов фактически совпадают и отличаются от полученных с помощью FDTD только масштабно, т.е. рассматриваемые в данной работе алгоритмы позволяют получать правильное распределение (среднеквадратичное отклонение с учётом масштаба составляет менее 2%), но несколько завышенное по амплитуде.

Такая амплитудная «завышенность» может быть связана с тем, что рассмотренные методы РЗ и ПВ не подразумевают наличия y-компоненты в дифракционной картине изначально х-поляризованного поля. Метод же FDTD показывает, что эта компонента присутствует.

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

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

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» (грант CRDF PG08-014-1), грантов РФФИ 10-07-00109-а, 10-07-00438-а и грантов Президента РФ поддержки ведущей научной школы НШ-7414.2010.9 и молодого кандидата наук (МК-64571.2010.2).

Литература

1. Martinez-Herrero, R. Vectorial structure of nonparaxial electromagnetic beams / R. Martinez-Herrero, P.M. Me-jias, S. Bosch, A. Carnicer // J. Opt. Soc. Am. A. - 2001. -Vol. 18. - P. 1678-1680.

2. Ciattoni, A. Vectorial analytical description of propagation of a highly nonparaxial beam / A. Ciattoni, B. Crosig-nani, and P. D. Porto // Opt. Commun. - 2002. - Vol. 202. - P. 17-20.

3. Guha, Sh. Description of light propagation through a circular aperture using nonparaxial vector diffraction theory / Shekhar Guha, Glen D. Gillen // Optics Express. - 2005. -Vol. 13, No. 5. - P. 1424-1447.

4. Guo, H. Vector plane wave spectrum of an arbitrary polarized electromagnetic wave / Hanming Guo, Jiabi Chen, Songlin Zhuang // Optics Express. - 2006. - Vol. 14, No. 6. - P. 2095-2100.

5. Deng, D. Analytical vectorial structure of radially polarized light beams / Dongmei Deng and Qi Guo // Optics Letters. - 2007. - Vol. 32, No. 18. - P. 2711-2713.

6. Anokhov, S.P. Plane wave diffraction by a perfectly transparent half-plane / Sergey P. Anokhov // J. Opt. Soc. Am. A. - 2007. - Vol. 24, No. 9. - P. 2493-2498.

7. Ковалёв, А.А. Непараксиальная векторная дифракция гауссового пучка на спиральной фазовой пластинке / А.А. Ковалёв, В.В. Котляр // Компьютерная оптика. -2007. - Том 31, № 4. - С. 19-22.

8. Wu, G. Analytical vectorial structure of hollow Gaussian beams in the far field / Guohua Wu, Qihong Lou, Jun Zhou // Optics Express. - 2008. - Vol. 16, No. 9. -P. 6417-6424.

9. Zhou, G. The analytical vectorial structure of a nonparax-ial Gaussian beam close to the source / Guoquan Zhou // Optics Express. - 2008. - Vol. 16, No. 6. - P. 3504-3514.

10. Delen, N. Verification and comparison of a fast Fourier transform-based full diffraction method for tilted and offset planes / Nuri Delen and Brian Hooker // Applied Optics. - 2001. - Vol. 40, No. 21. - P. 3525-3531.

11. Cooper, I.J. Numerical integration of diffraction integrals for a circular aperture / I.J. Cooper, C.J.R. Sheppard, M. Sharma // Optik. - 2002. - Vol. 113, No. 7. - P. 293298.

12. Duan, K. A comparison of the vectorial nonparaxial approach with Fresnel and Fraunhofer approximations / Kailiang Duan, Baida Lu // Optik. - 2004. - Vol. 115, No. 5. - P. 218-222.

13. Cooper, I.J. The numerical integration of fundamental diffraction integrals for converging polarized spherical waves using a two-dimensional form of Simpson's 1/3 Rule / I.J. Cooper, C.J.R. Sheppard, and M. Roy // Journal of Modern Optics. - 2005. - Vol. 52, No. 8. - P. 11231134.

14. Veerman, J.A.C. Calculation of the Rayleigh-Sommer-feld diffraction integral by exact integration of the fast oscillating factor / Jan A.C. Veerman, Jurgen J. Rusch, and H. Paul Urbach // J. Opt. Soc. Am. A. - 2005. - Vol. 22, No. 4. - P. 636-646.

15. Zhao, Zh. Focusing and diffraction by an optical lens and a small circular aperture / Zhiguo Zhao, Kailiang Duan, Baida Lu // Optik. - 2006. - V. 117. - P. 253-258.

16. Wang, X. Numerical calculation of a converging vector electromagnetic wave diffracted by an aperture by using Borgnis potentials. I. General theory / Xueen Wang, Zhaozhong Fan and Tiantong Tang // J. Opt. Soc. Am. A. - 2006. - Vol. 23, No. 4. - P. 872-877.

17. Shen, F. Fast-Fourier-transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula / Fabin Shen and Anbo Wang // Applied Optics. -2006. - Vol. 45, No. 6. - P. 1102-1110.

18. Kotlyar, V.V. Sharp focus area of radially-polarized Gaussian beam by propagation through an axicon / V.V. Kotlyar, A.A. Kovalev, S.S. Stafeev // Prog. In Electr. Res. C. - 2008. - Vol. 5. - P. 35-43.

19. Nascov, V. Fast computation algorithm for the Rayleigh-Sommerfeld diffraction formula using a type of scaled convolution / Victor Nascov and Petre Catalin Logofatu // Applied Optics. - 2009. - Vol. 48, No. 22. - P. 43104319.

20. Matsushima, K. Band-Limited Angular Spectrum Method for Numerical Simulation of Free-Space Propagation in Far and Near Fields / Kyoji Matsushima, Tomoyoshi Shi-mobaba // Optics Express. - 2009. - Vol. 17, No. 22. -P. 19662-19673.

21. Устинов, А. В. Быстрый способ вычисления интеграла Рэлея-Зоммерфельда первого типа // Компьютерная оптика. - 2009. - Т. 33, № 4. - С. 412-419

22. Osterberg, H. Closed Solutions of Rayleigh's Diffraction Integral for Axial Points / H. Osterberg and L.W. Smith // J. Opt. Soc. Am. - 1961. - Vol. 51. - P. 10501054.

23. Wolf, E. Comparison of the Kirchhoff and the Rayleigh-Sommerfeld Theories of Diffraction at an Aperture / E. Wolf, E.W. Marchand // J. Opt. Soc. Am. - 1964. -Vol. 54(5). - P. 587-594.

24. Gravelsaeter, T. Diffraction by circular apertures. 1: Method of linear phase and amplitude approximation / Tore Gravelsaeter and Jakob J. Stamnes // Applied Optics.

- 1982. - Vol. 21, No. 20. - P. 3644-3651.

25. Sheppard, C.J.R. Diffraction by a circular aperture: a generalization of Fresnel diffraction theory / C.J.R. Sheppard, M. Hrynevych // J. Opt. Soc. Am. A. - 1992. -Vol. 9, No. 2. - P. 274-281.

26. Mielenz, K.D. Optical Diffraction in Close Proximity to Plane Apertures. I. Boundary-Value Solutions for Circular Apertures and Slits / Klaus D. Mielenz // J. Res. Natl. Inst. Stand. Technol. - 2002. - Vol. 107. - P. 355-362.

27. Romero, J.A. Vectorial approach to Huygens's principle for plane waves: circular aperture and zone plates / Julio A. Romero and Luis Hernández // J. Opt. Soc. Am. A. -2006. - Vol. 23, No. 5. - P. 1141-1145.

28. Romero, J.A. Diffraction by a circular aperture: an application of the vectorial theory of Huygens's principle in the near field / J.A. Romero and L. Hernández // J. Opt. Soc. Am. A. - 2008. - Vol. 25, No. 8. - P. 2040-2043.

29. Li, J. The rigorous electromagnetic theory of the diffraction of vector beams by a circular aperture / Jianlong Li, Shifu Zhu, Baida Lu // Opt. Commun. - 2009. - Vol. 282.

- P. 4475-4480.

30. Born, M. Principles of Optics / M. Born, E. Wolf - 6th ed. - Pergamon, Oxford, 1980. - Chap. 8.3.

31. Andrews, C.L. Diffraction pattern in a circular aperture measured in the microwave region / C.L. Andrews // J. Appl. Phys. - 1950. - Vol. 21. - P. 761-767.

32. Silver, S. Microwave aperture antennas and diffraction theory / S. Silver // J. Opt. Soc. Am. - 1962. - Vol. 52. -P. 131-139.

33. Totzeck, M. Validity of the scalar Kirchhoff and Rayleigh-Sommerfeld diffraction theories in the near field of small phase objects / M. Totzeck // J. Opt. Soc. Am. A.

- 1991. - V. 8, No. 1. - P. 27-32.

34. Tsoy, V.I. The use of Kirchhoff approach for the calculation of the near field amplitudes of electromagnetic field / V.I. Tsoy, L.A. Melnikov // Optics Communications. -2005. - V. 256. - P. 1-9.

35. Luneburg, R.K. Mathematical Theory of Optics / R.K. Luneburg - University of California Press, Berkeley, California, 1966.

36. Carter, W.H. Electromagnetic Field of a Gaussian Beam with an Elliptical Cross Section / W.H. Carter // J. Opt. Soc. Am. A. - 1972. - Vol. 62, No. 10. - P. 1195-1201.

37. Agrawal, G.P. Gaussian beam propagation beyond the paraxial approximation / G.P. Agrawal, D.N. Pattanayak // J. Opt. Soc. Am. A. - 1979. - Vol. 69, No. 4. - P. 575578.

38. Marathay, A.S. On the usual approximation used in the Rayleigh-Sommerfeld diffraction theory / A.S. Marathay, J.F. McCalmont // J. Opt. Soc. Am. A. - 2004. - Vol. 21.

- P. 510-516.

39. Хонина, С.Н. Алгоритмы быстрого расчёта дифракции радиально-вихревых лазерных полей на микроапертуре / С.Н. Хонина, А.В. Устинов, С.Г. Волотов-ский, М.А. Ананьин // Известия Самарского научного центра РАН. - 2010. - № 12(3). - С. 15-25.

40. Goodman, J.W. Introduction to Fourier Optics / J.W. Goodman - McGraw-Hill, 1968. - Chap. 3.

41. Виноградова, М.Б. Теория волн / М.Б. Виноградова, О.В. Руденко, А.П. Сухоруков - М.: Наука. Главная редакция физико-математической литературы, 1979. -384 с.

42. Балалаев, С.А. Реализация быстрого алгоритма преобразования Кирхгофа на примере бесселевых пучков / С. А. Балалаев // Компьютерная оптика. - 2006. -№ 30. - С. 69-73.

43. Gradshteyn, S. Table of Integrals, Series, and Products / S. Gradshteyn and I.M. Ryzhik - Elsevier, 2007.

44. Zhang, Y. Vector propagation of radially polarized Gaussian beams diffracted by an axicon / Yaoju Zhang, Ling Wang, Chongwei Zheng // J. Opt. Soc. Am. A. - 2005. -Vol. 22, No. 11. - P. 2542-2542.

45. Mei, Z. Nonparaxial analysis of vectorial Laguerre-Bes-sel-Gaussian beams / Z. Mei and D. Zhao // Opt. Express.

- 2007. - Vol. 15. - P. 11942-11951.

46. Helseth, L.E. Optical vortices in focal regions / L.E. Hel-seth // Opt. Commun. - 2004. - Vol. 229. - P. 85-91.

References

1. Martinez-Herrero, R. Vectorial structure of nonparaxial electromagnetic beams / R. Martinez-Herrero, P.M. Me-jias, S. Bosch, A. Carnicer // J. Opt. Soc. Am. A. - 2001. -Vol. 18. - P. 1678-1680.

2. Ciattoni, A. Vectorial analytical description of propagation of a highly nonparaxial beam / A. Ciattoni, B. Crosig-nani, and P. D. Porto // Opt. Commun. - 2002. - Vol. 202.

- P. 17-20.

3. Guha, Sh. Description of light propagation through a circular aperture using nonparaxial vector diffraction theory / Shekhar Guha, Glen D. Gillen // Optics Express. - 2005. -Vol. 13, No. 5. - P. 1424-1447.

4. Guo, H. Vector plane wave spectrum of an arbitrary polarized electromagnetic wave / Hanming Guo, Jiabi Chen, Songlin Zhuang // Optics Express. - 2006. - Vol. 14, No. 6. - P. 2095-2100.

5. Deng, D. Analytical vectorial structure of radially polarized light beams / Dongmei Deng and Qi Guo // Optics Letters. - 2007. - Vol. 32, No. 18. - P. 2711-2713.

6. Anokhov, S.P. Plane wave diffraction by a perfectly transparent half-plane / Sergey P. Anokhov // J. Opt. Soc. Am. A. - 2007. - Vol. 24, No. 9. - P. 2493-2498.

7. Kovalev, A.A., Kotlyar, V.V. Nonparaxial vectorial diffraction of the Gaussian beam by a spiral phase plate // Computer Optics. - 2007. - Vol. 31, N 4. - P. 19-22.

8. Wu, G. Analytical vectorial structure of hollow Gaussian beams in the far field / Guohua Wu, Qihong Lou, Jun Zhou // Optics Express. - 2008. - Vol. 16, No. 9. -P. 6417-6424.

9. Zhou, G. The analytical vectorial structure of a nonparax-ial Gaussian beam close to the source / Guoquan Zhou // Optics Express. - 2008. - Vol. 16, No. 6. - P. 3504-3514.

10. Delen, N. Verification and comparison of a fast Fourier transform-based full diffraction method for tilted and offset planes / Nuri Delen and Brian Hooker // Applied Optics. - 2001. - Vol. 40, No. 21. - P. 3525-3531.

11. Cooper, I.J. Numerical integration of diffraction integrals for a circular aperture / I.J. Cooper, C.J.R. Sheppard, M. Sharma // Optik. - 2002. - Vol. 113, No. 7. - P. 293298.

12. Duan, K. A comparison of the vectorial nonparaxial approach with Fresnel and Fraunhofer approximations / Kailiang Duan, Baida Lu // Optik. - 2004. - Vol. 115, No. 5. - P. 218-222.

13. Cooper, I.J. The numerical integration of fundamental diffraction integrals for converging polarized spherical waves using a two-dimensional form of Simpson's 1/3 Rule / I.J. Cooper, C.J.R. Sheppard, and M. Roy // Journal of Modern Optics. - 2005. - Vol. 52, No. 8. - P. 11231134.

14. Veerman, J.A.C. Calculation of the Rayleigh-Sommer-feld diffraction integral by exact integration of the fast oscillating factor / Jan A.C. Veerman, Jurgen J. Rusch, and H. Paul Urbach // J. Opt. Soc. Am. A. - 2005. - Vol. 22, No. 4. - P. 636-646.

15. Zhao, Zh. Focusing and diffraction by an optical lens and a small circular aperture / Zhiguo Zhao, Kailiang Duan, Baida Lu // Optik. - 2006. - V. 117. - P. 253-258.

16. Wang, X. Numerical calculation of a converging vector electromagnetic wave diffracted by an aperture by using Borgnis potentials. I. General theory / Xueen Wang, Zhaozhong Fan and Tiantong Tang // J. Opt. Soc. Am. A.

- 2006. - Vol. 23, No. 4. - P. 872-877.

17. Shen, F. Fast-Fourier-transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula / Fabin Shen and Anbo Wang // Applied Optics. -2006. - Vol. 45, No. 6. - P. 1102-1110.

18. Kotlyar, V.V. Sharp focus area of radially-polarized Gaussian beam by propagation through an axicon / V.V. Kotlyar, A.A. Kovalev, S.S. Stafeev // Prog. In Electr. Res. C. - 2008. - Vol. 5. - P. 35-43.

19. Nascov, V. Fast computation algorithm for the Rayleigh-Sommerfeld diffraction formula using a type of scaled convolution / Victor Nascov and Petre Catalin Logofatu // Applied Optics. - 2009. - Vol. 48, No. 22. - P. 43104319.

20. Matsushima, K. Band-Limited Angular Spectrum Method for Numerical Simulation of Free-Space Propagation in Far and Near Fields, / Kyoji Matsushima, Tomoyoshi Shimobaba // Optics Express. - 2009. - Vol. 17, No. 22. -P. 19662-19673.

21. Ustinov, A.V. The fast way for calculation of first class Rayleigh-Sommerfeld integral // Computer Optics. - 2009.

- V. 33, N 4. - P. 412-419 [in Russian]

22. Osterberg, H. Closed Solutions of Rayleigh's Diffraction Integral for Axial Points / H. Osterberg and L.W. Smith // J. Opt. Soc. Am. - 1961. - Vol. 51. - P. 10501054.

23. Wolf, E. Comparison of the Kirchhoff and the Rayleigh-Sommerfeld Theories of Diffraction at an Aperture / E. Wolf, E.W. Marchand // J. Opt. Soc. Am. - 1964. -Vol. 54(5). - P. 587-594.

24. Gravelsaeter, T. Diffraction by circular apertures. 1: Method of linear phase and amplitude approximation / Tore Gravelsaeter and Jakob J. Stamnes // Applied Optics.

- 1982. - Vol. 21, No. 20. - P. 3644-3651.

25. Sheppard, C.J.R. Diffraction by a circular aperture: a generalization of Fresnel diffraction theory / C.J.R. Sheppard, M. Hrynevych // J. Opt. Soc. Am. A. - 1992. -Vol. 9, No. 2. - P. 274-281.

26. Mielenz, K.D. Optical Diffraction in Close Proximity to Plane Apertures. I. Boundary-Value Solutions for Circular Apertures and Slits / Klaus D. Mielenz // J. Res. Natl. Inst. Stand. Technol. - 2002. - Vol. 107. - P. 355-362.

27. Romero, J.A. Vectorial approach to Huygens's principle for plane waves: circular aperture and zone plates / Julio A. Romero and Luis Hernández // J. Opt. Soc. Am. A. -2006. - Vol. 23, No. 5. - P. 1141-1145.

28. Romero, J.A. Diffraction by a circular aperture: an application of the vectorial theory of Huygens's principle in the near field / J.A. Romero and L. Hernández // J. Opt. Soc. Am. A. - 2008. - Vol. 25, No. 8. - P. 2040-2043.

29. Li, J. The rigorous electromagnetic theory of the diffraction of vector beams by a circular aperture / Jianlong Li, Shifu Zhu, Baida Lu // Opt. Commun. - 2009. - Vol. 282.

- p. 4475-4480.

30. Born, M., Wolf, E. Principles of Optics, 6th ed. - Perga-mon, Oxford, 1980. - Chap. 8.3.

31. Andrews, C.L. Diffraction pattern in a circular aperture measured in the microwave region // J. Appl. Phys. -1950. - Vol. 21. - P. 761-767.

32. Silver, S. Microwave aperture antennas and diffraction theory // J. Opt. Soc. Am. - 1962. - Vol. 52. - P. 131-139.

33. Totzeck, M. Validity of the scalar Kirchhoff and Rayleigh-Sommerfeld diffraction theories in the near field of small phase objects / M. Totzeck // J. Opt. Soc. Am. A.

- 1991. - V. 8, No. 1. - P. 27-32.

34. Tsoy, V.I. The use of Kirchhoff approach for the calculation of the near field amplitudes of electromagnetic field / V.I. Tsoy, L.A. Melnikov // Optics Communications -2005. - V. 256. - P. 1-9.

35. Luneburg, R.K. Mathematical Theory of Optics - University of California Press, Berkeley, California, 1966.

36. Carter, W.H. Electromagnetic Field of a Gaussian Beam with an Elliptical Cross Section // J. Opt. Soc. Am. A. -1972. - Vol. 62, No. 10. - P. 1195-1201.

37. Agrawal, G.P. Gaussian beam propagation beyond the paraxial approximation / G.P. Agrawal, D.N. Pattanayak // J. Opt. Soc. Am. A. - 1979. - Vol. 69, No. 4. - P. 575-578.

38. Marathay, A.S. On the usual approximation used in the Rayleigh-Sommerfeld diffraction theory / A.S. Marathay, J.F. McCalmont // J. Opt. Soc. Am. A. - 2004. - Vol. 21.

- P. 510-516.

39. Khonina, S.N. Fast calculation algorithms for diffraction of radially-vortical laser fields on the microaperture / S.N. Khonina, A.V. Ustinov, S.G. Volotovsky, M.A. Ananin // Izvest. SNC RAS - 2010. - V. 12(3). - P. 1525. - (in Russian).

40. Goodman, J.W. Introduction to Fourier Optics -McGraw-Hill, 1968. - Chap. 3.

41. Vinogradova, M.B. Wave Theory, 2nd ed. / M.B. Vino-gradova, O.V. Rudenko, and A.P. Sukhorukov - Moscow, 'Nauka" Publisher, 1979. - (in Russian).

42. Balalayev, SA. Realisation of fast algorithm of Kirchhoff s diffraction integral on an example of Bessel modes / S.A. Balalayev, S.N. Khonina // Computer Optics.

- 2006. - V. 30. - P. 69-73. - (in Russian).

43. Gradshteyn, S. Table of Integrals, Series, and Products / S. Gradshteyn and I.M. Ryzhik - Elsevier, 2007.

44. Zhang, Y. Vector propagation of radially polarized Gaussian beams diffracted by an axicon / Yaoju Zhang, Ling Wang, Chongwei Zheng // J. Opt. Soc. Am. A. - 2005. -Vol. 22, No. 11. - P. 2542-2542.

45. Mei, Z. Nonparaxial analysis of vectorial Laguerre-Bes-sel-Gaussian beams / Z. Mei and D. Zhao // Opt. Express.

- 2007. - Vol. 15. - P. 11942-11951.

46. Helseth, L.E. Optical vortices in focal regions // Opt. Commun. - 2004. - Vol. 229. - P. 85-91.

PROPAGATION OF THE RADIALLY-LIMITED VORTICAL BEAM IN A NEAR ZONE.

PART I. CALCULATION ALGORITHMS

S.N. Khonina 12, A.V. Ustinov ', A.A. Kovalev 12, S.G.Volotovsky 1 'institution of Russian Academy of Sciences, Image Processing Systems Institute RAS, 2S.P. Korolyov Samara State Aerospace University

Abstract

On an example of plane wave diffraction by a circular aperture in a near zone (the order of several wavelengths) comparison of calculation algorithms such as vectorial Rayleigh-Sommerfeld diffraction integral (RS) and plane wave expansion (PWE) on accuracy and speed of calculations is executed.

In a scalar case that corresponds to calculation of a cross-section components of an electric field, results differ only in area very close to the aperture.

In a vector case at calculation of the longitudinal component in PWE method there is a singularity in the domain of spectral frequencies which radius is close to unit. Various variants of avoiding of this singularity are offered. On distance of several wavelengths results of two considered algorithms coincide and differ from finite-difference time domain (FDTD) method only in scale (the root-mean-square deviation with account of scale makes less than 2%). Thus, the algorithms considered in given work allow to receive for essentially smaller time structurally true (but a little overestimated on amplitude) a picture of diffraction in a near zone. Such "overestimation" of amplitude can be connected by that considered RS and PWE methods do not mean presence y-components in diffraction picture of initially x-polarized field. In the second part of the paper a modification of PWE method allowing to consider presence of all vector components is considered.

Key words: Rayleigh-Sommerfeld diffraction integral, plane wave expansion, diffraction by a circular aperture, vortical beam.

Сведения об авторах

Хонина Светлана Николаевна, доктор физико-математических наук, профессор Самарского государственного аэрокосмического университета имени академика С.П. Королёва; ведущий научный сотрудник Учреждения Российской академии наук Институт систем обработки изображений РАН. Область научных интересов: дифракционная оптика, сингулярная оптика, модовые и поляризационные преобразования, оптическое манипулирование, оптическая и цифровая обработка изображений. E-mail: khonina@smr.ru .

Svetlana Nikolaevna Khonina, Doctor of Physical and Mathematical Sciences; Professor of the Samara State Aerospace University named after S.P. Korolev. Leading researcher of the Image Processing Systems Institute of the RAS. Research interests: diffractive optics, singular optics, mode and polarization transformations, optical manipulating, optical and digital image processing.

Устинов Андрей Владимирович, 1968 года рождения, в 1991 году окончил Куйбышевский авиационный институт имени академика С.П. Королёва (КуАИ) по специальности «Прикладная математика», работает ведущим программистом в Учреждении Российской академии наук Институт систем обработки изображений РАН. Область научных интересов: разработка программ моделирования работы оптических элементов; обработка изображений, в частности, гидродинамических процессов и биомедицинских. E-mail: andr@smr.ru.

Andrey Vladimirovich Ustinov, (b. 1968) graduated from Kuibyshev Aviation Institute named after academician S. P. Korolev (KuAI) on a specialty "Applied mathematics", works as the leading programmer in the Image Processing Systems Institute of the RAS. Research interests: software design for modeling of optical elements operating; images processing, particularly images of hydrody-namic processes and biomedical images.

Ковалёв Алексей Андреевич, 1979 года рождения, в 2002 году окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва - СГАУ по специальности «Прикладная математика». Кандидат физико-математических наук (2005 год), работает научным сотрудником лаборатории лазерных измерений Института систем обработки изображений РАН (ИСОИ РАН), является докторантом кафедры технической кибернетики СГАУ. Ковалёв А.А. - специалист в области дифракционной оптики и нанофотоники. В списке научных работ 50 статей. Область научных интересов: математическая теория дифракции, сингулярная оптика, фотонно-кристаллические устройства. E-mail: alanko@smr.ru.

Alexey Andreevich Kovalev (b. 1979), graduated (2002) from the S.P. Korolyov Samara State Aerospace University (SSAU)), majoring in Applied Mathematics. He received his Candidate in Physics & Maths degree (2002). He is a researcher of Laser Measurements laboratory at the Image Processing Systems Institute of the Russian Academy of Sciences (IPSI RAS), holding a part-time position of assistant at SSAU's Technical Cybernetics sub-department. He is a specialist in such areas as diffractive optics and nanophotonics. He is co-author of 50 scientific papers. Research interests are mathematical diffraction theory, singular optics, and photonic crystal devices.

Волотовский Сергей Геннадьевич, 1959 года рождения, в 1984 году окончил Куйбышевский авиационный институт имени академика С.П. Королёва (КуАИ) по специальности «Прикладная математика», работает ведущим программистом в Учреждении Российской академии наук Институт систем обработки изображений РАН.

Область научных интересов: разработка программного обеспечения расчёта и моделирования работы элементов дифракционной оптики. E-mail: sv@smr.ru.

Sergey Gennadjevich Volotovsky (b. 1959) graduated from Kuibyshev Aviation Institute named after academician S. P. Korolev (KuAI) on a specialty "Applied mathematics", works as the leading programmer in the Image Processing Systems Institute of the RAS. Research interests: software design, modeling of systems with diffractive optical elements.

Поступила в редакцию 18 июля 2010 г.