На днях задумался над некоторыми проблемами геометрической оптики со стороны инженерных задач. Литература по расчёту оптических приборов в основном апеллирует к центрированным системам с аксиальной симметрией [1, 2], хотя, как мне кажется, ценой этого упрощения теряется изрядная масса творческой инженерной свободы. Я захотел взглянуть, как можно описать оптичекие процессы в геометрическом приближении с некоторой более общей точки зрения: меня ведь не страшит объём вычислений и аналитическая сложность формул появляющихся при таком подходе.

Математическое описание

Для начала рассмотрим преобразование Хаусхолдера (с матрицей отражения, и да, это точно фамилия, а не плохой перевод англоязычной номенклатуры), чтобы понять, как вообще аналитически следует описать процесс отражения лучей с точки зрения геометрической оптики (здесь мы не интересуемся эффектами ослабления интенсивности луча или его волновой природой): на некоторой неявно заданной $S(\vec{r}) = 0$ поверхности направление отражённых лучей $\hat{d^\mathcal{s}}$ зависит от направления падающих лучей $\hat{d^\mathcal{i}}$ следующим образом: $$\hat{d^\mathcal{s}} = \mathbf{R} \cdot \hat{d^\mathcal{i}},$$ где матрицу Хаусхолдера $\mathbf{R}$ можно записать через градиент (естественно задающий нормаль $\hat{d^\mathcal{n}} = \nabla S(\vec{r})/|\nabla S(\vec{r})|$, [3]): \begin{equation} \mathbf{R} = \mathbf{I} - 2 \hat{d^\mathcal{n}} \hat{d^\mathcal{n}} = \mathbf{I} - \frac{2}{|\nabla S(\vec{r})|^2} \cdot \nabla S(\vec{r}) \cdot \nabla S(\vec{r})^T. \label{eq:refLawGeneric} \end{equation}

Обозначив $$ \mathcal{S}_{ij} = \frac{\partial S(\vec{r})}{\partial r_i} \frac{\partial S(\vec{r})}{\partial r_j} |\nabla S(\vec{r})|^{-2},$$ перепишем уравнение на векторное поле направлений отражённых лучей $\hat{d^{\mathcal{s}}}$, заданное на поверхности зеркала формы $S(\vec{r}) = 0$: \begin{equation} \hat{d^{\mathcal{s}}_j} = \big[ \mathbf{I}_{ij} - 2 \cdot \mathcal{S}_{ij} \big] \hat{d^{\mathcal{i}}_{i}}. \label{eq:refLaw} \end{equation}

Для обуславливания счётной процедуры, помимо, собственно, самой поверхности $S(\vec{r})$ понадобится задать набор её частных производных ($\partial S(\vec{r})/\partial r_i$), которые, впрочем, могут быть расчитаны и численным методом. Для непосредственного счёта же $\hat{d^{\mathcal{s}}_j}$ — нужно определить точку на поверхности и направление падающего луча $\hat{d^{\mathcal{s}}_i}$.

Счёт

Реализовать формулу \eqref{eq:refLaw} можно в виде набора простейших процедур на Си, и на данном этапе возникает соблазн изобразить векторное поле для какого-нибудь простенького в смысле исходных данных, случая: так, например, можно проверить свойства параболоида собирать пучок параллельных лучей в фокусе. Возьмём его в следующем виде:

\begin{equation} S(x, y, z) = 2 - (x^2 + y^2 + z) = 0, \\ \vec{\nabla S} = \{ -2x; -2y; -1 \}. \end{equation}

Легко видеть, что матрица Хаусхолдера $\mathbf{R}$ симметрична, и для её хранения хватит массива из шести элементов, и никакие другие матрицы в этом алгоритме не встречаются. Хотя я хочу пока посмотреть простой центированный параболоид, универсальность метода наверняка ещё позабавит меня, и целесообразно организовать код с возможностью повторного использования. Собираясь писать на чистом Си, я определил структуры для трёхмерного вектора, направления в сферической системе координат и симметричной матрицы. Также, определил тип коллбэка на компоненты градиента (каждый из которых можно определить как скалярное поле), и реализовал несколько функций. Люблю выразительное удобство выбора между именованием векторов по компонентам и целочисленным индексом; подобное совмещение в Си можно организовать через union и вычислительно оно ничего не стоит.

/* 3-dimensional vector */
typedef union {
    struct { Real_t x, y, z; };
    Real_t r[3];
} Vector3;

/* 9-component matrix which is effectively 6-component (symmetrical) */
typedef struct {
    Real_t el[6];
} SymMatrix;

/* Type of function pointer defining scalar field */
typedef Real_t (*ScalarField)( const Vector3 * r );

/* Type of fully-conditioned surface */
typedef union {
    struct { ScalarField dSx, dSy, dSz; };
    ScalarField dR[3];
} Grad;

Функция вычисления градиента в точке:

/* Calculates $v = \vec{nabla V(\vec{r})}$ */
void
nabla_at( const Vector3 * r, const Grad * g
        , Vector3 * v ) {
    for( Index_t i = 0; i < 3; ++i ) {
        v->r[i] = g->dR[i](r);
    }
}

Функцию внешнего произведения (outer product) для градиентов реализуем отмечая в комментариях индексы матрицы чтобы не запутаться:

/* Calculates open product
 * $\nabla V \otimes \nabla V^T$ and returns $|\nabla V|^2$ */
Real_t
oprod_at( const Vector3 * r, const Grad * g
        , SymMatrix * sm) {
    Vector3 n;
    nabla_at( r, g, &n );
    sm->el[0] = n.x*n.x; /*1, 1*/
    sm->el[1] = n.y*n.y; /*2, 2*/
    sm->el[2] = n.z*n.z; /*3, 3*/
    sm->el[3] = n.x*n.y; /*1, 2*/
    sm->el[4] = n.y*n.z; /*2, 3*/
    sm->el[5] = n.x*n.z; /*1, 3*/
    return sm->el[0] + sm->el[1] + sm->el[2];
}

Тогда расчёт матрицы Хаусхолдера:

/* Calculates Householder matrix $R$ at point $\vec{r}$ given by gradient $g$ */
void
householder_m_at( const Vector3 * r, const Grad * g
               , SymMatrix * sm ) {
    Real_t norm = oprod_at( r, g, sm );
    sm->el[0] = 1. - 2*sm->el[0]/norm;
    sm->el[1] = 1. - 2*sm->el[1]/norm;
    sm->el[2] = 1. - 2*sm->el[2]/norm;
    sm->el[3] =    - 2*sm->el[3]/norm;
    sm->el[4] =    - 2*sm->el[4]/norm;
    sm->el[5] =    - 2*sm->el[5]/norm;
}

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

Пока меня не слишком интересует (и не уверен, заинтересует ли достаточно в будущем) raycasting, который потребует проективной сетки лучей задаваемой от плоскости камеры с последующим отысканием точек пересечения луча с поверхностью. Я буду оперировать направлениями выраженными в виде сферических углов, которых достаточно для задания единичных ортов $\hat{d^\mathcal{i}}$, $\hat{d^\mathcal{s}}$. Понятно, что для этого, наверное, придётся воспользоваться несколько избыточными преобразованиями прямоугольных координат в сферические, но пока сгодится:

/* 2-angle direction in spherical space */
typedef union {
    struct { Real_t theta, phi; };
    Real_t angles[2];
} Direction;

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

set pm3d at s
set palette grey
set style line 1 black
set xrange [-3:3]
set yrange [-3:3]
set zrange [-8:2.5]
unset colorbox
set view 65, 115
splot 2 - x**2 - y**2 notitle lt black \
    , "fields.txt" u ($1-$7):($2-$8):($3-$9):7:8:9 w vectors lw 2 title 'Incident' \
    , "fields.txt" u 1:2:3:4:5:6 w vectors lw 2 title 'Scattered'

Результаты визуализации:

Ну и, как водится, полный gist:

Хотя центрированные фигуры вращения с коллинеарным равномерным пучком — в целом очень примитивный результат, думаю, в дальнейшем эти соображения могут для чего-нибудь пригодиться. Представляет, в частности, интерес получить на основе \eqref{eq:refLaw} дифференциальное уравнение, позволяющее для заданных полей падающего и отражённого света, отыскивать гладкую форму отражающей поверхности, хотя это может (а может и нет?) привести к довольно вычурным вещам, вроде функционалов, тензорной алгебры и вариационных методов: всё зависит от того, в каком виде задавать, собственно, падающее световое поле.

Литература

Для будущих размышлений:

• "Методичка от НГУ" быстро и доходчиво излагает идею фазовых портретов дифференциальных систем.
• Например ещё на сайте eqWorld содержится информация о решении различных классов уравнений, в т.ч. и всяких там систем нелинейных ДУ в ЧП.

Ссылки

[1] Математические модели оптических поверхностей при автоматизированном моделировании / С.А. Родионов, А. А. Ашехонин

[2] ПО Кафедры прикладной оптики ИТМО

[3] "Gradients and normals", Kenneth I. Joy