Быстрый способ вычисления интеграла Релея-Зоммерфельда первого типа Текст научной статьи по специальности «Математика»

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

Андрей Владимирович Устинов (ведущий программист; e-mail: andr@smr.ru) Учреждение Российской академии наук Институт систем обработки изображений РАН

Аннотация

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

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

Введение

Интеграл типа Кирхгофа (в зависимости от граничных условий различают дифракционные интегралы Кирхгофа или Рэлея-Зоммерфельда) в скалярной или векторной форме позволяет достаточно точно описывать дифракцию световых полей в ближней зоне (<< длины волны) на малых фазовых объектах [1], на малых отверстиях в металлическом и черном экране [2].

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

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

2. Входное поле должно убывать до нуля при приближении к краям области вычислений.

3. Наиболее важно, что при увеличении расстояния начинает появляться эффект перекрытия спектров. Это в принципе можно преодолеть, применяя в пространстве частот более густую сетку и интерполируя значения спектра Фурье входного поля на новой сетке. Это, однако, довольно гро-моздко и будет приводить к неприемлемому увеличению времени вычисления и объёма исполь -зуемой памяти.

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

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

1. Общие соотношения

Интеграл Кирхгофа выражает решение дифракционных задач при помощи теоремы [4], согласно которой амплитуда электромагнитного поля U (х, у, z) на выходе выражается через её значение

U0( x, у,0) на входе:

ikz г г e I i |

U (x, y, z) = - — jjU0(u, v,0) — ll + — I dudv , (1)

2p

R

i

kR

где R = ^(м^у^+Т^ - расстояние между точками на входной и выходной плоскостях, к = 2р/1 - волновое число, 1 - длина волны падающего излучения, £ - область входной плоскости, в которой задано начальное поле.

В дальнейшем изложении учтём два обстоятельства: входная амплитуда задана не аналитическим выражением, а набором точек, получаемых при равномерной дискретизации по координатам с кусочно-постоянной аппроксимацией. Кроме того, временно не будем учитывать слагаемое / /(kR), о котором будет отдельно сказано позже. Так как в пределах каждой ячейки разбиения входная амплитуда является постоянной, фактические вычисления проводятся по формуле:

U (x ^z) = - Yu оК , vn]х

2p m,n

exp

ik*J (u - x)2 + (v - y)2 + z2

(2)

(u - x)2 + (v - y)2 + z2

dudv.

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

и0 [м , V ]

_01 т? п -I_

U (x, y, z) = - — Y- 2 2 2

2p mn (um - x)2 + (v„ - y)2 + z2

m+1 V„+i ._

j j exp i(u - x)2 + (v - y)2 +:

dudv

где (u m, vn) - центр ячейки.

S

u

v

2

u... v

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

2. Прямой метод В данном методе используется вычисление интегралов, входящих в формулу (3), непосредственно по квадратурным формулам (используется формула центральных прямоугольников). Отметим, что применено некоторое усовершенствование по сравнению с [5]. Оно состоит в том, что шаг интегрирования выбирается не на основе общей оценки для всей области, а выбирается для каждой ячейки, хотя, если мы находимся в условиях дальней зоны, выигрыш будет невелик. Примем в качестве шага интегрирования величину Т/8, Т - период экспоненциальной подынтегральной функции. Значение знаменателя 8 выбрано как наименьшее число дроблений, при котором период описывается достаточно точно. Поэтому при применении данного метода (и других тоже) необходимо вычислить или хотя бы оценить период функции, которая не является точной синусоидой.

Будем рассуждать следующим образом. Для точной периодической функции ехр(/'кх) частота равна

к и период 2 я / к . Получить частоту можно также следующим образом. Имеем равенство

ехр(гкх) ='к ехр(/'кх),

dx

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

d

— ехр

du

='к

(гку!(и - х)2 + (V - у)2 + г2) =

•у/(и - х)2 + (V - у)2 + 22 хехр (¡к^(и - х)2 + (V - у)2 + г2).

Из него следует, что величину

и - х

к ■ 1 1 у] (и - х)2 + (V - у)2 + 22

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

Т„;„ =

2я

(и - х)

2 . / \2 . 2 + (V - у) + 2

к\и - х

(4)

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

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

3. Переход к полярным координатам

Очевидный способ повышения быстродействия состоит в использовании того факта, что после замены переменных р = и - х; q = V - у подынтегральная функция в (3) становится радиально симметричной. В этом случае целесообразно сделать переход к полярным координатам, что позволит аналитически взять интеграл по одной из переменных, что существенно повысит скорость вычислений, так как квадратурная формула будет применяться лишь для одномерного интеграла. Опишем это более подробно.

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

11 /(Р2 + q 2)dpdq,

Р1 = ит - х Р2 = ит

ql = vn- у q2 = Vn+1

-у

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

Ц /(гф, в котором область Б - некоторый

в

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

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

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

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

4) прямоугольник соприкасается с обеими осями. Варианты 1 и 4 также разделяются на два подва-

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

- х

и - х

тах

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

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

С В

f arccos( pi / r) Л

42

ч i

D

гс

Га

Гв

Ч>

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

Преобразование границ описывается соотношениями:

AB : p = p2 ® A' B': r cos j = p2 BC : q = q2 ® B' C': r sin j = q2 CD : p = p1 ® C'D': r cos j = p1 DA : q = q1 ® D' A': r sin j = q1 Рисунок является схематичным: выпуклость или вогнутость границ не анализировалась.

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

j j f (p2 + q2)dpdq = jj f (r )rdrd j =

pi q W

= jj (.) + jj (.)+jj (.X

Wi

rB f arcsin(q2 / r) Л

jj f (r)rdrd j= j f (r )r j d j

W1 rC y arccos(p2/ r) у

rB

= j f (r )r (arcsin — - arccos—)dr r

rC

rC f arccos

dr =

r

(pi/r) Л

jj f (r)rdrd j= j f (r )r j d j

W2 rA У arccos(p2 / r)

rC

= j f (r)r (arccos — - arccos—)dr,

j 1* 1*

dr =

jj f (r)rdrd j= j f(r )r j d j

W3 rD y arcsin(qi /r)

rA

= j f (r)r(arccos— - arcsin-qL)dr.

j 1* 1*

dr =

Если теперь использовать формулу агссо8 х = р /2 - агсБт х , раскрыть скобки и объединить интегралы с одинаковой подынтегральной функцией, то получим, что

jj I (г)гdгd ф =

= j I (г )г а1тат —2 dг + j I (г )г а1тат —2 dг -

гС Г гА Г

гС гА

- j I (г )г а1тат Р dг - j I (г )г атсБШ —1 dг -

гв г гв г

р гВ р гА

- - j I(г )^г + - j /(г)гdг.

гС гв

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

'2

j f (r)rdr и

r

r2 a

j f (r)r arcsin —dr .

(5)

(6)

В формулах (5) и (6) г1, г2 - полярные радиусы некоторых вершин прямоугольника, а - одно из чисел р1, р2, —х, q2, причём гарантировано, что аргумент арксинуса не превосходит единицы.

Теперь в написанные выше выражения подста-

¡кл]г2 +

вим функцию из (3) f (r) = exp

и полу-

чим, что результат является линейной комбинацией интегралов

j exp z'WÍ

Ir2 + z2

j exp z'kv/j

lr2 + z2

rdr и

r arcsin— dr.

r

(7)

(8)

Интеграл в формуле (7) берётся в явном виде:

r2 |- ,-

j exp zk4r2 +

rdr =

■■-i(exp z'Wr2 + z2 -(i-ik*Jr2 + z2))

(9)

Интеграл (8) вычисляется по квадратурной формуле, при этом шаг выбирается равным 1/8 периода, так же как и в параграфе 1. Значение периода находится аналогично выводу формулы (4) и равно:

A

W

B

B

W

W

A

Т„;„ =

2р

г2 + т 2

тах

Ю„

кг

(10)

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

. а

том, что в этой ситуации множитель г агсвт — из-

г

меняется на отрезке интегрирования быстрее экспоненты. Поэтому сделаем следующее преобразование:

| ехр 1к-у/

г агсБШ— ёг ■■ г

■■ ехр

] • |

(11)

г агсБШ— ёг. г

где г = (г + г2)/2 .

Интеграл в правой части (11) берётся в явном виде:

| г агсвт — ёг = [а / г = вт х] = -а2|

- ёх =

2 вт х 2

- + — = а

вт х агсвт(— / г) 2вт2(агсвт(а / г))

(12)

а г а а 1~2 2

I--ctg(arcsin(— / г)) = — агсвт—I—л/г - а .

2 2 г 2

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

| ехр 'к

¡г2 + т2

г агсБШ— ёг ■■ г

ехр

('^г* + т2) I Г— агсзша + 2

. (13)

Теперь выясним, когда возможно применение преобразования (11). Можно доказать, используя выражение для погрешности квадратурных формул [6] и считая, что шаг интегрирования равен Т/8,

что, если отбросить общий множитель -, отно-

16 • 24

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

■ а 0,679 , Л А

функции г агсвт— - величине -- (оценка), А -

г р гАю2

размер прямоугольника, ю= 2я / Т. Таким образом, если выполняется неравенство

1,258

2 <-

(14)

яг1Аю2 '

то преобразование (11) возможно. Отметим, что правая часть неравенства (14) оценена довольно грубо. Если г| велико, то в числителе будет число большее, чем 1,258. Но в этом случае знаменатель тоже увеличится - см. формулу (10).

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

4. Использование асимптотической формулы

в полярных координатах При использовании алгоритма, описанного в пункте 3, наибольшие затраты времени происходят при вычислении на ячейках, далёких от начала координат. Это очевидно, так как шаг интегрирования на них меньше; в то же время их вклад в общую сумму сравнительно мал. Но из этого противоречия есть удачный выход - применение асимптотической формулы. Строго говоря, это название является условным, так как это не есть асимптотическое разложение в привычном его понимании. Мы учитываем то, что на удалённых ячейках арксинус мало изменяется:

| ехр

г атсБШ— ёг г

■ атсБш— I ехр

V 3

гёг =

1 . а —агсБШ— х

к г.

(15)

ехр ('к^/г22 +; ехр ('к^2 +; г, = (г! + г2) /2.

•(1 -'кТТ2

•(1 - 'кЛ1

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

г2

Для интеграла |е'югёг относительная погреш-

г

ность вычисления по квадратурной формуле с ша-, Т р

гом п = — = — независимо от значений ю, г,, г2

8 4ю 12

Р2 я

равна . Для интеграла от арксинуса при замене

х

а

=а

, г

2

а

2

подынтегральной

функции

константой

Гагсвта ёг » (г2 - г1) агсвша относительная по-

грешность, вычисленная по формуле (17), примерно

равна

а(г2 - г)2 2г/ - а2

. а г (г - а2)3 атсБш— ^1 7

. Если она более чем в

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

а(г2 - г)2 2г/ - а2

. а г (г - а2)3 агсвш— , 4, 7

384

(16)

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

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

5. Запоминание интегралов по ячейкам в матрице с её последующим использованием

Ещё раз запишем сумму из формулы (3):

и0[Мт , V, ]

(и т - х)2 + (V, - у)2 + 2

| | ехр 'ку[

222 I р + q + 2

dpdq,

р = и - х, q = V - у.

Возможность повышения быстродействия обусловлена тем, что частичные интегралы в этой сумме повторяются многократно с одинаковыми значениями рт, рт+1, qn, qn+1 для различных точек (х, у) выходной плоскости. Это означает, что если нужно вычислить одно значение и(х, у,т), то факт повторения использовать невозможно. Но, если требуется найти целое поле и(х, у,т) и точки (х, у) расположены на дискретной сетке, то можно использовать тот факт, что величины рт, qn принимают конечное число значений. Это позволяет заранее вычислить для них интегралы, запомнить их в массиве и затем использовать их при вычислении вышеуказанной суммы.

Идея такого подхода взята из статьи [5], но в ней авторы не учли (во всяком случае, в статье об этом явно не сказано), что многократные повторения будут лишь в том случае, когда шаги дискретизации входного и выходного поля одинаковы. Если это ус-

ловие не выполнено, то повторений будет меньше и соответственно размер матрицы будет возрастать, что приведёт к уменьшению ускорения. Поэтому далее будем считать, что данное условие выполнено, тем более, что оно не является обременительным. (В принципе, можно использовать матрицу и в случае, когда шаги кратны друг другу, но при этом ускорение уменьшается и возрастает сложность программы, так как для разных соотношений шагов алгоритм несколько отличается.) Если входное поле имеет N х N точек, а выходное - Ы2 х Ы2 точек, то можно доказать, что размер матрицы будет N + ы2 -1) х (N + ы2 -1).

Данный подход был программно реализован, было получено существенное ускорение, особенно после программной оптимизации цикла. До использования матрицы оптимизация не производилась, так как основное время тратилось на вычисление частичных интегралов. Было обнаружено, что почти всё время (даже после оптимизации) затрачивается на непосредственное вычисление выходного поля, а не на формирование матрицы. Однако это не значит, что преобразования, описанные в пунктах 3 и 4, оказываются бесполезными. Благодаря им получено ускорение большее, чем в [5], где частичные интегралы вычислялись по двумерным квадратурным формулам и получили упомянутое соотношение времён формирования матрицы и вычисления поля. С учётом этого соотношения, в отличие от [5], мы не стали использовать возможность вычислять не всю матрицу, а только её восьмую часть.

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

6. Результаты вычислений

Здесь приводятся результаты вычислений, выполненных программой, реализующей описанные выше алгоритмы. Все эксперименты проводились при следующих значениях параметров: длина волны излучения 633 нм, размеры входной и выходной областей 4х4 мм с разбиением 51 х51 точек. Конечно, такой шаг дискретизации в ряде случаев велик, а выходная область обычно берётся больше входной, но приводимая таблица предназначена только для сравнения времени выполнения различных алгоритмов. Входная функция - круглое отверстие. Расстояние между плоскостями является переменным. В таблице 1 приводится время выполнения программ и величина ускорения. Вычисления производились на компьютере, имеющем 512 МБайт оперативной памяти и тактовую частоту 2,01 ГГц.

Мы видим, что по сравнению с прямым методом алгоритм с использованием матрицы даёт ус-

г

г

V

г

Л

2

<

г

Л

т ,п

q

корение как минимум в 34 раза, и более чем в 70000 раз в ближней зоне. В работе [5] получено для сопоставимых условий ускорение в 4 и более 1440 раз. У нас получено существенно большее ускорение благодаря использованию аналитических преобразований и, возможно, лучшей программной реализации. Также был повторён эксперимент, описанный в работе [3]. В этой статье сравнивался метод, предлагаемый авторами, и быстрые алгоритмы метода распространения плоских волн и аппроксимации Френеля. У нас было затрачено время 3,6 секунды, что существенно меньше полученного в [3] (5 минут для предлагаемого метода и 10 секунд для быстрых алгоритмов). Отметим, что эксперименты, описанные в [3], производились на более мощном компьютере: 4 ГБайт оперативной памяти и тактовая частота 2,8 ГГц.

Таблица 1. Время вычисления для области 4х4 мм с разбиением 51х51 точек

7, мм Прямой метод (п.2) Полярный 1 (п.3) Полярный 2 (п.4) Матрица (п.5)

2000 12 с 40с 36с 0,34с

1000 46 с 44с 36с 0,34с

500 3м03с 58с 35с 0,34с

250 12м10с 1м29с 35с 0,34с

125 =50м 2м33с 36с 0,34с

70 =400м 4м15с 37с 0,34с

30 =1200м 9м23с 41с 0,34с

10 — 30м13с 55с 0,35с

1 — 156м 3м38с 0,43с

Точность работы программы оценивалась визуально следующим образом. Было проведено сравнение результатов, полученных матричным методом и прямым методом, реализованным в программе SimuLight.exe. Во всех случаях падающая волна была плоской. Примеры работы программы показаны на рисунках 2-4.

а) _б) I

Рис.2. Входное поле: Я=1 мм, 100 отсчетов, а) интенсивность, б)фаза

ш

Рис 3. Выходное поле: расстояние г=300 мм, К=2 мм, 200 отсчетов, а) интенсивность (негатив), б)фаза. Использован матричный метод

Рис.4. Выходное поле: расстояние z=1000 мм, К=2 мм, 200 отсчетов, а) интенсивность (негатив), б)фаза.

Использован матричный метод

Обе программы дают результаты, примерно одинаковые визуально. Среднеквадратическое отклонение интенсивности равно 6,6% на расстоянии 300 мм и 4% на расстоянии 1000мм. (Визуально такой же результат даёт программа Fresnel.exe, использующая приближение Френеля.) Другой способ сравнения следующий. Если расстояние между плоскостями очень мало, то выходное поле примерно равно входному. Такой пример показан на рисунках 5, 6. Наличие чёрной полосы на рисунке 6 а предсказывается теоретически.

Рис.5. Входное поле: скачок фазы, 200 отсчётов, а)интенсивность, б)фаза

Рис.6. Выходное поле: расстояние z=0,01 мм, 200 отсчётов, а)интенсивность, б)фаза

7. О возможности считать постоянной величину 1/К2 + 1/Ш3 Выше, в формуле (3) мы предположили, что ве-

личина

1/((и - х)2 + (V - у)2 + z ) + I / к ((и - х)2 + (V - у)2 + z2)3/2

в пределах

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

В самом начале было сказано, что слагаемое I /(кК) в подынтегральном выражении формулы (1) не учитывается. Условие, когда это можно делать,

достаточно очевидно - кЯ >> 1 [4]. Но, как мы увидим, это слагаемое целесообразно сохранять во всех случаях. Так сделано при программной реализации в методе с использованием матрицы (в остальных методах слагаемое ' /(кЯ) не учитывается). Условия вынесения за знак интеграла рассмотрим вначале для каждого слагаемого в отдельности.

Пусть дана функция одной переменной /(г). Если мы запишем следующее приближённое равен-

ство:

'2

|/(г)ёг » /(г,)(г2 -г1), г, = (г + г2)/2 г

относительная погрешность примерно равна

/V, )(г2 - г1)2

5 = -

24/ (г.)

то его

(17)

В нашем случае /(г) = —-2 , поэтому, исходя

г + г

из (17), имеем:

5 =

(3г/ - 22)(г2 - г1)2 12(г2 + 22)2

(18)

Разность г2 - г1 < А при любом расположении ячейки (равенство может достигаться, если сторона ячейки лежит на координатной оси). Функция

1 2 2 3г. - 2

рассматриваемая как функция одной пе-

(г.2 + 2 2)2

ременной г., принимает наибольшее (по модулю)

I I А2

значение при г, = 0. Таким образом, 5 <-т-.

12т2

Сравнивая это значение с величиной , получа-

ем, что вынесение функции

—-2 за интеграл не

г + г

приводит к ухудшению точности, если выполняется неравенство:

г >^31 А» 1,801А . (19)

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

/(г) = "Л—12 3/2 . Исходя из (17), имеем:

5 =

(г 2 + 2 2)3

(4г,2 - г2)(г2 - г)2 8(г2 + т2)2

(20)

Анализируя правую часть (20), получим, что |5 , откуда получаем, что вынесение функции

1

(г2 + 2 2)3/2

точности, если выполняется неравенство:

2 А»2,205А .

за интеграл не приводит к ухудшению

Теперь рассмотрим оба слагаемых вместе: вычислим погрешность для функции

1

г2 + 2 2 к (г2 + 2 2)3'2 Аналогично (18) и (20), имеем:

(22)

5 =

3г 2 - 2 2

3' 4г2 - г2

'(г2 + т2)3/2 к (г? + т2)2 24 {(г,2 + т2)1/2 +' / к}

(гг - П)

(23)

Разность г2 - г1 < А . Найдём максимум модуля оставшейся функции, обозначим её g(г,). Не будем выписывать его явный вид и приравнивать к нулю производную. При рассмотрении слагаемых по отдельности было доказано, что значение g(г,) при г, = 0 больше, чем при нуле производной. Поэтому предположим, что максимум модуля функции g (г,) будет при г, = 0 .

1

5

тах| g (г, ^ =|g (0)| = + т 2 к 2 + Г

Таким образом,

5 <

А2

24т2

4+

т2к2 +1

(24)

Отсюда видно, что погрешность лежит между значениями, соответствующими каждому слагаемо-

А2

му. При большом к получим |5| < ^ 2 , а при малом

12г2

А2

кг будет 5 <—- . В обоих случаях получим форму-8т

лу для превалирующего слагаемого. Если нет явного преобладания одного из слагаемых, точное граничное значение для т получается из (24) решением уравнения 6 степени, которое сводится к кубическому. Можно также использовать оценку сверху -формулу (21). (Фактически в программе граница была увеличена до 10А .)

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

( 1 ^ ^-г +

| ехр 'к

¡г2 + т2

| ехр 'к

¡г2 + т2

к (г2 + т2)3/2 1

гёг и

(25)

+

V к (г2 + т2)3/2 ,

гагсБт— ёг . (26)

г

Вычислим интеграл (25).

5

1

г 2 + 2 2

j exp ikj,

r

t = WT

2 , 2 r + z

1I-2, 2 \3/2

k (r + z )

rdr =

/

r2 + z2

4

= j—dt + — k j—dt = t t kit2

t2 it

=j e~dt+ij 7dt=j 7d

(27)

t2 ea

t2- eit

dt +

+i

e"i e"2 % e"

---+ il — dt

t t t '1 '2 t 1

Л f =i

/

v

ti

2 У

где t1 = k^2 + z2, t2 = kj

r22 + z2

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

Для интеграла (26), если выполняется условие (16), можно применить приём, аналогичный формуле (15): ( 1 ^

j exp ik-v/i

lr2 + z2

k (r2 + z2)3

(28)

f it1

xrdr » arcsin— i r

e

V ti

jt2 Л

e

t

a

arcsin—

2

Заключение

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

Литература

1. Totzeck, M. Validity of the scalar Kirchhoff and Rayleigh-Sommerfeld diffraction theories in the near field of small phase obiects // J. Opt. Soc. Am. - 1991.- Vol.8, N1. - P.27-32

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

3. Veerman, J.AC.- Calculation of the Rayleigh-Sommerfeld diffraction integral by exact integration of the fast oscillating factor / J.A.C.Veerman, J.J.Rusch, HP.Urbach // J. Opt Soc. Am. - 2005.- Vol.22, N4. - P.636-646

4. Виноградова, М.Б. Теория волн / Виноградова М.Б., Руденко О.В., Сухоруков А.П.- М.: Наука, 1979. - 383 с.

5. Балалаев, С.А. Реализация быстрого алгоритма преобразования Кирхгофа на примере Бесселевых пучков / С.А.Балалаев, С.Н.Хонина. // Компьютерная оптика

- 2006. - Т.30 - С.69-73. - ISSN 0134-2452.

6. Ланцош, К. Практические методы прикладного анализа (глава VI) / К. Ланцош; пер. с англ. - М.: Государственное издательство физико-математической литературы, 1961. - 524 с. (Cornelius Lanczos. Applied Analysis, Prentice Hall, Inc., 1956).

References

1. Totzeck, M. Validity of the scalar Kirchhoff and Rayleigh-Sommerfeld diffraction theories in the near field of small phase obiects // J. Opt. Soc. Am. - 1991.- Vol.8, N1. - P.27-32

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

3. Veerman, J.A.C.- Calculation of the Rayleigh-Sommerfeld diffraction integral by exact integration of the fast oscillating factor / J.A.C.Veerman, J.J.Rusch, H.P.Urbach // J. Opt. Soc. Am. - 2005.- Vol.22, N4. -P.636-646

4. Vinogradova, M.B. Waves theory / Vinogradova M.B., Rudenko O.V., Suhorukov A.P.- Moscow: Science, 1979.

- 383 p. - (in Russian).

5. Balalayev, S.A Realization of fast algorithm of Kirchhoffs diffraction integral on an example of Bessel modes / S.A.Balalayev, S.N.Khonina. // Computer Optics. - 2006. -Vol.30 - P.69-73. - ISSN 0134-2452. - (in Russian)/

6. Lanczos, C. Applied Analysis / Prentice Hall, Inc., 1956.

- 524 p.

(

\

1

Л

e

+

r2 + z2

X

r

s

the fast way for calculation of first class rayleigh-sommerfeld integral

Andrey Vladimirovich Ustinov (leadingprogrammer; e-mail: andr@smr.ru) RAS Image Processing Systems Institute

Abstract

Ways for acceleration of Kirchhoff integral owing to using of analytical transformations and the effective program realization are described in present article. Consecutive steps of the speed increase after every transformation are shown. Application conditions of different approximations are obtained. It should be noted that after an application of all transformations the time of program performance weakly depends on the distance between input and output planes.

Key words: Kirchhoff integral, paraxiality, quadrature formula, asymptotic expansion.

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