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

КОМПЬЮТЕРНАЯ ОПТИКА И ОБРАБОТКА ИЗОБРАЖЕНИЙ

УДК 621.38

СИНТЕЗ ВОЛНОВЫХ ПОЛЕЙ В РАМКАХ ЭЛЕКТРОМАГНИТНОЙ ТЕОРИИ

© 1999 П.Г. Серафимович, С.И. Харитонов Институт систем обработки изображений РАН, г. Самара

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

1. Алгоритм вычисления электромагнитного поля

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

Систему уравнений Максвелла запишем в следующем виде:

* ,ЕХ = - * х к

* -Е> = -к * у

{д хИу - * уИх )^+ ікНу

(дхНу - *уНх )1-ікНх

*гНх = І- д х (дхЕу - *уЕх ) - і-еЕу д,Ну = -і- д у (дхЕу - *уЕх ) + і-еЕх

(1)

Данное представление удобно тем, что в нем используется только 2 компоненты векторов электрического и магнитного поля. Представим вектора электрического и магнитного поля и функцию диэлектрической проницаемости в виде разложения в интеграл Фурье по поперечным координатам[1-2]:

да

Е (х, у, г) = | Е (а, р, г )ехр (гк (а х + р у ))(1а<1р

—да

да

Н (х, у, г) = | Н (а, р, г )ехр (гк (а х + р у ))йа<1р

—да

да

е (х, у, г) = | е(а, р, г )ехр (гк (ах + р у ))йа<1р (2)

—да

да

е—1(х, у, г) = | е—1(а, р, г )ехр (гк (ах + р у ))(1а<1р

—да

Здесь символ е—1 надо рассматривать,

- —1 1

как единый символ т.е. е ф —.

е

Функции Е(а, р, 2) и Н(а, р, 2) в дальнейшем будем называть поперечными пространственно-частотными компонентами.

Подставляя эти разложения в уравнения Максвелла, записанные в форме (2), получаем уравнения Максвелла в импульсном (поперечном пространственно-частотном) представлении:

ГЩНу (Щ ,юу, 2) -

Я- ЩуНх (Щ ,Щу > 2)у

*А (а, р, г) = -іка |

-к

X е -1 (а - щр - Щу \іюх(1юу + ікНу (ЩХ Щу,г)

а хИу (а х, а у, 2) -

- ауИх (ах, ау, 2)у

дгЕу (а, р, г) = —1кр |

—да\

х е —1 (а — а х, р — а у )йа хйа у + 1кИ х (а х, а у, г) дИ х (а, р, г) = 1к а (аЕу (а, р, г) - рЕх (а, р, г))-

да

- 1к | е(а - ах, р - ау)Еу (ах, ау, г)йахйау

—да

(3)

дИу (а, р, г) = 1кр(аЕу (а, р, г)—рЕх (а, р, г))+

да

+й | е(а—ах, р—ау )Ех (ах ,ау, г)йахйау

-да

Приведенные формулы имеют громоздкий вид. Для дальнейшего анализа приведенных уравнений запишем их в операторноматричной форме. Операторная форма записи уравнений позволит легко производить требуемые преобразования. В операторной записи система уравнений Максвелла имеет вид:

д2 w = , (4)

где W - матрица-столбец из 4-х компо-

г Е, ^

нент,

W =

Нх

V Ну У

, И-блочно-матричный диф-

ференциальный оператор: И

і (д. 0л

А =

( А 0 Л ( 0 Л

V 0 Ву V і 0у

(5)

і 0 Л (в-1 0 Л (- ду д. ^ (( 0 1^^

у +їк

к V 0 дУ V V 0 е-1^ V- ду дх ) V-1 0,

В =

(д. 0^ 0 д

(- д, д.^ (е 0^(0 -1

'у /V- ду дх У

+їк

е 0

0 е/

VI 0 у

го

W(x, у.г) = | W(a, Д z)exp(k(ax + ру^ойр.

(6)

Система уравнений в пространственночастотном представлении имеет вид:

д2 W (а, р, 2) = Н(а, р, 2) W (а, р, 2),

И(о, р, г) =

Ґ 0 1л

А(о ,р, г) 0

V 0 В(о, р, г)У V 1 0/

(7)

где двумерные матричные операторы действуют на двумерные матрицы-столбцы по формулам

А(о, р) " 9" = ік

у,

ґо 0Л

0 р

го

| е -1(о - ах , Р - ау , 2)

—го

( 9 (ах ,ау,2) л

V у(аху, 2)

ГЮу -

Vа у

, (8)

( 0 1' ^ 9''

йайа у + ік

ху V- 1 0, У У

В(о, Р)

9Л

VУy

= їк

ґо 0Л

Р -оІ9Л

V0 ркр -оку)

+ їк | Є(0 -ах , Р -ау ) 1

0 -1|( 9Іах ,ау )Л

0

/V

У(ах ,ау )у

сіах(іаУ

х у

Представление (2) в виде разложения по поперечным пространственно-частотным компонентам в матричной форме имеет вид:

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

дЕ (о, р, 2) = -іко(оИу (о, р, 2) - р Их (о, р, 2))+ +ікИу (о, р, 2)

д2Еу (о, р, 2) = -ікр(оИу ^ р, 2) - рИх (о, р, 2))- ікИх (о, р, 2)

д Дх (а, р, 2) = іко{аЕу (о, р, 2) - рЕх (а, р, 2) )-- ікЕу (о,р, 2)

д2Иу (о, р, 2) = ікр(оЕу (о, р, 2) - рЕх (о, р, 2))+

+ ікЕх (о, р, 2)

(9)

Общее решение данной системы уравнений имеет вид:

W (о , р , 2 ) =

х ехр

1 - о

Е + (о , р ) W + е (о , р) +^| + И + (о , р) W +11 (о , р) ~ГГрТ2 )+

У

х ехр

Е - (о , р) W -е (о , р) + ^ + И - (о , р) W -1 (о , р)

(-ік 4-

1 - о

2 - р 2 2

)

м W *• (о, р) =

+ д/1 - о 2 - р 2 о + УІ1 - о 2 - р 2 р р

- о

- р ' о

+ УІ1 - о2 - р 2 о + д/1 - о 2 - р 2 р

(10)

где М = ,/(2

- о -

р2)(о2 + р2) •

о

р

л/е г ,

м (о, р, г) = WSP,

^Е+(о, р)Л

Н + (о, р) Е - (о, р) V Н (о, р)

S =

ехр(ік^/1 - о2 - р2

ехр(- ік-)] 1 - о2 - р2

'2 У

W=(^+е (а,р) Ш+к (а,р) (а,р) (а,р))

Г1 0^

Е =

I о 1У.

Связь между пространственно-частотными компонентами в различных плоскостях имеет вид:

Ж(а, р, 2) = WSM -1WW(а, р,о), (12)

где W -матрица эрмитово-сопряженная

к W,

М =

V м1 е2

м1 = м+ем-е =

(о2 + р2)2

м2

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

р

Используя вышеприведенные выражения, подставляя в (6), получим выражение оператора распространения электромагнитного поля в свободном пространстве (про-пагатора):

к2

Ж(х, у, 2) = ГЖ(х, у,0) =—2 х 4р

ио

$Ж(х©, у©,0)х

г Ґ

хехр

ік

(11)

Перепишем полученное выражение в блочно-матричной форме

V V

о(х - х©)

р(у - у©)

\\

dx©dy©

dоdр

(10)

(13)

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

2. Определение псевдоскалярного произведения и унитарность оператора распространения

Введем в пространстве четырехкомпонентных векторов функций W понятие скалярного произведения. Введем скалярное произведение таким образом, чтобы оно не зависело от координаты z, и в тоже время произведение вектора самого на себя было пропорционально потоку вектора Умова-Пойтинга. Для этого запишем уравнения Максвелла в обычной форме для поля Е1 Н1 и для комплексно сопряженного поля

Е*2 Н*2:

гоїЕ1 = ікН1 гоїН1 = -ікЕ1 гоїЕ 2 = -ікН 2 гоїН 2 = ікЕ 2

(14)

вычитаем второе уравнение из первого, и, используя известную формулу векторного анализа

И* 2 го1Ех — Ехго1Н*2 = 1кНх Н*2 — 1кЕх Е *2, Е*2гоШ1 — Н1го1 Е*2 = —1кЕхЕ*2 + гкН1 Н*2,

Шу[а,Ь] = Ьго1а - аго№, (15)

получаем следующее выражение:

(Ну Ех Н*2 — (Ну Нх Е*2 = (16)

или (у([ Е1, Н *2 ] + [ Е *2, Н1]) = 0.

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

Я

ґ[Е1(х, у, 2Х\ И *2(x, у, 21)] +Л

Я

+ [Е 2(х, у,2Х\ И 1(x, у,21)]

' *

[Е1(х, у, 22ІИ 2(х, у, г2)] +

+ [ Е *2( х, у, 2 2\И 1( x, у, 2 2)] /

dxdy =

Л . (17)

dxdy =

Запишем последнее выражение в следующем виде:

Ц Ат ((х, у, г1))ОВ*(х, у, г1)йхйу =

|| Ат ((х, у, г 2)) О В *( х, у, г2) йхйу , (18)

где

о =

( 0 0 0 1>

0 0 -1 0

0 -1 0 0

V1 0 0 0у

или

||< А(х,у, г2),В(х,у, г2) > йхйу = ||< А(х,у, гх),В(х,у, гх) > йхйу , (19)

где < А,В >= (Ат) ОВ - псевдоскалярное произведение в пространстве четырехкомпонентных матриц столбцов.

Формула (19) означает, что оператор распространения (13) сохраняет скалярное произведение, т.е. является в пространстве с данным скалярным произведением унитарным оператором. Это свойство сохранения скалярного произведения можно использовать для решения обратных задач дифракции. Это в значительной мере упрощает вычисление градиента функционала невязки по аналогии с тем, как это наблюдается в случае скалярного приближения [3-4]. Кроме того, наличие сохраняющейся величины можно использовать для контроля правильности решения прямой задачи дифракции.

3. Пример численного моделирования

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

м(х У,0) = 2^11о(х у) ехр(і9(у))

0

0

V1 у

, (20)

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

б)

а)

г) в)

Рис. 1. Распределение интенсивности от геометрооптических фокусаторов

метры: размер ДОЭ по горизонтали и вертикали, соответственно: О х = 256 1 ,

Оу = 2561 , параметр гауссова пучка

а = 641, фокусное расстояние / = 5121.

Кроме того, для Рис. 1: (а) шаг дискретизации Д = 1, радиус кольца г0 = 321; (б) шаг дискретизации Д = 41, радиус кольца г0 = 321; (в) шаг дискретизации Д = 1, размер квадрата а х Ь = 321 х 321; (г) шаг дискретизации Д = 41, размер квадрата

а х Ь = 321 х 321.

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

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

шаг дискретизации на элементе Д = 1 , размеры апертуры элемента О х = 64 1 ,

Оу = 641, фокусное расстояние / = 641, размер квадрата в плоскости фокусировки, а х Ь = 161 х 161 . Кроме того, для рис.2: (а)

а = 161 - параметр гауссова пучка, дифракционная эффективность еГГ = 93%, среднеквадратичное отклонение рассчитанного распределения интенсивности от заданного гш8 = 20% (электромагнитная теория); (б) а = 161, ей- = 94%, 1Ш8 = 17% (скалярная теория); (в) а = 321, ейй = 96%, 1Ш8 = 49% (электромагнитная теория); (г) а = 321, ейй = 91%, 1Ш8 = 26% (скалярная теория).

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

4. Градиентный метод синтеза волновых полей

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

Пусть 10(х,у) - распределение интенсивности в плоскости z = 0, I (х, у)

-распре-

9 ■1

а) б)

0 ■

в) г)

Рис. 2. Распределение интенсивности от фокусатора в квадрат при различных параметрах

(21)

деление интенсивности в плоскости z = й. Требуется найти волновое поле Ж(х, у, г), удовлетворяющее следующим условиям:

Ж (х, у, / ) = ТЖ(х, у ,0)

10(х,у) =< Ж(х,у,0),Ж(х,у ,0) > ,

1 (х,у) =< Ж(х,7,/X Ж(х,у,/) >,

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

Ш1П е[Ж(х, у,/)], (22)

2

е = |(I(х,у)-д/< Ж(х,у,/),Ж(х,у,/) >) Жсйуу

Ж(х, у, /) = ТЖ(х, у ,0),

10(х,у) =< Ж(х,у,0),Ж(х,у ,0) > .

Для простоты будем искать минимум в классе функций, имеющих следующий вид:

Ж(x,у,0) = :2л/1сСх^у) ехР(ф(x,у))е(x,у,0),

(23)

где

Е (х, у,0) =

ео8 а (х, у)

- віп а (х, у) ехр(г'£ (х, у))

- віп а (х, у) ехр(г£ (х, у))

сов а(х, у)

8е=0.5

1(Х)

Ж(Х,/,

)Ж(х,/),ж(х,/)

Ж(Х,/),Ж(Х.

1(Х)

І\)Ж(Х,/),Ж(Х,/І --0.5^ 8Ж(Х,0),у(Х,0^ + )у(Х,0),Ж(Х,0))

У (х ,0) = Т-1

I (х)

'(Ж(х, / ),Ж(х, / )>

Ж(х, /)

5Ж ( х, /) = Т (5Ж ( х ,0)).

При выводе формулы (24) было использовано свойство сохранения псевдоскалярного произведения. Переходя к дискретному представлению и используя вышеприведенные выражения, получим выражение для градиента функционала в каждой точке (и, V) по аналогии с [1-2].

Уе(х,0) = (Цх,0),у(Х,0)) + [у(х,0),Цх,0)), (25) —

Отметим, что при таком выборе условие (21) выполняется автоматически. Этот вид представления поля выбран потому, что при данной форме записи каждый из параметров в приближении геометрической оптики имеет простой физический смысл. Функция ф (х, у) определяет фазу (или эйконал) волны в плоскости z = 0, а (х, у),5 (х, у) определяют поляризацию волны.

Найдем градиент функционала

где Цх,0) = — ехр(1ф(х,0))Е(х,0) .

Далее, используя выражение для градиента, можно построить итерационную процедуру минимизации функционала

фп+1 (х,0) = фп (х,0) - Nе, 1-шаг градиентного алгоритма.

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

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

5. Пример численного решения задачи синтеза

В качестве примера приведем расчет поля на выходе дифракционного оптического элемента (ДОЭ), фокусирующего в квадрат, расположенный на расстоянии / от плоскости ДОЭ (рис3). Шаг дискретизации

а) б) в)

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

на элементе А = 1, размеры апертуры элемента Б х = 64 1 , Бу = 641, / = 641 -фокусное расстояние, размер квадрата в плоскости фокусировки, а х Ь = 161х161 , а = 321 - параметр гауссова пучка.

На рисунке 3 а изображено распределение интенсивности от фокусатора гауссова пучка в квадрат, рассчитанное в приближении Френеля. Расчет поля также приведен в приближении Френеля. Дифракционная эффективность ейй=94%, среднеквадратичное отклонение гшв=16%.

На рис 3б изображено распределение интенсивности от фокусатора гауссова пучка в квадрат, но расчет поля проводился в рамках электромагнитной теории, ейй=96%, гшв=64%.

На рис 3в изображено распределение интенсивности от фокусатора гауссова пучка

в квадрат, рассчитанное в рамках электромагнитной теории, eff=97%, rms=22%.

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

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

1. Electromagnetic Theory of Gratings: Topics in current physics, v.22, Ed. by R.Petit, N.Y.: Springer-Verlag, 1980.

2. B.A. Зверев. Радиооптика. М. Советское Радио, 1975.

3. Воронцов M.A., Корябин. A.B., Шмальгау-зен В.И. Управляемые оптические системы. М.:Наука, 1988.- 272c.

4. Воронцов M.A., Шмальгаузен В.И. Принципы адаптивной оптики.- М.: Наука, 1985.-335c.

SYNTHESIS OF WAVE FIELDS IN FRAMEWORK OF ELECTROMAGNETIC

THEORY

© 1999 P.G. Serafimovich, S.I. Kharitonov Image Processing System Institute of Russian Academy of Sciences, Samara

In the last decade the issue of wave field determination on the base of intensity measurements in two planes is at the centre of attention. Such inverse problems arise in radio-astronomy; during calculation of diffraction optical elements, which form a specified intensity distribution; and for laser radiation focusing device. The latter is often used for dealing with technological problems, which are connected with laser treatment of materials. However, the problems of wave fields determination are solved in the framework of scalar theory mostly. In this paper the problem is solved in the framework of electromagnetic theory. The results of direct and inverse problems solving for scalar approximation and within the framework of electromagnetic theory were numerically compared in the paper.