«Физика для программистов» — как физтехи применяют её в приложениях. Дифракция. Интеграл Френеля23.05.2024 01:00

Введение

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

Исходная задача

Как и в прошлый раз, в основу статьи лягут несколько задач на моделирование из программы курса «Цифровизация физических процессов» на факультете ВШПИ МФТИ.

М1 Численно вычислить интеграл Френеля для симметричных отверстий различной формы (эллипс, квадрат, прямоугольник, ромб — по выбору) для точки, расположенной на оси симметрии отверстия на различных расстояниях от отверстия. Выбрать оптимальный способ расчета исходя из формы отверстия.

М2. Вычислить значение интеграла Френеля для точек, смещённых от оси круглого отверстия.

М3. Численно вычислить дифракционную картину по схеме Юнга исходя из спектра падающего света, углового размера источника, геометрических размеров щелей. Предусмотреть возможность добавить входную щель.

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

Теория

Иллюстрация к формуле

Для понимания того, что мы будем дальше делать, стоит немного поговорить о теории. Согласно теории дифракции Рэлея-Зоммерфельда, дифракционная картина электрического поля в точке (x, y, z) определяется следующим решением уравнения Гельмгольца:

$E(x,y,z) = \frac{1}{i\lambda} \int{\int_{-\infty}^{+\infty}{E(x',y',0)\frac{e^{ikr}z}{r^{2}}(1+\frac{i}{kr})dx'}dy'} (1),$

где

E (x', y', 0) — значения поля в точке плоскости «с которой светим»,
λ‎ — длина волны,

$k = \frac{2\pi}{\lambda}, r = \sqrt{(x-x')^{2} + (y-y')^{2} + z^{2}}$

Пусть

$\rho^{2} = (x-x')^{2}+(y-y')^2$

Применим аппроксимацию Френеля:

$r = \sqrt{(x-x')^{2} + (y-y')^{2} + z^{2}} =\sqrt{\rho^{2}+z^{2}} = z\sqrt{1+\frac{\rho^{2}}{z^{2}}} = z(1 + \frac{\rho^{2}}{2z^{2}} + ...) \approx$ $\approx z + \frac{\rho^{2}}{2z}$

Положим также λ << ρ << z , тогда из формулы (1) получим:

$E(x,y,z) = \frac{e^{ikz}}{i\lambda z} \int{\int_{-\infty}^{+\infty}{E(x',y',0)e^{\frac{ik}{2z}\rho^{2}}dx'}dy'} (2)$

Возможно, Вы заметите сходство полученного результата с довольно важным объектом гармонического анализа — преобразованием Фурье, и действительно, этот интеграл может быть довольно легко переписан в таком виде, можете попробовать, это несложно.

Современные вычислительные мощности, конечно, позволяют вычислить и исходный интеграл, однако в данном случае мы будем моделировать именно по формуле (2), чтобы продемонстрировать результат, который выдаёт интеграл Френеля.

Моделирование

Выбор технологий

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

Pipeline программы

Заметим, что из теории следует, что каждая следующая перегородка использует только выходной слой предыдущей для того, чтобы получить картинку за собой, таким образом, появляется естественная идея представить программу в виде некоторого конвейера, который на вход получает input-слой и cover-слой, а на выход выдаёт output-слой — то, что получается на экране.

Ускорение с помощью torch

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

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

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

За параллельные вычисления в решении отвечала библиотека Torch, соответственно, все вычисления происходили в матричном виде. Даже распараллелив только преобразований над точками, удалось достичь ускорения в сотню раз, уменьшив время на один слой примерно до 10 минут (при удовлетворительной дискретизации) вместо ~10–15 часов при полностью последовательном вычислении.

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

def func(self, x, y, z, px, py, c):
  k = 2 * pi / CONFIG["lmbd"]
  kf = exp(complex(0, 1) * k * z) / (complex(0, 1) * CONFIG["lmbd"] * z)
  mult = torch.exp(torch.mul(torch.square(torch_l2(px, py, x, y)), (complex(0, 1) * k / 2 / z)))
  ind = c
  return torch.mul(torch.mul(kf, mult), ind)

Информирование

Для отслеживания прогресса выполнения расчётов использовалась библиотека tqdm.

Примеры использования

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

Везде далее λ = 500 нм.

Круглая щель

На примере этой задачи будет продемонстрирован полный pipeline, далее будут приводиться только результаты на последнем экране.

$r_{\text{щели}} = 0.001\text{м}, z = 2\text{м}$ pipeline приложения на примере круглой щели

pipeline приложения на примере круглой щели

Двухщелевой эксперимент

$d_{между} = 0.005\text{м}, r_{\text{щели}} = 0.001\text{м}, z = 5\text{м}$

Определите расстояние до первого пика. (Ответ d = 0.0005 совпадает с моделью)

Изображения на экранах для различного приближения

Графики значений на срезе

Пятна Пуассона

Пусть

$\frac{d^{2}}{l\lambda} \gtrsim 1,$

где d — диаметр объекта с круглой тенью, l — расстояние до экрана, λ — длина волны, тогда в центре тени, отбрасываемой объектом, должна наблюдаться светлое пятно — пятно Пуассона (или Араго, или Френеля).

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

$r_{\text{диска}} = 0.001\text{м},0.002\text{м},0.004\text{м}, z = 1\text{м}$ Изображения на экранах для различных радиусов

Изображения на экранах для различных радиусов

Графики значений на срезе

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

Заключение

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

Ссылки

GitHub проекта
Статья про маятники
Статья про броски объектов под углом к горизонту

Habrahabr.ru прочитано 4483 раза