Алгоритмы быстрого расчета дифракции радиально вихревых лазерных полей на микроапертуре Текст научной статьи по специальности «Физика»

УДК 541.144

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

© 2010 С.Н. Хонина1, А.В. Устинов1, С.Г. Волотовский1, М.А. Ананьин2

'Институт систем обработки изображений РАН, г. Самара 2Самарский государственный аэрокосмический университет

Поступила в редакцию 22.04.2010

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

Ключевые слова: непараксиальные дифракционные операторы распространения, фазовая вихревая сингулярность, ближняя зона дифракции

ВВЕДЕНИЕ

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

Непараксиальная модель, основанная на теории Релея-Зоммерфельда, позволяет получать согласующиеся с экспериментами результаты на очень близких расстояниях [22, 23]. Однако при уменьшении поперечных размеров светового поля (или его характерных деталей) до размеров порядка длины волны также необходимо учитывать векторный характер светового поля. В этом случае можно воспользоваться векторным вариантом интегралов Релея-Зоммерфельда [24, 2, 25] или методом разложением по плоским волнам [26, 27].

В первом случае сложно воспользоваться какой-либо симметрией световых полей, поэтому для упрощения расчетов или получения аналитических выражений приходится использовать аппроксимации [12, 28]. Во втором случае можно воспользоваться радиальной симметрией или разделением входного поля на радиальную и угловую составляющие и свести двумерные интегралы к одномерным [29, 30].

Хонина Светлана Николаевна, доктор физико-математических наук, ведущий научный сотрудник. E-mail: khonina@smr.ru.

Устинов Андрей Владимирович, ведущий программист, E-mail: andr@smr.ru.

Волотовский Сергей Геннадьевич, ведущий программист. E-mail: sv@smr.ru;

Ананьин Михаил Александрович, аспирант. E-mail: ananin@smr.ru.

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

Среди рассмотренных алгоритмов: непараксиальная аппроксимация дифракционных интегралов Релея-Зоммерфельда в векторной форме при радиально-вихревой симметрии; обобщение быстрого алгоритма вычисления интегрального преобразования Релея-Зоммерфельда для произвольных полей [21]; быстрый алгоритм, основанный на разложении по плоским волнам при радиально-вихревой симметрии, позволяющий применять преобразование Ханкеля [32] и реализация этого алгоритма для произвольных полей на основе быстрого преобразования Фурье (БПФ). Обсуждаются особенности реализации известных алгоритмов, а также достоинства, недостатки и область применимости предложенных.

1. НЕПАРАКСИАЛЬНАЯ СКАЛЯРНАЯ МОДЕЛЬ

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

Дифракционный интеграл Релея-Зоммер-фельда первого типа [33] в декартовых координатах имеет следующий вид:

Е{и, V, г) = -П Ц Ео у) { й - ^Мхф. (1)

где E0(x,y) - входное поле, 1 (ы-х)2 +(у-у)2 +22, Е0 - область, в которой задано входное поле, к = 2п / Л - волновое число, 1 - длина волны.

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

Ш(и,\, г) =Л ||Е)(х, у)х

Л

х| Цехр(¡к^ 1-|2 -г)2) ехр(-х) +г(У-yDdr|ldxdy ( )

Е0 (г, р) = ехр (-¡ка0г + ¡тр(,

(7)

где : ст1 +ц2 <а2 - область учитываемых пространственных частот. При а1 = 0, а2 = 1 рассматриваются только распространяющиеся волны, а при а1 = 1, а2 > 1 - только затухающие волны.

Вычисление выражения (2) может быть реализовано через БПФ [10, 17].

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

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

1, г < Г0

0, г > г

(3)

В общем случае, когда входное поле предста-вимо в радиально-вихревом виде:

Е0 (г, р) = Е0 (г) ехр(т р), (4)

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

и0

Е ( р,в, г ) = - ¡т+1 к ехр( ¡тв) | ехр (•>/1 -а2 )

|Е0(г)■1т (каг)г^ ■1т (кар)а<1

(5)

где г0 - радиус ограничивающей апертуры, а0 -радиус учитываемых пространственных частот, m -порядок фазовой вихревой сингулярности.

При численном интегрировании дискретизация входного поля А определяет (по теореме Найквиста) максимальное значение а0:

а0 < ^ • (6)

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

2. ДИФРАКЦИЯ РАДИАЛЬНО-ВИХРЕВОГО ПОЛЯ В РАМКАХ НЕПАРАКСИАЛЬНОЙ СКАЛЯРНОЙ МОДЕЛИ

Рассмотрим дифракцию на круглом микроотверстии конической волны с вихревой составляющей [35]:

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

В табл. 1 приведены результаты численного моделирования с помощью выражения (5) дифракции волны (7) на отверстии радиусом 10Л при а0 =0.5 и а0 =0.9 для осесимметричного (т=0) и вихревого (т=1) случаев. Размер апертуры выбирался из соображений исследовать область, где еще работает скалярная модель, но уже существенно влияние поляризации.

Параметры расчета: г0 = 10Л, шаг дискретизации Аг = 0,01Л, а0 = 3 , шаг дискретизации Аа = 0,01. Картины интенсивности в поперечных плоскостях показаны в области радиусом р = 2 Л.

По результатам моделирования, приведенным в табл. 1 видно, что при увеличении параметра конической волны а0 для осесимметричного случая (т=0) происходит уменьшение размера центрального светового пятна. При а0 =0.9 скалярная теория предсказывает преодоление дифракционного предела по полуспаду интенсивности (РШИМ=0,391), равный для линзы 0,511.

При наличии вихревой составляющей на оптической оси будет нулевое значение интенсивности.

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

3. НЕПАРАКСИАЛЬНАЯ ВЕКТОРНАЯ МОДЕЛЬ

3.1. Алгоритм реализации векторного варианта интегральной теоремы Рэлея-Зоммерфельда

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

Ех(и V г) = -"2ПЯЕ0х(ХУ) е~[к--]\ЫУ, (8а) Еу(и Vг)=-2п Л Е0у(х у) ^т {¡к -'1^dxdy, (8б)

Е (м,v, г) = 2-{{[Е0х(x, У)(и -х)+Е0у(x, УУУ-У^

2п

х [ ¡к - - I dxdy

(8в)

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

Алгоритм быстрого вычисления выражения (1) был описан в работе [21]. Его можно приме-

х

Таблица 1. Дифракция конической вихревой волны в рамках непараксиальной скалярной модели.

и

и

о Я

m Й

^ К

CD g

а0 =0.5

-А-

а0 =0.9

Н

ill V._

г / \

/ N \

1- \

4 6 10 12 14 16 18 20 22 24 26 28 30

§ «

m ю

К I I

н =3

z=rn, FWHM = 0,66А.

12 345 67 89

z=m

Z=2,5^, FWHM = 0,39A.

Z=2,5^

нить для вычисления выражений (8а) и (8б), однако для вычисления выражения (8в) требуется модификация, которая и рассмотрена в данном разделе.

В дальнейшем изложении учтём, что входная амплитуда задана не аналитическим выражением, а дискретным набором точек, и сделаем замену переменных р = и - х; д = V - у. Так как в пределах каждой ячейки разбиения входная амплитуда является постоянной, то вычисления проводятся по дискретным формулам:

ются в [21]), то можно считать функции h(p, q), ph(p, q), qh(p, q) в пределах ячейки постоянными. В этом случае упомянутые функции заменяются своими значениями в центре ячейки и выносятся за знак интеграла. При этом формулы (9) примут следующий вид:

z)="2лZEx[mn]'hpqO J ] exp[k/p2+q +z1\pq

mn p„ q

EyQcyz)=-kZEy[m>n]•h(p»,q) J j exp^i^p2+q +z l^dq

Ex(x,y,z)=--—ZE>x[m,n] J J exp|k/p2 +q2 +z2~li.p,q)dpdq

2n m n pm qn

kz Pm+1 qn+1 p _

Ey(x,y,z)=—•ZE0y[m>n] J J ®Ф| ifajp1 +q +z2 kpqyipdq

mn pm q

ik pm+1 ^n+l Г I-"1

Ez(x,y,z)=--ZE0x[m,n] J J exxpl ikjp2 +q +z \ph(p,q)dpdq-

2n mn 1

Ez (x У, z) ^Ж^' pm +U0y[m,n]'qn)h'(Pm,q„) ZJI mn

(10)

Pm+ Qt+l |- -■

<J J e^pM/p2 +q2 + £ Ipdq

ik pm+i qn+i г I-"I

—270y[m>n] J J epi^/p2 +q2 qipqypdq

(9)

q)=■

p2 +q2 + z2 Mjf +q +z2)3'

где т, п - номера дискретного набора точек, в которых задано входное поле.

Если размер ячейки достаточно мал, в частности А < 0,2 г (подробнее условия обсужда-

где (~Рт, дп) - центр ячейки.

Отсюда следует, что вычисления ничем не отличаются от скалярного варианта [21].

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

Обобщение алгоритма на векторный случай

Также как и в скалярном варианте делается переход к полярным координатам, хотя подын-

20

тегральная функция не является радиально симметричной. Тем не менее, по угловой переменной интеграл берётся аналитически, поэтому в результате придём к одномерному интегралу.

Нам необходимо вычислить интегралы следующего вида:

Р2 «2

Р2 «2

JJ f (Р2 + q2) pdpdq, JJ f (p2 + q2 )qdpdq

Р «i

pi «i

Pi = Um - * Р2 = Um+1 - ^

«i =v»- у q2 =

n+1

AB: p = p2 ^ A' B': r cos ф = p2 BC : q = q2 ^ B' C': r sin ф = q2 CD: p = p1 ^ C' D': r cos ф = p1 DA : q = q1 ^ D' A': r sin ф = q1

Рис. 1 является схематичным: выпуклость или вогнутость границ не анализировалась.

Преобразование интеграла выражается формулами:

Р2 «2

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

jj f (r)r2 cos фdrdф и Ц f (r)r2 sin фdrdф, в кото-

D D

рых область D - некоторый криволинейный четырёхугольник. Интеграл по переменной ф берётся аналитически, но, так как стороны четырёхугольника описываются разными уравнениями, исходный интеграл превращается в сумму интегралов по переменной r. Пределы интегрирования зависят от того, как расположен исходный прямоугольник. Есть четыре варианта:

1) прямоугольник не соприкасается с осями координат;

2) прямоугольник соприкасается с осью Oq;

3) прямоугольник соприкасается с осью Op;

4) прямоугольник соприкасается с обеими осями.

В отличие от скалярного случая [21] варианты 1 и 4 не разделяются на два подварианта. (Более точно, ход вычислений немного различается, но конечный результат получается одинаковым.) Не будем описывать все варианты, рассмотрим только вариант 1 для интеграла

jJ f (r)r2 cos фdrdф в случае, когда правая ниж-

D

няя вершина ближе левой верхней к началу координат. При вычислениях используются формулы sin(arcsin х) = х, sin(arccos х) = V1 - х2 ,

cos(arccos х) = х, cos(arcsin х) = ■>/1 - х2 .

При замене переменных p = r cosф; q = r si^ прямоугольник ABCD преобразуется в криволинейный четырёхугольник A'B'C'D' (см. рис.1). Преобразование границ описывается соотношениями:

J J f (Р2 + «2)Pdpdq =

«1

JJ f (r)r2 cosфdrdф = JJ (.) + JJ (.) + JJ (.)

Р1 «1

arcsin(q2l r)

JJ f (r )r2 cos фdrdф = J f (r )r 2dr J cos

□ rC arccos( p2/r)

rB I 2"

= J f (r )r 2( -A 1 - p~)dr

rC arccos( p1 r)

JJ f (r )r2 cos фdrdф = J f (r )r 2dr J cos

arccos( p2 / r)

= J. f (r )r - Р--J 1 - "p^)dr

rA ' '

rA arccos( plr)

JJ f (r )r2 cos фdrdф = J f (r )r 2dr J cos

arcsin(«11 r)

'A ¡ 2

J f (r)r2U1 - 4 - qq-)dr

Если теперь раскрыть скобки и объединить интегралы с одинаковой подынтегральной функцией, то получим, что

,г

___с

pi g3

Рис. 1. Пояснения к преобразованию границ области при замене переменных

C

r

D

Г3

---D

qi

В

Р

Ф

Л / (г)

п

- )/ (г)

г 2 008

) / (г )г ^ - рг *

1 -

Р2

'В 'А

•йг + ^2 ) /(г)г<$г - ) /(г)гйг

'2

) / (г)гйг ,

г1

г2

) / (г )г 2йг

)/(г )г V1 - О2*

(11)

(12)

(13)

В формулах (11)-(13) г1, г2 - полярные радиусы некоторых вершин прямоугольника, a - одно из чисел р1, р2, , причём гарантировано, что подкоренное выражение неотрицательно. Интеграл вида (5) появляется при соприкосновении области интегрирования с осью координат, что при переходе к полярной системе даёт постоянный предел интегрирования (0 или п /2 ).

Подставим в формулы (11)-(13) функцию

-V 1 Л

к/

г2 + г2

/ (г) еХр - ^ ^ 1 . 2 2\3/2

к (г + г ) у Интеграл (11) вычисляется в явном виде: он ра-

вен

е

V

где

2 У

2 , 2 + г

'2 '2 ) /(г)г2йг « г, •) /(г)гйг

= г • ¡

С е^

V 7

А ^

2 У

,(14)

где г, = (г1 + г2)/2 . Оно применимо, если выполняется неравенство

2

г2 - г1 2г

п

384

Аналогичное преобразование интеграла (13) имеет вид:

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

)/(г)г2М - -"2йг = )/(г)гл/г2 - а2 йг «

г1 г г1

С^ Л ,(16)

•) / (г)гйТ = ^

, 2 2 = * г - а • ¡

V ^2 у

которое применимо при выполнении неравенства

(17)

а'(г2 - г1)2. <.П"

(г,2 - а2)2

384

3.2. Аппроксимация интеграла Рэлея-Зоммерфельда для микроапертур

Хотя алгоритм, рассмотренный в предыдущем разделе, значительно уменьшает время расчета выражений (8) по сравнению с прямыми квадратурными формулами, затраты эти все же существенные (см. далее результаты моделирования).

Рассмотрим аппроксимацию интегралов (8) при условии, что расстояние до точки наблюдения достаточно велико по сравнению с размерами апертуры [7, 36]. Такая аппроксимация является более точной, чем хорошо известная параксиальная аппроксимация, рассматривающая лишь область возле оптической оси.

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

1« Я +

г2 - 2(их + уу) 2Я

(18)

■2 --2 + -2 ~=4~х2 --2

= к^г22 + г2 (см. описание скалярного варианта в [21]). Интегралы (12) и (13) вычисляются или по квадратурным формулам, или с применением асимптотических формул, аналогичных скалярному варианту.

Для интеграла (12) преобразование имеет вид:

где Я = 4и2 + V2 + г2 , г = ^X + у2 , а в остальные места 1« Я, а также пренебрегая слагаемым

1

с Я3' получаем:

Ех (р,в, г) = -

¡г ехр (¡кЯ) % С, г2 Л -^—- ехр ¡к— ;

ЛЯ2 { \ 2Яу

Гг

) Е0х (г,Р)ехР

¡кгР008(Р-в)

Я

(19а)

йр\гФ,

¡г ехр (¡кЯ) % С г2 Л

Е (рвг>=- —Ле 0ехр [ ¡к 2К>

¿.л

) Е0у (г,р)ехр

, гроо&(р-в)

-¡к-

Я

(19б)

йф\гйг,

г.

о

i exp (MR) 1

Ez ^ z) = "^^J eXP

Mil ^ 2R у

J |E0x (r,<)(pcos0-r cos<) + E0y (r, <) (psin 0 - r sin <) ] x

exp

-ik

rpcos(<-0) R

d<>rdr

(19в)

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

Z Cq eXP(imq<)

exp [ia cos(<-0)]d< =

= 2nZCqim exP(imqQ)Jmq (a).

q

Выражения (19) можно упростить:

(20)

im+1zk exp (ikR )

Ex (p, 0 z) =--^-- exp(im0) x

'0

<J exp

ik-r—

V 2R у

R

E0 x (r )Jm i dr,

(21а)

im+1zk exp (ikR )

Ey (p, 0z) =--—1 exp(im0) x

0

<J exp

ik r

V 2R у

R

E0 y (r ) Jm ( krPrdr,

(21б)

R

Ez(pA z) =-

im+1k exp (ik

R2

expim0) x

jppexp \ikJfExWco^+EyWs^J/m {kpPrdr -

A

r2 ^

2 Jexp I ik^rR EJf)

krp R

-e-0

krp R

■r1 (fr -

(21в)

-2Ь (mЬr) ^eJ. ГMJ ^

Под интенсивностью здесь и далее понимается сумма модулей в квадрате всех компонент электрического поля: |Е|2 = |Ех |2 +| Еу| + \Ег |2.

3.3. Быстрый алгоритм расчета распространения радиально-вихревых полей на основе разложения по плоским волнам

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

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

E( x, y, z) =

I Ex (x, y, z) ^ Ey (x, y, z) Ez (x, y, z)

= JJ M(£n)

\ Fx (4,п) } \ Fy (4,П

xexp

exp[ik(£x+ny)]d^dn

(22)

где

Г Fx Ы Л

i E0x (x, y) ^

Я2 iJV E0 y (x, y)

F JJ

exp[-ik(4x+ny)]dxdy (23)

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

Me (4,п) =

1 0

4

0

1 п

V1-(4+п) V1-(4+п2).

.(24)

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

1 o1 2п

E(p, 0, z) = — JJ P(o, ф) exp [ikop cos(0 - ф)] x

x exp

Из выражений (21) видно, что при учете векторного характера электромагнитного поля для \т = 1 возможно ненулевое значение интенсивности на оптической оси, причем за счет 7-ком-поненты. Данный факт также отмечался при рассмотрении высокоапертурных фокусирующих систем [31].

ikz V1-c

o d o d ф

, (25)

P (o, ф) =

1

0

o cos ф

0 1

o sin ф

л/1 -o2 л/1-

o

I Fx (o, ф) ^ Fy (o, ф)

у

(26)

o 0

exp[[racos^^)]r<5r5^, (27)

\Fy (цф)^ j {Цу (г,ф)^

Учитывая радиально-вихревой вид компонент входного поля:

ГЕх(а, ф^ ,

vF (а ф) j

=i ех

E0x (r) ^

о vЦy(r) J

JJkro)rdr =

= 2ro'm ехр(/тф)

ГЕх (а) ^

U (а) J

(28)

выражения (16)-(18) также можно упростить:

Е(р, 0, z) = | у | i2m exp(im0) х

а2 р _-,

J ехр ikz\l 1 -а2 Qm (£ар, 0) а d а'

(29)

где

Qm(t,0) =

Jm(t)

о

;Q(i,0)

о

Jm() а

Sm(, 0)

ч^У(а)у

,(30)

где Ст (г,в) = 2[ев/т+1(г) - е-вт_,(Г)] ,

^т (г, в) = 2[/1С) + е/^)], г = кар .

В дальней зоне дифракции можно получить аналитический вид выражений (29) методом стационарной фазы [37], но в ближней зоне их использовать некорректно.

В табл. 2 приведены сравнительные результаты численного моделирования с помощью векторных дифракционных интегралов Релея-Зом-мерфельда (8), аппроксимирующих формул (21) и разложения по плоским волнам с учетом ра-диально-вихревой симметрии (28)-(30) дифракции линейно-поляризованной волны (7) на отверстии радиусом 10Л при а0 =0.5 для осесим-метричного (т=0) и вихревого (т=1) случаев. Параметры расчета для выражений (8): шаг дискретизации на входе и выходе Ах = Ау = 0,1Л ; для выражений (12): шаг дискретизации на входе и выходе Аг = 0,01Л ; для выражений (19)-(21) шаг дискретизации в объектной плоскости Аг = 0,01Л , шаг дискретизации в спектральной области Аа = 0,01, а0 = 3 . Картины интенсивно-

сти в поперечных плоскостях показаны в области радиусом p = 21.

Время расчета приведено для персонального компьютера Pentium®4, CPU 2,40GHz.

Как видно из табл. 2, применение аппрокси-мационных формул (21) в ближней зоне дифракции некорректно, в то время как выражения (28)-(30) позволяют получить с меньшими временными затратами результат, практически совпадающий с вычислениями по формулам Ре-лея-Зоммерфельда.

Также видно, что учет поляризационных эффектов сказывается на картине дифракции: за счет вклада продольной компоненты теряется радиальная симметрия. При этом в случае отсутствия фазовой вихревой сингулярности (m=0) центральное световое пятно вытягивается (уширяется) вдоль оси поляризации, а при наличии фазовой вихревой сингулярности первого порядка (|m|=1) появляется ненулевое значение интенсивности на оптической оси. При увеличении параметра конической волны а0 вклад продольной компоненты будет расти, и описанные явления будут проявляться с большей очевидностью.

В табл. 3 приведены сравнительные результаты численного моделирования с помощью векторных дифракционных интегралов Релея-Зом-мерфельда (8) и разложения по плоским волнам с учетом радиально-вихревой симметрии (28)-(30) дифракции линейно-поляризованной волны (7) на отверстии радиусом 101 при а0 =0,9 для осесимметричного (m=0) и вихревого (m=1) случаев. Параметры расчета такие же, как в предыдущем эксперименте.

На рис. 2 показаны сравнительные графики радиальных сечений распределений суммарной интенсивности в поперечных плоскостях, соответствующих максимальному значению интенсивности на оптической оси. Среднеквадратичное отклонение полученных результатов при m=0 составило 7,4%, а при m=1 - 3,5%.

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

Аналогичные результаты получаются при реализации формул (22)-(24) с использованием алгоритма БПФ. В этом случае при меньших временных затратах (несколько секунд) требуются значительные ресурсы оперативной памяти (при рассмотренных параметрах понадоби-

Таблица 2. Дифракция линейно-поляризованной конической вихревой волны при а0 =0.5 в рамках непараксиальной векторной модели (интенсивность х-компоненты показана светло-серым цветом, 7-компонента - тонкой линией,

суммарная интенсивность - толстой)

л о

О

л о

о

оо сч

л о

о

Л о

о

Л ! О

О

00

Л о

о

Распределение интенсивности вдоль оптической оси

Интенсивность в плоскости 2=11Л, ре [0,2Я]

ч N

N —V si

н -ч

(\ \

Ф У \

2 4 6 8

2 0 22 24 2 6 28 30

V

\, /- \ V

4 \ 1 \ ч

2468

20 22 24 26 28 3

N

«й т —Л

25 30

=1 5=

к 1 ГГ" л-

17 V

5 10 15 20 25 30

-1

4 *

-/ £ £

V 1 0 \ /- -V

2 4 6 8 1 0 1 2 1 4 1 6 1 8 20 22 24 26 28

Время расчета

17 мин

2 мин

41 сек

17 мин

2,5 мин

43 сек

50

14

лось более 400 Мб). Этот недостаток компенсируется универсальностью, т.е. возможностью применения данного алгоритма для задач, не обладающих радиальной симметрией.

В табл. 4 приведены сравнительные результаты численного моделирования с помощью векторных дифракционных интегралов Релея-Зом-мерфельда (8) и разложения по плоским волнам

(22)-(24) с использованием БПФ дифракции радиально-поляризованной волны (7) на отверстии радиусом 10Л при а0 =0,9 для осесиммет-ричного (т=0) и вихревого (т=1) случаев.

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

Таблица 3. Дифракция линейно-поляризованной конической вихревой волны при а0 =0,9 в рамках непараксиальной векторной модели, рассчитанная по выражениям (8) (тонкая линия) и (28)-(30) (толстая линия)

Интенсивность на

Распределение в плоскости 20, ре [0,2Л]

450 400 350 300 250 200 150 100 50

300 250

200

150

100

50

(а)

(б)

Рис. 2. Сравнительные графики радиальных сечений распределений суммарной интенсивности в поперечных плоскостях на расстоянии 70=2,6Л для т=0 (а) и расстоянии 70=2,25Л для т=1 (б): тонкая линия соответствует расчетам с помощью выражения (8), а толстая - с помощью (28)-(30)

Таблица 4. Дифракция радиально-поляризованной конической вихревой волны при а0 =0,9 в рамках непараксиальной векторной модели, рассчитанная по выражениям (8) и с использованием БПФ в формулах (22)-(24)

о

N

Распределение в плоскости 20, ре [0,2Л]

Е

аГё Ех

Е

Е, а^ Е2

ГЖНМ = 0,39Л

ГШЫМ = 0,41 Л

_1 ю II

FWЫM (-) = 0,32Л

сивности на оптической оси. Среднеквадратичное отклонение полученных результатов при т=0 составило 4,9%, а при т=1 - 7,5%.

Заметим, что в отличие от скалярного случая, рассмотренного в [10, 17, 19, 20], в векторном варианте разложения по плоским волнам

2

0

2

(а) (б)

Рис. 3. Сравнительные графики радиальных сечений распределений суммарной интенсивности в поперечных плоскостях на расстоянии 70=2,151 для т=0 (а) и расстоянии 70=2,61 для т=1 (б): тонкая линия соответствует расчетам с помощью выражения (8); толстая - с использованием БПФ в формулах (22)-(24).

имеется особенность при вычислении продольной компоненты. Как следует из поляризационной матрицы (24) и соответствующей матрицы (26) при достижении спектральных частот радиуса, близкого к единице + п2 = о2 = 1, что соответствует для конической волны а0 = 1 предельному наклону лучей, возникает особенность. Данная особенность не соответствует входному распределению, т.к. приводит к резкому росту продольной компоненты при 7=0, поэтому в разработанных алгоритмах в спектральной плоскости игнорируется узкая кольцевая область, средний радиус которой о2 = 1, а ширина Ао = 0,01.

ЗАКЛЮЧЕНИЕ

В данной работе рассмотрены:

- обобщение быстрого алгоритма вычисления интегрального преобразования Релея-Зоммер-фельда первого типа в векторной форме;

- непараксиальная аппроксимация дифракционных интегралов Релея-Зоммерфельда при радиально-вихревой симметрии;

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

- особенности реализации векторного оператора распространения на основе разложения по плоским волнам.

Проведено сравнение работы рассмотренных алгоритмов на примере моделирования дифракции вихревой конической волны на микроапертуре в ближней зоне. Показана эффективность в ближней зоне дифракции предложенного алгоритма непараксиального векторного распространения на основе разложения по плоским волнам с учетом радиально-вихревой симметрии по точности и экономии вычислительных ресурсов: при погрешности 4-8% обеспечивается сокращение времени вычислений в 7 раз.

Если в задаче не имеется радиальной симметрии, целесообразно использовать метод разложения по плоским волнам, реализованный через БПФ. Векторный вариант данного метода име-

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

БЛАГОДАРНОСТИ

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

СПИСОК ЛИТЕРАТУРЫ

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

2. CiattoniA, CrosignaniB, andPorto P.D. Vectorial analytical description of propagation of a highly nonparaxial beam // Opt. Commun. 2002. Vol. 202. Pp. 17-20 .

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

4. Hanming Guo, Jiabi Chen, Songlin Zhuang. Vector plane wave spectrum of an arbitrary polarized electromagnetic wave // OPTICS EXPRESS. 2006. Vol. 14, No. 6.

5. Dongmei Deng and Qi Guo. Analytical vectorial structure of radially polarized light beams// OPTICS LETTERS. 2007. Vol. 32. No.18. Pp. 2711-2713.

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

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

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

9. Guoquan Zhou. The analytical vectorial structure of a nonparaxial Gaussian beam close to the source, OPTICS EXPRESS. 2008. Vol. 16. No. 6. Pp. 3504-3514.

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

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

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

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

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

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

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

17. Fabin Shen and Anbo Wang. Fast-Fourier-transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula, APPLIED OPTICS// Vol. 45. No. 6. Pp. 1102-1110.

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

19. VictorNascov and Petre Catalin Logofatu. Fast computation algorithm for the Rayleigh-Sommerfeld diffraction formula using a type of scaled convolution// APPLIED OPTICS. 2009. Vol. 48. No. 22. Pp. 4310-4319.

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

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

22. Totzeck M. Validity of the scalar Kirchhoff and Rayleigh-Sommerfeld diffraction theories in the near field of small phase objects // J. Opt. Soc. Am. A. 1991. V. 8. No. 1. Pp. 27-32.

23. Tsoy V.I., Melnikov L.A. The use of Kirchhof approach for the calculation of the near feld amplitudes of electromagnetic feld // Optics Communications. 2005. V. 256. Pp. 1-9.

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

25. Lu B.D., Duan K.L. Nonparaxial propagation of vectorial Gaussian beams diffracted at a circular aperture// Opt. Lett. 2003. No. 28. Pp. 2440-2442.

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

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

28. Ковалев А.А., Котляр В.В. Дифракция плоской волны на ограниченной спиральной фазовой пластинке: параксиальная векторная теория // Компьютерная оптика. 2007. T. 31. №2, Pp. 4-8.

29. Khonina S.N., Kotlyar V.V., Soifer V.A., Shinkaryev M.V., Uspleniev G.V. /Trochoson// Optics Communications 1992. No. 91 (3-4), Pp.158-162 .

30. Павельев В.С., Хонина С.Н. Быстрый итерационный расчет фазовых формирователей мод Гаусса-Лагер-ра // Компьютерная оптика. 1997. Вып. 17. С. 15-20.

31. Хонина С.Н., Волотовский С.Г. Исследование применения аксиконов в высокоапертурной фокусирующей системе // Компьютерная оптика. 2010. Т. 34. №1. С. 35-51.

32. Khonina S.N., Kotlyar V.V., Soifer V.A. Fast Hankel transform for focusator synthesis // Optik. 1991. No. 88 (4). Pp. 182-184.

33. Goodman J.W. Introduction to Fourier Optics (McGraw-Hill, 1968), Chap. 3.

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

35. Victor V. Kotlyar, Alexey A. Kovalev, Svetlana N. Khonina, Roman V. Skidanov, Victor A. Soifer, Henna Elfstrom, Noora Tossavainen, andJari Turunen. Diffraction of conic and Gaussian beams by a spiral phase plate// APPLIED OPTICS. 2006. Vol. 45. No.12. Pp. 2656-2665.

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

37. Guoquan Zhou. Analytical vectorial structure of controllable dark-hollow beams in the far field //J. Opt. Soc. Am. A. 2009. Vol. 26, No. 7. 1654-1660.

FAST CALCULATION ALGORITHMS FOR DIFFRACTION OF RADIALLY-VORTICAL LASER FIELDS ON THE MICROAPERTURE

© 2010 S.N. Khonina1, A.V. Ustinov1, S.G. Volotovsky1, M.A. Ananin2

1 Image Processing Systems Institute of Russian Academy of Sciences, Samara 2 Samara State Aerospace University

Algorithms based on Rayleigh-Sommerfeld diffractive integrals of the first type and plane wave decomposition for nonparaxial vector simulation of propagation of the laser fields limited by the aperture are considered. Comparison of the considered algorithms on an example of diffraction of a vortical conic wave on the microaperture in a near zone is accomplished. Efficiency of the offered nonparaxial vector algorithm based on plane wave decomposition taking into account radially-vortical symmetry is shown on accuracy and economy of computing resources. Keywords: nonparaxial diffractive propagation operators, phase vortical singularity, a near zone of diffraction

Svetlana Khonina, Doctor of Physics and Mathematics, Leading Research Fellow. E-mail: khonina@smr.ru. Andrey Ustinov, Leading Programmer. E-mail: andr@smr.ru. Sergey Volotovsky, Leading Programmer. E-mail:sv@smr.ru. Michael Ananin, Graduate Student. E-mail: ananin@smr.ru