CC BY
29
ВАРИАЦИОННЫЙ ПОДХОД К РАСЧЁТУ ФУНКЦИИ ЭЙКОНАЛА
12 2 12 12 Л.Л. Досколович ' , А.А. Мингазов , Д.А. Быков ' , Е.С. Андреев '
1 Самарский национальный исследовательский университет имени академика С.П. Королева, 443086, Россия, Самарская область, г. Самара, Московское шоссе, д. 34, 2 ИСОИ РАН - филиал ФНИЦ «Кристаллография и фотоника» РАН, 443001, Россия, Самарская область, г. Самара, ул. Молодогвардейская, д. 151
Аннотация
Задача расчёта функции эйконала из условия фокусировки в заданную область сформулирована как вариационная задача и как задача Монжа-Канторовича о перемещении масс. Получено, что функция стоимости в задаче Монжа-Канторовича соответствует расстоянию между точкой исходной области, в которой задана функция эйконала, и точкой области фокусировки. Предложенный в работе формализм позволяет свести расчёт функции эйконала к решению задачи линейного программирования. При этом расчёт «лучевого отображения», соответствующего функции эйконала, сводится к решению линейной задачи о назначениях. Предложенные вариационные подходы проиллюстрированы на примере расчёта оптических элементов для фокусировки пучка круглого сечения в прямоугольную область.
Ключевые слова: геометрическая оптика, функция эйконала, вариационная задача, задача Монжа-Канторовича.
Цитирование: Досколович, Л.Л. Вариационный подход к расчёту функции эйконала / Л.Л. Досколович, А.А. Мингазов, Д.А. Быков, Е.С. Андреев // Компьютерная оптика. -2018. - Т. 42, № 4. - С. 557-567. - Б01: 10.18287/2412-6179-2018-42-4-557-567.
Введение
Задача расчёта оптического элемента из условия формирования заданного распределения освещённости в некоторой плоскости относится к классу обратных задач неизображающей оптики и является крайне сложной. В большинстве случаев данная задача сводится к решению нелинейного дифференциального уравнения эллиптического типа. Методы расчёта оптических поверхностей, основанные на прямом численном решении уравнения данного типа, появились только в последние годы [1 - 4]. Общим недостатком указанных методов является высокая вычислительная сложность. Поэтому для решения обратных задач не-изображающий оптики широко используются различные итерационные методы. Одним из универсальных и широко используемых итерационных методов является метод согласованных квадрик (МСК, supporting quadric method) [5 - 12]. Наиболее известен вариант метода для расчёта зеркал или преломляющих оптических элементов, формирующих дискретные распределения интенсивности (освещённости) в виде набора точек. В зависимости от задачи в качестве квадрик используются параболоиды, эллипсоиды или гиперболоиды. В частности, в задаче расчёта зеркал поверхность зеркала представляется в виде набора сегментов параболоидов (задача формирования распределения интенсивности в дальней зоне) или эллипсоидов (задача фокусировки в набор точек в ближней зоне) с определёнными параметрами. При этом число сегментов равно числу точек фокусировки. Расчёт параметров параболоидов (эллипсоидов) осуществляется итерационным методом, при этом сходимость метода строго доказана [6]. В работе [12] предложена модификация МСК для расчёта функции эйконала светового поля из условия фокусировки в набор точек. Данная задача является важной при рас-
чёте дифракционных оптических элементов, в особенности преобразователей формы пучка [13]. Кроме того, по функции эйконала может быть восстановлена рефракционная оптическая поверхность, что позволяет применять предложенный метод для расчёта преломляющих оптических элементов [14].
Основными достоинствами МСК являются его универсальность и возможность формирования сложных дискретных изображений. В то же время при решении светотехнических задач обычно требуется формирование непрерывных распределений в областях относительно простой формы (прямоугольник, эллипс, треугольник и т.п.). В базовых теоретических работах [11, 15] показано, что проблема расчёта зеркала для формирования заданного непрерывного распределения интенсивности в дальней зоне может быть сформулирована как вариационная задача минимизации функционала и как задача Монжа-Канто-ровича о перемещении масс со специальной функцией стоимости. Важным практическим результатом указанных работ является сведение вариационной задачи к задаче линейного программирования (ЗЛП). В силу специфического вида матрицы ограничений, для решения соответствующей ЗЛП существуют эффективные алгоритмы, обладающие полиномиальной сложностью [16]. Расчёт зеркал на основе решения ЗЛП рассмотрен в [17, 18].
В настоящей работе вариационный подход из работ [11, 15] впервые применён к решению задачи расчёта функции эйконала светового поля. Показано, что задача расчёта функции эйконала из условия фокусировки в заданную область также может быть сформулирована как вариационная задача и как задача Мон-жа-Канторовича о перемещении масс. Отметим, что в отличие от работ [11, 15], в задаче расчёта эйконала не требуется выполнение условия дальней зоны. Под
массами в задаче Монжа-Канторовича понимаются исходное распределение освещённости и требуемое распределение освещённости в области фокусировки. При этом получено, что функция стоимости соответствует расстоянию между точкой исходной плоскости (в которой задана функция эйконала) и точкой области фокусировки. Данный результат показывает, что решение обратной задачи расчёта эйконала соответствует отображению, для которого минимизируется суммарное расстояние между точками исходной плоскости и области фокусировки. Предложенный формализм позволяет свести расчёт функции эйконала к ЗЛП. При этом расчёт отображения, соответствующего функции эйконала, может быть сведён к решению линейной задачи о назначениях [19]. В качестве примеров рассчитаны функции эйконала из условия фокусировки пучка круглого сечения в прямоугольную область. По функциям эйконала рассчитаны преломляющие оптические элементы.
Текст статьи организован следующим образом. В параграфе 1 приведена постановка задачи. В параграфе 2 определяется так называемое слабое (обобщённое) решение задачи. В параграфе 3 сформулирована соответствующая вариационная задача. В параграфе 4 рассмотрена связь вариационной задачи с задачей Монжа-Канторовича. В параграфе 5 рассмотрен расчёт функции эйконала из условия фокусировки пучка круглого сечения в прямоугольную область. Расчёт проведен на основе решения ЗЛП и задачи о назначениях. Для удобства читателя доказательства лемм и теорем вынесены в приложения.
1. Постановка задачи
Пусть в области О, расположенной в плоскости г = 0, задано световое поле с распределением освещённости Е0(и!, и2) и функцией эйконала Ф0 («ь и2), где («ь и2) - декартовы координаты. Область О будем называть апертурой. Обозначим
Г («1, И2) = ^Ф „г Ф и2^1 -(УФ)2 ^
единичный вектор луча для среды с показателем преломления п = 1, где Ф«. = дФ/ ди, I = 1, 2, УФ = (Фи1, Ф%). Рассмотрим «лучевое отображение» уФ области О в область на плоскости г =/> 0, задаваемое функцией эйконала Ф(иь «2). В данном отображении точка (хь х2) = уФ(«ь «2) определяется как точка пересечения луча с плоскостью г =/ (рис. 1):
Х1 = «1 + Ф
Х2 = «2 +Фи
, //^ -(УФ)2
2 //л/1 -(УФ)2.
(1)
Для удобства записи далее будем использовать следующие обозначения. Точки («ь «2) и (хь х2) будем обозначать и и х соответственно. Интеграл по ^ будет обозначать двойной интеграл по du\du2
Распределение освещённости в области фокусировки определяется из закона сохранения энергии вдоль лучевых трубок в виде:
Е (х) = Е (уф (и)) = Ео (и ) / 3 (уф (и)) :
(2)
где ,/(уФ(и)) = дх! / д«• дх2/ ди2- дх2/ дщ-дх! / ди2 - якобиан преобразования координат (1).
Щ и)/
и2
щ
Уф
Е(х)
Х2
XI
У
ъ
2=0 ' г=/
Рис. 1. Геометрия задачи расчёта функции эйконала Соотношение (2) может быть представлено в интегральном виде:
| Е0(и) du = |Е(х)dx,
(3)
где ю с О - любая область (любое измеримое подмножество). Закон сохранения энергии может быть также записан в следующем виде (аналогично см. лемму 4.9 в [11])
]И(уФ (и)) Е0 (и) d« = |И( х) Е (х) dx, (4)
О Б
где И (х) - любая непрерывная функция.
Сформулируем обратную задачу расчёта эйконала светового поля из условия формирования заданного распределения освещённости. Требуется рассчитать функцию эйконала Ф (и), иеО, индуцирующую биективное лучевое отображение уФ: О ^ Б, которое формирует в области Б плоскости г =/ заданное распределение освещённости Е(х), хеБ (рис. 1).
2. Слабое решение В работах [11, 15], посвященных расчёту зеркала для формирования заданного распределения интенсивности в дальней зоне, вводится понятие так называемого слабого (обобщённого) решения. В данном параграфе мы введём аналогичное понятие для задачи расчёта функции эйконала.
Подставляя (1) в (2), можно видеть, что расчёт функции Ф(и) сводится к решению нелинейного дифференциального уравнения в частных производных эллиптического типа. Требование дифференци-руемости функции эйконала Ф(и) и отображения уФ(«) существенно ограничивает класс решаемых задач (класс формируемых распределений освещённости). В связи с этим введём понятие слабого (обобщённого) решения, позволяющего расширить класс допустимых функций Ф(«). Для определения слабого решения рассмотрим метод построения функции эйконала Ф(и) из эйконалов линз с фокусами в точках хеБ [12]. Под функцией эйконала линзы Ф? (и; х) понимается такая функция эйконала в плоскости г = 0, при которой все лучи приходят в точку х в плоскости г = / Из (1) несложно получить, что функция эйконала линзы имеет вид:
Ф/ (и; х) = ¥( х )-р(и, х), (5)
ю
где р(и, х) = у /2 + (х - и)2 - расстояние между точками с координатами (иь и2, 0) и (х1, х2, /). Действительно, подставляя (5) в (1), получим, что х = уФ(и). Кроме того, заметим, что уравнение (5) при г = 0 соответствует эйконалу сходящегося сферического
пучка Тъ (и, г) = -у](/ - г)2 + (х - и)2 + Т(х) с фокусом в точке хеБ, записанному при г = 0. При этом константа Т(х) соответствует значению эйконала в точке фокусировки. Функция эйконала для фокусировки в заданную область Б может быть представлена в виде огибающей семейства функций (5) по параметрам х = (х1, х2)еБ [12]. Действительно, по определению, огибающая поверхность касается каждой функции эйконала в семействе (5) в некоторой точке. При этом частные производные огибающей поверхности совпадают в точке касания с частными производными эйконала линзы. Направление луча определяется частными производными функции эйконала. Поэтому лучи, соответствующие эйконалу в виде огибающей поверхности, будут направлены на точки области Б, которые являются фокусами линз в (5). Уравнение огибающей поверхности для семейства эйконалов линз (5) имеет вид [12]:
ФI (и; х) = х )-р(и, х),
д
< — ((х )-р(и, х) ) = 0, (6)
д
—((х )-р(и, х) ) = 0.
Система уравнений (6) содержит эйконал линзы (первое уравнение) и производные данного уравнения по параметрам (х1, х2) (второе и третье уравнения). Отметим, что второе и третье уравнения в системе (6) являются необходимым условием экстремума функции (5) по переменным (х1, х2) при фиксированном значении и = (и1, и2). Данное условие соответствует принципу Ферма и позволяет представить функцию эйконала в переменных (и1, и2) в виде:
Ф(и) = тах (Т( х )-р(и, х)), (7)
хеБ У У ' '
или
Ф(и) = тт (т( х)-р(и, х)). (8)
хеБ У У ' '
Заметим, что максимум в формуле (7) для Ф(и0) достигается в точке х0 такой, что х0 = у(и0). Действительно, пусть максимум достигается в точке х0, тогда функция эйконала линзы Т(х0) - р(и, х0) касается функции Ф(х0) в точке и0. Поскольку обе функции дифференцируемы, то производные Ф в этой точке совпадают с производными функции эйконала линзы, а отображение у задаётся производными функциями, то уФ(и0) совпадает с образом точки и0 под действием отображения у, индуцированного функцией эйконала линзы с фокусом в х0. Аналогично для функции Ф, определяемой формулой (8).
Представления (5) зависят от функции ¥(х), которая определяет распределение энергии в области фокусировки. По построению функция ¥(х) соответствует эйконалу в области D и также может быть представлена в виде огибающей семейства эйконалов линз, заданных в области D и фокусирующих в точки и eG. Функция эйконала линзы, обеспечивающая фокусировку из области D в точку ие G, имеет вид:
¥I (х;и) = Ф(и) + р(и, х), (9)
где константа Ф(и) соответствует значению эйконала в точке фокусировки. Отметим, что в отличие от (5), величина р(и, х) входит в (9) с противоположным знаком. Это связано с тем, что в данном случае лучевое отображение определяется «инвертированным» единичным вектором луча
r( х1, хг) = - (Y х!, ¥ х2^1 -(VT)2 j,
который имеет отрицательную компоненту z. Кроме того, заметим, что уравнение (9) соответствует эйконалу сферического пучка
T sb (х, z ) = yj z 2 +(х - и )2 + ¥ (х),
расходящегося из точки (и, 0), записанному при z = f Огибающая семейства (9) по параметрам и = (и1, и2) имеет вид:
TI (х;и) = Ф(и) + р (и,х),
д
—(ф(и ) + р (и, х)) = 0, (10)
д
— (ф(и ) + р(и, х) ) = 0.
ди2 4 '
Второе и третье уравнения в системе (10) являются необходимым условием экстремума функции (9) по переменным (иь и2) при фиксированном значении х = (хь х2). Это позволяет представить функцию в виде:
¥( х) = max (ф( и ) + р (и, х)), (11)
neG
или
¥( х) = min (ф(и ) + р (и, х)). (12)
neG
Покажем, что если прямое отображение уФ: G ^ D определяется функцией (7), то обратное отображение у-1: D ^ G определяется функцией (12). Действительно, пусть х0 = уФ(и0), где
х0 = arg max (t(х)-р(и0, х)).
xeD
Покажем, что в этом случае и0 = у¥(х0). Докажем это от противного. Предположим, что и! = у¥(х0), где и! Ф и0. В этом случае, согласно (12), будем иметь
T(х0) = Ф(и1 ) + р(и1,х0) < Ф(и0) + р(и0,х0)
и Ф(и0) > ¥(х0) - р (и0, х0). Последнее неравенство противоречит исходному условию х0 = уФ (и0), согласно
которому Ф(и0) > ¥(х0) - р (u0, х0). Таким образом, ут = уф1. Аналогично можно показать, что если прямое отображение уФ: G ^ D определяется функцией (8), то обратное отображение уф1: D ^ G определяется функцией (11). Далее будем предполагать, что функции эйконала на апертуре G и в области фокусировки D связаны соотношениями:
Ф(и) = max х )-р (u, х)),
t , ч X (13)
Y (х) = min (Ф(и ) + р (u, х)).
ueG
При этом оба экстремума в (13) достигаются, когда u и х связаны соотношением х = уФ(и) (u = Уф1( х) = ут (х)). Из уравнений в (13) легко получить, что функция эйконала Ф (u) является «огибающей сверху», поскольку для любой функции эйконала линзы Ф? (u; х), входящей в семейство (6), выполняется условие Ф (u) > Ф? (u; х).
Отметим, что функции Ф (u) и Т(х) связаны неравенством Т(х) - Ф (х) < р (u, х) которое превращается в равенство в случае, когда u и х связаны соотношением уФ (u) = х. Это позволяет переопределить отображение уФ следующим образом:
УФ (u) = {х e D | х) -O(u) = р (u, х)}. (14)
Обобщим понятие решения исходной задачи формирования заданного распределения освещённости Е(х), xeD, не предполагая дифференцируемости функции Ф (u) и Т(х).
Определение 1. «Хорошей парой» будем называть пару функций (а (u), ß (х)), таких что a(u) - непрерывна на G, ß (х) — непрерывна на D, и функции а (u) и ß (х) связаны соотношениями
a(u) = max {ß (х) - р (u, х)},
xeD
ß (х) = min {a(u) + р (u, х)}.
ueG
(15)
Отметим, что на функции а(и) и р (х) не накладывается требование дифференцируемости. При этом Р (х) - а (и)< р (и,х) УиеО, УхеБ. В приложении 1 показано, что функции а (и), р (х) являются липшицевы-ми в областях О и Б соответственно. Согласно теореме Радемахера липшицева функция в области является дифференцируемой почти всюду в данной области (то есть мера множества точек, в которых функция не дифференцируема, равна нулю). Значит, для хорошей пары функции а (и) и р (х) почти всюду дифференцируемы и соответствующее отображение уФ (см. (14)) однозначно в точках дифференцируемости.
Теперь сформулируем, что следует понимать под слабым (обобщённым) решением.
Определение 2. Слабым решением называется «хорошая пара» функций (Ф (и), Т(х)), удовлетворяющая условию сохранения энергии (3) для отображения, определённого формулой (14).
Отметим, что в слабом решении производные функции Ф (и) и, соответственно, отображение (1) определены почти всюду. При этом множество точек недифференцируемости не влияет на величину интегралов в (3).
3. Вариационная формулировка задачи
Задача поиска слабого решения, рассмотренного в параграфе 2, может быть сформулирована как вариационная задача. Рассмотрим множество пар непрерывных функций
О = {(а,Р)е С (О) х С (Б)|
Р (х) -а(и) < р(и, х) У и е О, Ух е Б }.
(16)
Определим на О следующий функционал
I(а,Р)= |р(х)Е(х)йх- ]а(и)Е0(и)&и. (17)
Б О
Очевидно, что I (а, Р) является непрерывным относительно нормы ||(а, Р)||п = шах{||а||ю ||Р||Л, где норма || - ||к определяется как максимум модуля функции по всей области.
Теорема 1. Следующие условия эквивалентны:
1) Пара (Ф, Т) является слабым решением.
2) Пара (Ф, Т) максимизирует функционал (17).
Доказательство теоремы приведено в приложении 2.
В дискретном варианте задача максимизации
функционала (17) может быть записана как задача линейного программирования (ЗЛП). Действительно, возьмём наборы точек {и, е О}г=1"^ и {хj е Б}-=1М,
являющиеся узлами прямоугольных решёток Л и © на О и Б. Обозначим е0 и ej средние значения функций Е0(и) и Е(х) в ячейках решёток. Кроме того, обозначим р,- = р (щ, х,). В результате е0, ej, р,- - известные числа, а а, = а(и,), р,- = Р(х,) - неизвестные переменные.
В дискретном варианте задача максимизации функционала (17) на множестве О сводится к следующей ЗЛП:
М N
IЛ'© (а,,р,) = ХР]е] | т,- | -Ха,е0 | т, шах,
р=1 ,=1 (18)
Р- - а, < р,, , = 1,...,N, р = 1,...,М, где |т,- |, |т,| - площади соответствующих ячеек решёток. Отметим, что функция 1Л©(а,, р,) в (18) соответствует интегральной сумме для непрерывного функционала (17). Естественно ожидать, что при измельчении решёток решения дискретной задачи (18) сходятся к решению исходной непрерывной задачи. Отметим, что матрица ограничений для ЗЛП (18) имеет специфический вид: элементы матрицы принимают значения 0 или ±1. Для решения таких ЗЛП существуют алгоритмы с полиномиальной сложностью [16].
4. Связь с задачей перемещения масс
Рассмотрим связь вариационной задачи максимизации функционала (17) с транспортной задачей пе-
ремещения массы из области О в область Б с функцией стоимости р (и, х). При этом распределения масс в областях описываются функциями Е0(и), иеО и Е(х), хеБ.
Планом будем называть отображение Р: О^Б, сохраняющее энергию (меру), то есть для любого измеримого подмножества ю с О выполнено
J E0(u )du = Je (x) dx.
P-(<o) Ю
(19)
Допустимо, чтобы отображение Р(и) было определено почти всюду. Определим функционал на множестве планов
C (P) = Jp(u, P(u))E0(u)du .
(20)
Теорема 2. Пусть (Ф, Т) — слабое решение проблемы. Тогда отображение уФ минимизирует функционал (20). Более того, если Р - другой оптимальный план, то он совпадает с уФ почти всюду.
Доказательство теоремы приведено в приложении 3.
Теорема 2 представляет интересный результат, который может быть рассмотрен как обобщение принципа Ферма. Действительно, среди всех отображений О^-Б, сохраняющих энергию, слабое решение обратной задачи расчёта эйконала Ф(и) соответствует отображению, для которого минимизируется суммарное расстояние с весом Е0(и) между точками областей О и Б.
Рассмотрим расчёт функции эйконала Ф(и) по отображению х = уФ(и). Для реализации отображения уФ(и) частные производные функции Ф(и), определяющие направление лучей, должны иметь вид:
УФ =
((, фи2 )
p(u, x)
(21)
Функция эйконала может быть восстановлена по своим частным производным на основе следующей известной формулы:
ф( И1, и2 )=Ф( 0,0 ) +
и1 и2 . (22) + (/,0)<1/ +|Фи2 (щ^)<к 0 0 Действительно, в силу определения слабого решения (см. параграф 2) функция Ф(и), иеО является липшицевой. Несложно показать, что функция Ф(и) будет также липшицевой вдоль любой кривой и(т) с натуральной параметризацией. Таким образом, функция Е(т) = Ф(и(т)) является абсолютно непрерывной функцией натурального параметра т. В силу свойств абсолютно непрерывных функций функция Е(т) однозначно восстанавливается по своей производной, существующей почти всюду, интегрированием вдоль кривой. В частном случае, когда кривая и(т) является ломаной со звеньями, параллельными осям координат, приходим к формуле восстановления (22).
В дискретном варианте задача перемещения масс может быть сформулирована как линейная задача о
назначениях (ЛЗН) [19]. Действительно, предположим, что области О и Б разбиты на N ячеек (или аппроксимированы N ячейками) с одинаковой энергией, то есть для любой пары ячеек gi с О, ^ с Б выполняется равенство
JE0 (u) du = JE(x) dx = e.
В этом случае все отображения G ^ D, сохраняющие энергию, можно описать перестановками из N чисел (i1, ..., iN), которые определяют, в какие ячейки dj области D отображаются ячейки g1, ..., gN. В рамках ЛЗН отображение ячейки gi в dj можно интерпретировать как выполнение i-м работником j-й работы. Согласно (21), стоимость работ описывается матрицей Cy = р (ui, xj), i,j = 1, ..., N, где u, xj - центры ячеек gi в dj. Таким образом, для численного построения отображения уФ, соответствующего слабому решению, можно решать следующую задачу о назначениях
C(j1,...., jN) = Xp(u/, xji) ^ min, (23)
i
где минимум вычисляется по всем перестановкам j1, ...,jN. Для решения ЛЗН предложены эффективные алгоритмы, решающие задачу за полиномиальное время [19].
5. Расчётные примеры
Рассмотрим задачу расчёта функции эйконала из условия фокусировки пучка круглого сечения с постоянной освещённостью в прямоугольник с постоянной освещённостью. Расчёт проведём для следующих параметров: радиус пучка (радиус области G) R = 3 мм, расстояние до плоскости фокусировки f= 50 мм, размер сторон прямоугольника (области D) dx1 = 12 мм, dx2 = 4 мм. В силу симметрии задачи, функцию эйконала достаточно рассчитать в первом квадранте из условия фокусировки в прямоугольник D1 6^2 мм2, являющийся одной четвёртой частью прямоугольника Dy расположенной в первом квадранте. Для записи ЗЛП (19) были взяты прямоугольные решётки Л и © на G и D с размерами ячеек 0,23x0,23 мм2 и 0,375x0,222 мм2. В этом случае функция Ф^) в первом квадранте и функция T(x) в прямоугольнике D1 представляются в узлах решёток наборами из N=M= 144 значений. ЗЛП в данном случае содержит 288 переменных и N2 = 20736 ограничений-неравенств. Для решения ЗЛП использовалась программа linprog.m из библиотеки Matlab. Рассчитанная функция эйконала Ф (u) приведена на рис. 2.
По эйконалу Ф (u) может быть рассчитан преломляющий оптический элемент, расположенный непосредственно перед плоскостью z = 0 и формирующий данное распределение эйконала в области G. Рассмотрим плоский освещающий пучок, распространяющийся в направлении оси z. В качестве первой поверхности оптического элемента будем использовать плоскую поверхность, перпендикулярную падающему пучку. Вторая поверхность может быть рассчитана из условия формирования заданного эй-
конала в области О плоскости г = 0 с использованием принципа Ферма [14]. В то же время при указанных «параксиальных параметрах» высота преломляющей поверхности может быть рассчитана в приближении тонкого оптического элемента [13]:
г(и) = Ф(и)/(п -1), (24)
где п - показатель преломления материала оптического элемента.
Ui, мм
О О
U2, ММ
Рис. 2. Функция эйконала в первом квадранте, рассчитанная на основе решения ЗЛП
Для проверки правильности полученных в работе теоретических результатов, по функции эйконала Ф(и) на рис. 2 была рассчитана преломляющая поверхность (24) и затем было проведено моделирование работы полученного оптического элемента в программе для светотехнических расчётов TracePro® с использованием метода трассировки лучей [20]. Для этого поверхность (24) была аппроксимирована системой неоднородных рациональных сплайнов Безье (NURBS) в программе автоматизированного проектирования Rhinoceros® [21]. На рис. 3 показаны преломляющий оптический элемент (n = 1,5) и рассчитанное в программе TracePro распределение освещённости, которое формируется при освещении оптического элемента коллимированным пучком с постоянной освещённостью.
Рис. 3 а, б показывают фокусировку в прямоугольник заданного размера. Среднеквадратическое отклонение (СКО) полученного распределения освещённости от постоянного значения составляет 16,2 %. Неравномерность распределения энергии в точках фокусировки объясняется недостаточным числом точек N = 144, использованным для представления функции эйконала (рис. 2). В то же время увеличение числа точек N приводит к сильному увеличению вычислительной сложности ЗЛП, поскольку размерность матрицы ограничений, определяющая сложность задачи, равна N2. Время решения ЗЛП для N = 144 на стандартном персональном компьютере (Intel® Core™ i7-2600 CPU @3.4 GHz) составило около 5 часов.
Существенно более эффективный метод решения задачи состоит в расчёте отображения х=уФ(м) с использованием ЛЗН (24). Для расчёта отображения уФ в рамках той же задачи фокусировки апертура G была аппроксимирована набором из N = 1060 квадратных ячеек gi с размером 0,162*0,162 мм2. Прямоугольник D также был представлен в виде объединения N = 1060 прямоугольных ячеек d с размером 0,231*0,21 мм2.
4 6
Хь Х2, ММ
Рис. 3. Оптический элемент для фокусировки пучка круглого сечения в прямоугольник 12^4 мм2, рассчитанный
на основе решения ЗЛП. Формируемое распределение освещённости (а) и сечения распределения освещённости по осям координат (б)
Для решения ЛЗН использовалась программа типкгеБ.ш [22], реализующая венгерский алгоритм в среде МаНаЪ [19]. В этом случае время решения ЛЗН при прочих равных условиях составило менее одной секунды. По рассчитанному отображению по формуле (21) были вычислены производные функции эйконала, и затем по формуле (22) была восстановлена функция эйконала численным интегрированием. Далее с использованием (24) был построен оптический элемент, фокусирующий в прямоугольник. Полученный оптический элемент и рассчитанное в программе ТгасеРго распределение освещённости, формируемое оптическим элементом, приведено на рис. 4. В отличие от рис. 3, распределение освещённости на рис. 4 является значительно более равномерным, СКО полученного распределения освещённости от постоянного значения составляет всего 5,6 %.
Представленные примеры показывают, что расчёт функции эйконала на основе решения ЛЗН (23) является оптимальным как с точки зрения времени расчёта, так и качества формируемого распределения.
Заключение
В работе показано, что задача расчёта функции эйконала из условия фокусировки в заданную область может быть сформулирована как вариационная задача максимизации функционала (17) и как задача Монжа-Канторовича о перемещении масс. Получено, что функция стоимости в задаче Монжа-Канторовича соответствует расстоянию между точкой области задания эйконала и точкой области фокусировки.
4 6
XI, Х2, ММ
Рис. 4. Оптический элемент для фокусировки пучка круглого сечения в прямоугольник 12*4 мм2, рассчитанный
на основе решения ЛЗН. Формируемое распределение освещённости (а) и сечения распределения освещённости по осям координат (б)
Данный результат показывает, что решение обратной задачи расчёта эйконала соответствует отображению, для которого минимизируется суммарное расстояние между точками исходной плоскости и области фокусировки. Предложенный в работе формализм позволил свести расчёт функции эйконала к решению ЗЛП (18). При этом расчёт отображения, соответствующего функции эйконала, сводится к решению ЛЗН (23). Представленные примеры расчёта оптических элементов показывают, что метод расчёта эйконала на основе решения ЛЗН (23) требует существенно меньших вычислительных затрат по сравнению с методом на основе решения ЗЛП (18).
Отметим, что подход, основанный на использовании ЛЗН, может быть также использован для расчёта зеркал, формирующих заданные распределения интенсивности. Отличие состоит только в виде функции стоимости [11, 15]. С точки зрения вычислительных затрат данный подход является существенно более эффективным, чем расчёт зеркал на основе решения ЗЛП, рассмотренный в [17, 18].
Благодарности Работа выполнена при поддержке гранта РФФИ 1807-00982 в части формулировки и решения указанной задачи как задачи о перемещении масс с использованием линейной задачи о назначениях и при поддержке Федерального агентства научных организаций (соглашение № 007-ГЗ/Ч3363/26) в части вариационной формулировки обратной задачи расчёта функции эйконала из условия фокусировки в заданную область.
Литература
1. Wu, R. Freeform illumination design: a nonlinear boundary problem for the elliptic Monge-Ampere equation / R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, X. Liu // Optics Letters. - 2013. - Vol. 38, Issue 2. - P. 229-231. - DOI: 10.1364/OL.38.000229.
2. Wu, R. Influence of the characteristics of a light source and target on the Monge-Ampere equation method in freeform optics design / R. Wu, P. Benitez, Z. Yaqin, J.C. Minano // Optics Letters. - 2014. - Vol. 39, Issue 3. - P. 634-637. -DOI: 10.1364/OL.39.000634.
3. Ma, Y. Hybrid method of free-form lens design for arbitrary illumination target / Y. Ma, H. Zhang, Z. Su, Y. He, L. Xu, X. Lui, H. Li // Applied Optics. - 2015. - Vol. 54, Issue 14. -P. 4503-4508. - DOI: 10.1364/AO.54.004503.
4. Bosel, C. Ray mapping approach for the efficient design of continuous freeform surfaces / C. Bosel, H. Gross // Optics Express. - 2016. - Vol. 24, Issue 13. - P. 14271-14282. -DOI: 10.1364ЮЕ.24.014271.
5. Oliker, V.I. Mathematical aspects of design of beam shaping surfaces in geometrical optics / V.I. Oliker. - In Book: Trends in nonlinear analysis / V.I. Oliker, M. Kirkilionis, S. Kromker, R. Rannacher, F. Tomi. - Chapter 4. - Berlin, Heidelberg: Springer, 2003. - P. 193-224. - DOI: 10.1007/978-3-662-05281-5_4.
6. Kochengin, S.A. Computational algorithms for constructing reflectors / S.A. Kochengin, V.I. Oliker // Computing and Visualization in Science. - 2003. - Vol. 6, Issue 1. - P. 1521. - DOI: 10.1007/s00791-003-0103-2.
7. Oliker, V. Controlling light with freeform multifocal lens designed with supporting quadric method (SQM) / V. Oliker // Optics Express. - 2017. - Vol. 25, Issue 4. - P. A58-A72. -DOI: 10.1364/OE.25.000A58.
8. Oliker, V. Supporting quadric method in optical design of freeform lenses for illumination control of a collimated light / V. Oliker, J. Rubinstein, G. Wolansky // Advances in Applied Mathematics - 2015. - Vol. 62 - P. 160-183. - DOI: 10.1016/j.aam.2014.09.009.
9. Fournier, F.R. Fast freeform reflector generation using source-target maps / F.R. Fournier, W.J. Cassarly, J.P. Rolland // Optics Express. - 2010. - Vol. 18, Issue 5. - P. 52955304. - DOI: 10.1364/OE.18.005295.
10. Michaelis, D. Cartesian oval representation of freeform optics in illumination systems / D. Michaelis, P. Schreiber, A. Bauer // Optics Letters. - 2011. - Vol. 36, Issue 6. -P. 918-920. - DOI: 10.1364/OL.36.000918.
11. Glimm, T. Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem / T. Glimm, V. Oliker // Journal of Mathematical Sciences. - 2003. -Vol. 117, Issue 3. - P. 4096-4108. - DOI: 10.1023/A:1024856201493.
12. Doskolovich, L.L. On the use of the supporting quadric method in the problem of the light field eikonal calculation / L.L. Doskolovich, M.A. Moiseev, E.A. Bezus, V. Oliker // Optics Express. - 2015. - Vol. 23, Issue 15. - P. 1960519617. - DOI: 10.1364/OE.23.019605.
13. Soifer, V. Iterative methods for diffractive optical elements computation / V. Soifer, V. Kotlyar, L. Doskolovich. -London: Taylor & Francis Ltd., 1997. - 244 p. - ISBN: 07484-0634-4.
14. Doskolovich, L.L. Analytic design of optical elements generating a line focus / L.L Doskolovich, A.Yu. Dmitriev, S.I. Kharitonov // Optical Engineering. - 2013. - Vol. 52, Issue 13. - 091707. - DOI: 10.1117/1.OE.52.9.091707.
15. Wang, X.J. On the design of a reflector antenna II / X.J. Wang // Calculus of Variations and Partial Differential Equations. - 2004. - Vol. 20, Issue 3. - P. 329-341. - DOI: 10.1007/s00526-003-0239-4.
16. Tardos, E. A strongly polynomial algorithm to solve combinatorial linear programs / Operation Research. - 1986. -Vol. 34, Issue 2. - P. 250-256. - DOI: 10.1287/opre.34.2.250.
17. Canavesi, C. Direct calculation algorithm for two-dimensional reflector design / C. Canavesi, W.J. Cassarly, J.P. Rolland // Optics Letters. - 2012. - Vol. 37, Issue 18. -P. 3852-3854. - DOI: 10.1364/OL.37.003852.
18. Canavesi, C. Observations on the linear programming formulation of the single reflector design problem / C. Canavesi, W.J. Cassarly, J.P. Rolland // Optics Express. - 2012. - Vol. 20, Issue 4. - P. 4050-4055. - DOI: 10.1364/OE.20.004050.
19. Munkres, J. Algorithms for the assignment and transportation problems / J. Munkres // Journal of the Society for Industrial and Applied Mathematics. - 1957. - Vol. 5, No. 1. -P. 32-38. - DOI: 10.1137/0105003.
20. Trace Pro [Electronical Resource]. - URL: https://www.lambdares.com/tracepro/ (date request 12.03.2018).
21. Rhinoceros®: design, model, present, analyze, realize... [Electronical Resource]. - URL: http://www.rhino3d.com (date request 12.03.2018).
22. MathWorks. File exchange. Munkres assignment algorithm [Electronical Resource]. - URL: https://www.mathworks.com/matlabcentral/fileexchange/20 328-munkres-assignment-algorithm (date request 12.03.2018).
(25)
Приложение 1 Лемма. Пусть (а(и), Р(х)) - хорошая пара. Тогда функция а(и) (функция Р(х)) является липшицевой на О (на Б), то есть существует константа Ь, для которой верно:
|а(и') -а(и ")|< Ь • |и'- и "|
для любых и', и'' е О.
Возьмём произвольные и', и''еО. Будем полагать, что а(и'' ) > а(и') (случай а(и'' ) > а(и') рассматривается аналогично). Пусть х е Б - точка, для которой выполняется равенство Р(х'' ) - а(и'' ) = р(и '' , х'' ). Тогда а(и'' ) -а(и' ) =
= Р (х'' ) -р (и'' , х'' ) - а (и' ) <р (и' , х'' ) -р (и'' , х'' ).
Согласно формуле конечных приращений Ла-гранжа, правую часть последнего неравенства можно записать в виде следующего скалярного произведения:
р(и', х '' ) - р(и'', х'' ) = Ур (и, х'' )-(и'- и") , (26)
где й е О - некоторая точка, лежащая на отрезке, соединяющем точки и и и , а градиент функции р(и, х) вычисляется по переменным (иь и2). Так как левая часть неравенства (25) неотрицательна, её модуль с учётом (26) и неравенства Коши-Шварца может быть оценен как
|а(и'' ) - а(и')| < |Ур(и, х'')| • |и'- и "| < Ь • |и'- и "|,
где
L = max
ueG,xeD
|Ур(м, х)| .
Таким образом, условие липшицевости выполняется с указанной константой Ь. Доказательство лип-шицевости для функции Р(х) аналогично.
Приложение 2
Доказательство теоремы 1 Докажем сначала, что из условия 1 следует условие 2. Пусть (Ф, Т) - слабое решение проблемы. Кроме того, возьмём произвольную пару (а, Р)еО. В частности,
Р(УФ (и)) - а(и) < р(и, уф (и)) = Т(уф (и)) - Ф(и).
Умножим левую и правую часть неравенства на Е0 (и) и проинтегрируем по В результате мы получим
|Р(уФ(u))E0(u)du - Ja(u)E0(u)du <
G G
< JV(уф (u))E0 (u)du - JO(u)E0 (u)du.
(27)
Поскольку отображение уФ сохраняет энергию, то из формулы (4) в основном тексте статьи получим следующие равенства:
J*P (уф (u)) E0(u)du = Jp (х) E (x)dx,
G D
J^(yo (u ))E0(u)du = JT( x) E (x)dx.
(28)
О Б
В результате неравенство (27) преобразуется к виду I (а, Р) < 1(Ф, Т), где I (а, Р) имеет вид (17). Это и означает, что максимальное значение функционал I (а, Р) принимает на паре (Ф, Т).
Докажем, что из условия 2 следует условие 1. Сначала докажем, что экстремальная пара (Ф, Т), максимизирующая значение функционала, является «хорошей», то есть выполнены соотношения
4>(u) = max{^ (x) - р (u, x)},
xeD
x) = min{0(u) +р (u, x)}.
ugG
(29)
Докажем это от противного. Если не выполнено первое соотношение в (29), то определим функцию Ф (и) в виде
Ф ' (и ) = шах{Т( х) -р (и, х)}. (30)
хеБ
Пара (Ф ', Т) также лежит в О по определению Ф ' (и). Возьмем произвольную точку и0еО и обозначим через х0еБ точку, для которой Т(х0) - Ф' (и0) = р (и0, х0). Теперь
Ф '(и0) = Т(х0) -р (и0, х0) < Ф(и0). (31)
Поскольку Е0 (и) > 0, то согласно (31)
Jo '(u)E0(u)du < JO(u)E0(u)du .
Из полученного неравенства следует, что I (Ф', Т) >/(Ф, Т). Причем если Ф'(и) ф Ф(и) на множестве положительной меры, то I(Ф', Т) > 1(Ф, Т), то есть пара (Ф, Т) не экстремальна, что противоречит предположению. Аналогичные рассуждения можно провести для Т(х) и убедиться в выполнении и второго соотношения в (29).
Теперь докажем выполнение условия сохранения энергии для отображения уФ, определённого уравнением (14) с функциями Ф(и), Т(х), максимизирующими функционал I.
Согласно необходимому условию (Эйлера-Лаг-ранжа) экстремума функционала 1(у) в точке экстремума первая вариация этого функционала равна нулю, то есть
Т" I ()
de
= 0 для у8 —> у при 8 — 0.
8=0
Мы выберем семейство функций к8, стремящихся к нулю по норме при 8—0, и покажем, что из условия Эйлера-Лагранжа (когда в качестве у взята экстремальная пара (Ф, Т)) следует условие сохранения энергии для отображения уФ). Пусть ю с Б - измеримое подмножество, хю - индикатор множества ю. Выберем последовательность непрерывных функций %„, поточечно сходящихся почти всюду к функции хю (существование такой последовательности следует из теоремы Фреше). Пусть (Ф, Т) - экстремальная пара. Возьмём достаточно малое 8 > 0 и определим «деформированную» пару (Ф 8, Т8) следующим образом
Т8 (х) = Т (х) + 8Х п (х),
Ф8 (и) = тах{Т8 (х) - р (и, х)} •
хеБ
(32)
Можно показать, что функция Ф8(и) непрерывна на О. Неравенство Т8(х) - Ф8(и) < р (и, х) выполнено по определению Т8(х), Ф8(и). Таким образом, пара (Ф8 , Т8) лежит в О. Зафиксируем точку и е О. Ближайшая цель - вычислить производную по 8 от функции Ф8(и). Пусть х8еБ - точка, для которой достигается равенство Т8(х 8) - Ф8(и) = р(и, х8), то есть точка, в которой достигается максимум из определения Ф8(и) в (2.7). Тогда Ф8 (и) -Ф(и) =
= Т8 (х8) -р (и, х8) -Ф(и) < (33)
<Т8 (х8)-Т( х8 ) = 8Х п (х8). Теперь пусть х = уФ(и), тогда Ф8 (и) -Ф(и) = = Ф8 (и) - Т( х) + р(и, х) = = Ф8 (и) - (Т(х) + 8Хп (х)) + р (и, х) +8Хп (х) — — 8Хп (х).
Используя (33), (34), запишем итоговое неравенство
8Хп (х) <Ф8(и)-Ф(и) <8Хп (хв). (35)
(34)
Вычтем из обеих частей (35) 8%п(х) и разделим на 8 > 0:
0 < Ф8 (и) -Ф(и) - 8Х п (х) < 8(хп (х8 ) - хп (х)) , Ф8 (и) -Ф(и)
--Хп (х)
< |Хп (х8)-Хп (х)| .
Заметим, что х8—х при 8—0, поэтому d
d8
Ф8 (Уф1(х))
= Х п (х) = х п (уф (и)).
(36)
8=0
Поскольку функционал имеет максимум в точке (Ф, Т), то выполнено условие Эйлера-Лагранжа, то есть
0 = АI (Ф8 ,Т8 )
d8 8=0 (37)
= |х п (х)Е (х) dx - |х п (У (и))Е0 (и )du.
Б О
Теперь переходя к пределу при п — ж, получаем условие сохранения энергии для отображения уФ:
\Е(х)6х = | Е0(и)du . (38)
ю Уф1 (ю)
Это завершает доказательство теоремы.
Приложение 3 Доказательство теоремы 2
Пусть Р - произвольный план. Имеет место неравенство р (и, Р(и)) > Т(Р(и)) - Ф(и), причем равенство достигается тогда и только тогда, когда Р (и) = уФ(и). Умножим это неравенство на Е0 (и) и проинтегрируем по du:
С (Р) = |р(и, Р(и)) Е0(и )du —
г О г (39)
— |Т (Р(и)) Е0 (u)du - |Ф(и) Е0 (и )du.
О О
Согласно формуле (4) основного текста статьи будем иметь
|Т( Р(и))Е0(и^и = |Т (х) Е (х)6х =
О Б
= |Т (уф (и)) Е0(и^и.
О
А значит,
|Т( Р(и)) Е0 (и) du - |ф(и) Е0 (u)du =
(40)
= |Т( уф (и)) Е0 (и )du - |ф(и) Е0 (и) du =
О О
= |(Т(уф (и)) -ф(и))Eo(u)du =
О
= |р(и, Уф (и ))Eo(u)du = С (уф).
Сравнивая (39), (40), получаем С (Р) > С (уФ), при этом, если Р (и) ф уФ(и) на множестве положительной меры, то неравенство строгое.
а
Сведения об авторах
Досколович Леонид Леонидович, в 1989 году с отличием окончил Куйбышевский авиационный институт (КуАИ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Прикладная математика». Доктор физико-математических наук (2001 год), профессор, работает заведующим лабораторией дифракционной оптики Института систем обработки изображений РАН (ИСОИ РАН) - филиала ФНИЦ «Кристаллография и фотоника» РАН, профессором кафедры технической кибернетики Самарского университета и ведущим научным сотрудником научно-исследовательской лаборатории прорывных технологий дистанционного зондирования Земли Самарского университета. Специалист в области дифракционной оптики, лазерных информационных технологий, нанофотоники. E-mail: leonid@,smr.ru .
Мингазов Альберт Айдарович, в 2010 году окончил Самарский государственный университет по специальности «Математика». В 2012 году окончил магистратуру НИУ ВШЭ по направлению «Математика», в 2015 окончил аспирантуру Санкт-Петербургского отделения математического института имени В.А. Стеклова. Кандидат физико-математических наук (2015), научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН. E-mail: mingazov88@gmail.com .
Быков Дмитрий Александрович, в 2009 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (СГАУ) по специальности «Прикладная математика и информатика». Доктор физико-математических наук (2017 г.), старший научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН (ИСОИ РАН) - филиала ФНИЦ «Кристаллография и фотоника» РАН и доцент Самарского университета. Области научных интересов: оптика резонансных дифракционных структур, магнитооптика, электромагнитная теория дифракции. E-mail: bykovd@gmail.com .
Андреев Евгений Сергеевич, в 2014 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (СГАУ) по направлению «Прикладные математика и физика». Аспирант по направлению 03.06.01 «Физика и астрономия», специальность «Оптика», Самарский университет. E-mail: gsomix@gmail.com .
ГРНТИ: 29.31.29
Поступила в редакцию 12 марта 2018 г. Окончательный вариант - 13 июня 2018 г.
VARIATIONAL APPROACH TO EIKONAL FUNCTION COMPUTATION
L.L. Doskolovich 1,2, A.A. Mingazov1, D.A. Bykovl2, E.S. Andreev1,2
1 Samara National Research University, 34, Moskovskoye shosse, Samara, 443086, Samara, Russia,
2IPSIRAS - Branch of the FSRC "Crystallography and Photonics" RAS, Molodogvardeyskaya 151, 443001, Samara, Russia
Abstract
The problem of calculating the eikonal function from the condition of focusing into a prescribed region is formulated as a variational problem and as a Monge-Kantorovich mass transportation problem. It is found that the cost function in the Monge-Kantorovich problem corresponds to the distance between a point of the original region (in which the eikonal function is defined) and a point of the focal region. The formalism proposed in this work makes it possible to reduce the calculation of the eikonal function to a linear programming problem. In this case, the calculation of the "ray mapping" corresponding to the eikonal function is reduced to the solution of a linear assignment problem. The proposed variational approaches are illustrated by examples of calculation of optical elements for focusing a circular beam into a rectangular region.
Keywords: geometrical optics, eikonal function, variational problem, Monge-Kantorovich mass transportation problem.
Citation: Doskolovich LL, Mingazov AA, Bykov DA, Andreev ES. Variational approach to eikonal function computation. Computer Optics 2018; 42(4): 557-567. DOI: 10.18287/2412-61792018-42-4-557-567.
Acknowledgements: The work was partially funded by the Russian Science Foundation (RSF) under grant #18-07-00982 (problem formulation and solution as a Monge-Kantorovich mass transportation problem) and by the Federal Agency for Scientific Organizations (FASO) under grant 007-ГЗ/Ч3363/26 (formulation of the problem of calculating the eikonal function from the condition of focusing into a prescribed region as a variational problem).
References
[1] Wu R, Xu L, Liu P, Zhang Y, Zheng Z, Li H, Liu H. Freeform illumination design: a nonlinear boundary problem for the elliptic Monge-Ampére equation. Opt Lett 2013; 38(2): 229-231. DOI: 10.1364/OL.38.000229.
[2] Wu R, Benitez P, Yaqin Z, Minano JC. Influence of the characteristics of a light source and target on the Monge-Ampere equation method in freeform optics design. Opt Lett 2014; 39(3): 634-637. DOI: 10.1364/OL.39.000634.
[3] Ma Y, Zhang H, Su Z, He Y, Xu L, Lui X, Li H. Hybrid method of free-form lens design for arbitrary illumination target. Appl Opt 2015; 54(14): 4503-4508. DOI: I0.1364/A0.54.004503.
[4] Bösel C, Gross H. Ray mapping approach for the efficient design of continuous freeform surfaces. Opt Express 2016; 24(13): 14271-14282. DOI: 10.1364/0E.24.014271.
[5] Oliker VI. Mathematical aspects of design of beam shaping surfaces in geometrical optics. In Book: Oliker VI, Kirkili-onis M, Krömker S, Rannacher R, Tomi F, eds. Trends in nonlinear analysis. Chap 4. Berlin, Heidelberg: Springer; 2003: 197-224. DOI: 10.1007/978-3-662-05281-5_4.
[6] Kochengin SA, Oliker VI. Computational algorithms for constructing reflectors. Computing and Visualization in Science 2003; 6(1): 15-21. DOI: 10.1007/s00791-003-0103-2.
[7] Oliker V. Controlling light with freeform multifocal lens designed with supporting quadric method (SQM). Opt Express 2017; 25(4): A58-A72. DOI: 10.1364/OE.25.000A58.
[8] Oliker V, Rubinstein J, Wolansky G. Supporting quadric method in optical design of freeform lenses for illumination control of a collimated light. Advances in Applied Mathematics 2015; 62: 160-183. DOI: 10.1016/j.aam.2014.09.009.
[9] Fournier FR, Cassarly WJ, Rolland JP. Fast freeform reflector generation using source-target maps. Opt Express 2010; 18(5): 5295-5304. DOI: 10.1364/OE.18.005295.
[10] Michaelis D, Schreiber P, Bäuer A. Cartesian oval representation of freeform optics in illumination systems. Opt Lett 2011; 36(6): 918-920. DOI: 10.1364/OL.36.000918.
[11] Glimm T, Oliker V. Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem. J Math Sci 2003; 117(3): 4096-4108. DOI: 10.1023/A:1024856201493.
[12] Doskolovich LL, Moiseev MA, Bezus EA, Oliker V. On the use of the supporting quadric method in the problem of the light field eikonal calculation. Opt Express 2015; 23(15): 19605-19617. DOI: 10.1364/OE.23.019605.
[13] Soifer V, Kotlyar V, Doskolovich L. Iterative methods for diffractive optical elements computation. London: Taylor & Francis Ltd; 1997. ISBN: 0-7484-0634-4.
[14] Doskolovich LL, Dmitriev AYu, Kharitonov SI. Analytic design of optical elements generating a line focus. Opt Eng 2013; 52(13): 091707. DOI: 10.1117/1.OE.52.9.091707.
[15] Wang XJ. On the design of a reflector antenna II. Calculus of Variations and Partial Differential Equations 2004; 20(3): 329-341. DOI: 10.1007/s00526-003-0239-4.
[16] Tardos E. A strongly polynomial algorithm to solve combinatorial linear programs. Operation Research 1986; 34(2): 250-256. DOI: 10.1287/opre.34.2.250.
[17] Canavesi C, Cassarly WJ, Rolland JP. Direct calculation algorithm for two-dimensional reflector design. Opt Lett 2012; 37(18): 3852-3854. DOI: 10.1364/OL.37.003852.
[18] Canavesi C, Cassarly WJ, Rolland JP. Observations on the linear programming formulation of the single reflector design problem. Opt Express 2012; 20(4): 4050-4055. DOI: 10.1364/OE.20.004050.
[19] Munkres J. Algorithms for the assignment and transportation problems. Journal of the Society for Industrial and Applied Mathematics 1957; 5(1): 32-38. DOI: 10.1137/0105003.
[20] Trace Pro. Source: (https://www.lambdares.com/tracepro/).
[21] Rhinoceros®: design, model, present, analyze, realize... Source: (http://www.rhino3d.com).
[22] MathWorks. File exchange. Munkres assignment algorithm. Source: (https://www.mathworks.com/matlabcentral/ fileexchange/20328-munkres-assignment-algorithm).
Author's information
Leonid Leonidovich Doskolovich graduated with honours (1989) from S.P. Korolyov Kuibyshev Aviation Institute (presently, Samara National Reseach University), majoring in Applied Mathematics. He received his Doctor in Physics & Maths (2001) degree from Samara State Aerospace University. Laboratory chief in the Image Processing Systems Institute RAS - Branch of the FSRC "Crystallography and Photonics" RAS, professor at Technical Cybernetics department of National Research University, the senior researcher at the Breakthrough Technologies for Earth's Remote Sensing laboratory at SSAU. His leading research interests include diffractive optics, laser information technologies, nanophotonics. E-mail: leonid@smr.ru .
Albert Aidarovich Mingazov graduated from Samara State University in 2010, majoring in Mathematics. In 2012 he got the master degree in Mathematics at National Research University Higher School of Economics, in 2015 he graduated from the post-graduate course of the St. Peterburg Department of Steklov Mathematical Institute of Russian Academy of Sciences. Candidate in Phisics and Mathematics (2015). Currently he is a researcher at Diffractive Optics laboratory of the Image Processing Systems Institute RAS - Branch of the FSRC "Crystallography and Photonics" RAS. E-mail: mingazov88@smail.com .
Dmitry Alexandrovich Bykov graduated with honors (2009) from Samara State Aerospace University (SSAU), majoring in Applied Mathematics and Computer Science. Doctor of Physics and Mathematics (2017). Senior researcher at the Samara University and at the Diffractive Optics Laboratory of the Image Processing Systems Institute RAS - Branch of the FSRC "Crystallography and Photonics" RAS. His current research interests include optics of resonant diffractive structures, magneto-optics of nanostructured materials, and electromagnetic diffraction theory. Email: bykovd@gmail.com .
Evgeniy Sergeevich Andreev graduated with honors (2014) from Samara State Aerospace University (SSAU), majoring in Applied Mathematics and Physics. Graduate student majoring in Optics at Samara University. E-mail:
gsomix@gmail.com .
Received March 12, 2018. The final version - June 13, 2018.
