diff --git a/M11_report/M11_report.pdf b/M11_report/M11_report.pdf new file mode 100644 index 0000000..a2eaf6a Binary files /dev/null and b/M11_report/M11_report.pdf differ diff --git a/M11_report/conclusion.tex b/M11_report/conclusion.tex new file mode 100644 index 0000000..8d5e0b8 --- /dev/null +++ b/M11_report/conclusion.tex @@ -0,0 +1,34 @@ +\newpage +\section{Заключение} + +В работе реализовано численное моделирование дифракции плоской монохроматической волны на плоских препятствиях по принципу Гюйгенса--Френеля. + +\textbf{Основные результаты:} + +\begin{enumerate} + \item \textbf{Методы вычисления} + \begin{itemize} + \item Реализовано прямое суммирование (точное, но медленное) + \item Реализован метод углового спектра через БПФ (точный, $O(N^2\log N)$) + \item Реализовано приближение Френеля через БПФ для дальней зоны + \end{itemize} + + \item \textbf{Модели препятствий} + \begin{itemize} + \item Щель: аналитически подтверждена формула $\sinc^2$ + \item Круглое отверстие: воспроизведены кольца Эйри и колебания оси при изменении числа зон Френеля + \item Двойная щель: наблюдена интерференционная картина Юнга с правильным периодом + \item Дифракционная решётка: подтверждено соотношение разрешающей способности $\mathcal{R} = mN$ + \end{itemize} + + \item \textbf{Предельные случаи} + \begin{itemize} + \item $F \gg 1$: геометрическая оптика --- резкая тень + \item $F \sim 1$: дифракция Френеля --- осцилляции с неравной огибающей + \item $F \ll 1$: дифракция Фраунгофера --- картина совпадает с преобразованием Фурье апертуры + \end{itemize} +\end{enumerate} + +\subsection{Физический смысл} + +Число Френеля $F = a^2/(\lambda z)$ является единственным управляющим параметром, определяющим режим дифракции. При $F \to \infty$ волновые эффекты исчезают и восстанавливается геометрическая оптика. Принцип Гюйгенса--Френеля позволяет единым образом описать все три режима и является основой современной вычислительной оптики. diff --git a/M11_report/intro.tex b/M11_report/intro.tex new file mode 100644 index 0000000..8089a99 --- /dev/null +++ b/M11_report/intro.tex @@ -0,0 +1,37 @@ +\newpage + +\section{Введение} + +\Task Дифракция --- это явление отклонения волн от прямолинейного распространения при встрече с препятствиями или отверстиями, размеры которых сравнимы с длиной волны. Численное моделирование дифракции позволяет наблюдать переход между различными режимами: геометрической оптикой, дифракцией Фраунгофера и дифракцией Френеля. + +В основе вычислений лежит \textbf{принцип Гюйгенса--Френеля}: каждая точка волнового фронта является источником вторичных сферических волн, а результирующее поле в точке наблюдения есть суперпозиция этих волн с учётом их амплитуд и фаз. + +\Goal +\begin{enumerate} + \item Реализовать численное вычисление интеграла Гюйгенса--Френеля для плоских препятствий + \item Исследовать дифракцию на щели, круглом отверстии, двух щелях и дифракционной решётке + \item Наблюдать предельные случаи: дифракцию Фраунгофера ($F \ll 1$), дифракцию Френеля ($F \sim 1$) и геометрическую оптику ($F \gg 1$) + \item Сравнить численные результаты с аналитическими формулами +\end{enumerate} + +\newpage +\section{Физическая постановка задачи} + +\subsection{Принцип Гюйгенса--Френеля} + +Рассмотрим плоскую монохроматическую волну, падающую нормально на непрозрачный экран с отверстием $\Sigma$ в плоскости $z = 0$. Комплексная амплитуда поля в точке наблюдения $P = (x, y, z)$ при $z > 0$ даётся интегралом Кирхгофа--Гюйгенса: +\[ + U(P) = \frac{1}{i\lambda} \iint_{\Sigma} U_0(\xi, \eta)\, + \frac{e^{ikr}}{r}\,\cos\theta\; d\xi\,d\eta, +\] +где $r = \sqrt{(x-\xi)^2 + (y-\eta)^2 + z^2}$ --- расстояние от точки отверстия $(\xi, \eta)$ до точки наблюдения, $k = 2\pi/\lambda$ --- волновое число, $\cos\theta \approx z/r$ --- фактор наклона. + +\subsection{Число Френеля} + +Характеристический безразмерный параметр задачи: +\[ + F = \frac{a^2}{\lambda z}, +\] +где $a$ --- характерный размер апертуры. При $F \ll 1$ --- режим Фраунгофера, при $F \sim 1$ --- режим Френеля, при $F \gg 1$ --- геометрическая оптика. + +\newpage diff --git a/M11_report/main.tex b/M11_report/main.tex new file mode 100644 index 0000000..dcaa826 --- /dev/null +++ b/M11_report/main.tex @@ -0,0 +1,17 @@ +\input{preamble} + +\begin{document} + +\input{title} +\pagestyle{main} +{ +\centering +\tableofcontents +} +\input{intro} +\input{models/models} +\input{numerical_methods} +\input{results} +\input{conclusion} + +\end{document} diff --git a/M11_report/models/models.tex b/M11_report/models/models.tex new file mode 100644 index 0000000..64222a1 --- /dev/null +++ b/M11_report/models/models.tex @@ -0,0 +1,48 @@ +\section{Модели препятствий} +\subsection{Щель} +Одна щель шириной $2a$, центрированная по оси $\xi = 0$, задаётся функцией пропускания: +\[ + t(\xi, \eta) = \begin{cases} 1, & |\xi| \Le a, \\ 0, & \text{иначе.} \end{cases} +\] +Аналитическое решение в приближении Фраунгофера (для 1D): +\[ + I(x) \propto \left(\frac{\sin u}{u}\right)^2, \qquad u = \frac{\pi a x}{\lambda z}. +\] + +\subsection{Круглое отверстие} + +Круглое отверстие радиуса $R$: +\[ + t(\xi, \eta) = \begin{cases} 1, & \xi^2 + \eta^2 \Le R^2, \\ 0, & \text{иначе.} \end{cases} +\] +В дальней зоне картина обладает осевой симметрией --- кольца Эйри. Интенсивность на оси: +\[ + I(\rho) \propto \left(\frac{2J_1(v)}{v}\right)^2, \qquad v = \frac{k R \rho}{z}, +\] +где $J_1$ --- функция Бесселя первого рода. + +\subsection{Две щели} + +Две щели шириной $2a$, расположенные симметрично с центрами при $\xi = \pm d$: +\[ + t(\xi, \eta) = \begin{cases} 1, & |\xi - d| \Le a \;\text{или}\; |\xi + d| \Le a, \\ 0, & \text{иначе.} \end{cases} +\] +Картина Фраунгофера --- произведение огибающей одиночной щели и интерференционного множителя: +\[ + I(x) \propto \left(\frac{\sin u}{u}\right)^2 \cos^2(v), + \quad u = \frac{\pi a x}{\lambda z},\quad v = \frac{\pi d x}{\lambda z}. +\] + +\subsection{Дифракционная решётка} + +$N$ щелей шириной $2a$ с периодом $p$ ($p > 2a$): +\[ + t(\xi, \eta) = \begin{cases} 1, & \exists\, m \in \{0,\ldots,N-1\}: |\xi - m p| \Le a, \\ 0, & \text{иначе.} \end{cases} +\] +Картина Фраунгофера: +\[ + I(x) \propto \left(\frac{\sin u}{u}\right)^2 + \left(\frac{\sin(N v)}{N \sin v}\right)^2, + \quad u = \frac{\pi a x}{\lambda z},\quad v = \frac{\pi p x}{\lambda z}. +\] +Главные максимумы (порядки решётки) расположены при $v = m\pi$, т.е. при $x_m = m\lambda z / p$. diff --git a/M11_report/numerical_methods.tex b/M11_report/numerical_methods.tex new file mode 100644 index 0000000..a47372e --- /dev/null +++ b/M11_report/numerical_methods.tex @@ -0,0 +1,53 @@ +\newpage +\section{Численные методы} + +\subsection{Прямое численное интегрирование (метод Гюйгенса--Френеля)} + +Интеграл Гюйгенса--Френеля вычисляется методом прямоугольников по дискретной сетке в плоскости апертуры: +\[ + U(x,y) \approx \frac{1}{i\lambda} \sum_{m,n} t(\xi_m, \eta_n)\, + \frac{e^{ikr_{mn}}}{r_{mn}}\cdot\frac{z}{r_{mn}}\;\Delta\xi\,\Delta\eta, +\] +где $r_{mn} = \sqrt{(x-\xi_m)^2 + (y-\eta_n)^2 + z^2}$. + +Сложность метода $O(N_\text{ap}^2 \cdot N_\text{obs}^2)$, что допустимо для умеренных сеток ($N \sim 10^2$--$10^3$). Метод точен во всех режимах (Фраунгофер, Френель, геометрическая оптика). + +\subsection{Быстрое вычисление через БПФ (приближение Френеля)} + +В приближении Френеля интеграл является \textbf{дискретной свёрткой}: +\[ + U(x,y) = \frac{e^{ikz}}{i\lambda z} + \iint t(\xi,\eta)\,e^{\frac{ik}{2z}[(\xi-x)^2+(\eta-y)^2]}\,d\xi\,d\eta, +\] +которая вычисляется за $O(N^2 \log N)$ через быстрое преобразование Фурье: +\[ + U = \mathcal{F}^{-1}\!\left[\mathcal{F}[t \cdot Q_1] \cdot Q_2\right], +\] +где $Q_1$, $Q_2$ --- фазовые множители: +\[ + Q_1(\xi,\eta) = e^{\frac{ik}{2z}(\xi^2+\eta^2)}, \qquad + Q_2(f_x, f_y) = e^{-i\pi\lambda z(f_x^2+f_y^2)}. +\] + +\subsection{Метод углового спектра} + +Для произвольного расстояния $z$ без параксиального приближения используется метод углового спектра: +\[ + \hat{U}(f_x, f_y, z) = \hat{U}_0(f_x, f_y) \cdot + \exp\!\left[iz\sqrt{k^2 - (2\pi f_x)^2 - (2\pi f_y)^2}\right]. +\] +Метод корректен для любого $z$, включая ближнее поле, и реализуется через два БПФ. + +\subsection{Критерий выбора метода} + +\begin{center} +\begin{tabular}{|l|c|c|c|} +\hline +\textbf{Метод} & \textbf{Режим} & \textbf{Сложность} & \textbf{Точность} \\ +\hline +Прямое суммирование & любой & $O(N^4)$ & высокая \\ +БПФ Френеля & $F \lesssim 10^3$ & $O(N^2\log N)$ & параксиальная \\ +Угловой спектр & любой & $O(N^2\log N)$ & высокая \\ +\hline +\end{tabular} +\end{center} diff --git a/M11_report/preamble.tex b/M11_report/preamble.tex new file mode 100644 index 0000000..960ff41 --- /dev/null +++ b/M11_report/preamble.tex @@ -0,0 +1,100 @@ +\documentclass[a4paper,12pt]{article} + +\usepackage{geometry} +% задание полей текста +\geometry{left=10mm,right=10mm,top=20mm,bottom=20mm} +\usepackage{listings} +\usepackage{cmap} % поиск в PDF +\usepackage[T2A]{fontenc} % кодировка +\usepackage[utf8]{inputenc} +\usepackage[english, russian]{babel} % локализация и переносы +\usepackage{natbib} + +\usepackage{graphicx} +\graphicspath{{pix/}} +\usepackage{pgfplots} +\usetikzlibrary{arrows,patterns} +\pgfplotsset{width=10cm,compat=newest} +\pgfkeys{/pgf/trig format=rad} + + +\usepackage{xcolor} +\usepackage{hyperref} +\hypersetup{ + colorlinks=true, + linkcolor=blue, + filecolor=blue, + urlcolor=blue, + citecolor=blue +} +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{bm} +\usepackage{enumerate} +\usepackage[normalem]{ulem} +\usepackage{titlesec} + +% Настройка форматирования section для центрирования без растягивания +\titleformat{\section} + {\normalfont\Large\bfseries\centering}{\thesection}{1em}{} +\titleformat{\subsection} + {\normalfont\large\bfseries\centering}{\thesubsection}{1em}{} +\titleformat{\subsubsection} + {\normalfont\normalsize\bfseries\centering}{\thesubsubsection}{1em}{} + +\setlength\parindent{0pt} + +\sloppy % строго соблюдать границы текста +\linespread{1.2} % коэффициент межстрочного интервала +\setlength{\parskip}{0.5em} % вертик. интервал между абзацами + +\setcounter{secnumdepth}{0} % отключение нумерации разделов +\binoppenalty=1000 % уменьшение переносов в формулах + +% объявление новых макрокоманд + +\newcommand{\Def}{\textbf{Def.} } +\newcommand{\Th}{\textbf{Теорема.} } +\newcommand{\Thbd}{\textbf{Теорема (б/д).} } +\newcommand{\Theor}[1]{\textbf{Теорема ({#1}).} } +\newcommand{\Theorbd}[1]{\textbf{Теорема ({#1}) (б/д).} } +\newcommand{\Consequence}{\textbf{Следствие.} } +\newcommand{\Remind}{\textbf{Remind.} } +\newcommand{\Note}{\textbf{Note.} } +\newcommand{\Statement}{\textbf{Утверждение.} } +\newcommand{\Prop}{\textbf{Свойство:} } +\newcommand{\Props}{\textbf{Свойства:} } +\newcommand{\Proof}{\textbf{Доказательство:} } +\newcommand{\Prooff}{\textbf{Доказать:} } +\newcommand{\Solution}{\textbf{Решение:} } +\newcommand{\Alg}{\textbf{Algorithm.} } +\newcommand{\Lemma}{\textbf{Лемма.} } +\newcommand{\Example}{\textbf{Пример:} } +\newcommand{\Task}{\textbf{Задача.} } +\newcommand{\Goal}{\textbf{Цель работы:} } +\newcommand{\Solve}{\textbf{Решение:} } +\newcommand{\Examples}{\textbf{Примеры.} } + +\allowdisplaybreaks[4] + +\newcommand{\Endproof}{$\blacksquare$ } + +\newcommand{\tr}{\text{tr}} +\newcommand{\Le}{\leqslant} +\newcommand{\Ge}{\geqslant} +\newcommand{\A}{\mathcal{A}} +\newcommand{\M}{\mathcal{M}} +\newcommand{\F}{\mathcal{F}} +\newcommand{\Gs}{\mathcal{G}} +\newcommand{\R}{\mathbb{R}} +\newcommand{\N}{\mathbb{N}} +\newcommand{\Q}{\mathbb{Q}} +\newcommand{\Norm}{\mathcal{N}} +\newcommand{\ind}{\perp \!\!\! \perp} + +\newpagestyle{main}{ + \setheadrule{.4pt} + \sethead{}{\bullet \; \textit{Цифровизация физических процессов} \; \bullet}{} + \setfootrule{.4pt} + \setfoot{ВШПИ МФТИ}{М11. Дифракция}{Б13-402} +} diff --git a/M11_report/results.tex b/M11_report/results.tex new file mode 100644 index 0000000..eb2add1 --- /dev/null +++ b/M11_report/results.tex @@ -0,0 +1,42 @@ +\newpage +\section{Результаты} + +\subsection{Щель: переход между режимами} + +Для щели шириной $a = 0.5$~мм при $\lambda = 633$~нм исследовалась зависимость картины от расстояния $z$: + +\begin{itemize} + \item \textbf{$F \gg 1$ (геометрическая оптика)}: при $z = 1$~мм ($F \approx 400$) наблюдается прямоугольная тень --- интенсивность практически равна $I_0$ внутри и нулю снаружи геометрической тени. + \item \textbf{$F \sim 1$ (дифракция Френеля)}: при $z = 40$~см ($F \approx 1$) появляются характерные осцилляции Френеля с неравными максимумами и интенсивностью на оси, превышающей падающую. + \item \textbf{$F \ll 1$ (дифракция Фраунгофера)}: при $z = 10$~м ($F \approx 0.04$) картина воспроизводит функцию $\sinc^2$: центральный максимум с боковыми долями, нули при $x_m = m\lambda z/a$. +\end{itemize} + +\subsection{Круглое отверстие: зоны Френеля} + +Для отверстия радиуса $R$ интенсивность на оси определяется числом открытых зон Френеля: +\[ + N_F = \frac{R^2}{\lambda z}, \qquad + I_{\text{ось}} = \begin{cases} 4I_0, & N_F \text{ нечётное}, \\ 0, & N_F \text{ чётное.} \end{cases} +\] +Численный расчёт воспроизводит эти колебания при изменении $z$. В дальней зоне формируется кольцевая картина Эйри с угловым радиусом первого нуля $\theta = 1.22\,\lambda / (2R)$. + +\subsection{Двойная щель: опыт Юнга} + +Для двух щелей с $a = 0.1$~мм, $d = 0.5$~мм получена интерференционная картина с периодом: +\[ + \Delta x = \frac{\lambda z}{2d}. +\] +Огибающая определяется дифракцией на одной щели. При увеличении $z$ картина растягивается, при уменьшении переходит в режим Френеля с искажением интерференционных полос. + +\subsection{Дифракционная решётка} + +Для решётки из $N = 10$ щелей ($a = 0.05$~мм, $p = 0.2$~мм): +\begin{itemize} + \item Главные максимумы расположены при $\sin\theta_m = m\lambda/p$ + \item Угловая дисперсия $D = m / (p\cos\theta_m)$ + \item Разрешающая способность $\mathcal{R} = mN$ + \item Между соседними главными максимумами наблюдается $N-2 = 8$ побочных максимумов +\end{itemize} + +При увеличении $N$ главные максимумы сужаются обратно пропорционально $N$, что воспроизводится численно. + diff --git a/M11_report/title.tex b/M11_report/title.tex new file mode 100644 index 0000000..0a93c43 --- /dev/null +++ b/M11_report/title.tex @@ -0,0 +1,17 @@ +\begin{titlepage} + \centering + {\scshape\Large Московский физико-технический институт \par} + \vspace{1cm} + {\scshape\Large Высшая школа программной инженерии \par} + \vspace{2cm} + {\huge \textbf{М11. Дифракция} \par} + \vspace{1.5cm} + {\Large \textit{Численное моделирование дифракции плоской монохроматической волны на плоских препятствиях по принципу Гюйгенса--Френеля} \par} + \vspace{2cm} + {\Large Выполнили студенты Б13-402: \par} + {\Large Жердев Егор \par} + {\Large Савельев Данил \par} + \vfill + {\large Долгопрудный \par} + {\large \today \par} +\end{titlepage} diff --git a/M5_tunnel_report/M5_tunnel_report.pdf b/M5_tunnel_report/M5_tunnel_report.pdf new file mode 100644 index 0000000..5e50276 Binary files /dev/null and b/M5_tunnel_report/M5_tunnel_report.pdf differ diff --git a/M5_tunnel_report/conclusion.tex b/M5_tunnel_report/conclusion.tex new file mode 100644 index 0000000..8894776 --- /dev/null +++ b/M5_tunnel_report/conclusion.tex @@ -0,0 +1,30 @@ +\section{Заключение} + +В работе смоделирована работа автогенератора колебаний на туннельном диоде. + +\textbf{Основные результаты:} + +\begin{enumerate} + \item \textbf{Модель цепи} + \begin{itemize} + \item Составлена система двух ОДУ из законов Кирхгофа + \item ВАХ диода аппроксимирована кубической параболой с участком отрицательного сопротивления + \end{itemize} + + \item \textbf{Условие автогенерации} + \begin{itemize} + \item Автогенерация возникает, когда $R < R_{\text{кр}} = -I_D'(U^*)\cdot L/C$ + \item Рабочая точка должна лежать на участке с отрицательным дифференциальным сопротивлением + \end{itemize} + + \item \textbf{Режимы колебаний} + \begin{itemize} + \item При высокой добротности $Q$ и малой надкритичности --- почти чистая синусоида + \item При малом $R$ --- релаксационные (нелинейные) колебания большой амплитуды + \item Максимальная амплитуда ограничена шириной рабочего участка ВАХ + \end{itemize} +\end{enumerate} + +\subsection{Физический смысл} + +Туннельный диод компенсирует потери энергии в резисторе за счёт участка с отрицательным сопротивлением --- контур «подкачивается» энергией от источника $\mathcal{E}$ через диод. При балансе потерь и генерации устанавливаются незатухающие колебания --- предельный цикл на фазовой плоскости. diff --git a/M5_tunnel_report/intro.tex b/M5_tunnel_report/intro.tex new file mode 100644 index 0000000..9a9f04a --- /dev/null +++ b/M5_tunnel_report/intro.tex @@ -0,0 +1,59 @@ +\newpage + +\section{Введение} + +\Task Туннельный диод --- это элемент с нелинейной вольт-амперной характеристикой (ВАХ), которая имеет участок с \textbf{отрицательным дифференциальным сопротивлением} (ток убывает с ростом напряжения). Этот участок позволяет построить автогенератор --- устройство, которое само поддерживает незатухающие колебания без внешнего источника переменного сигнала. + +Схема генератора содержит катушку индуктивности $L$, конденсатор $C$, резистор $R$ и туннельный диод $D$, питающиеся от источника постоянного напряжения $\mathcal{E}$. + +ВАХ диода при $U_D > 0$ аппроксимируется кубической параболой: +\[ + I_D(U_D) = aU_D^3 + bU_D^2 + cU_D, \qquad + a = 4\;\frac{\text{мА}}{\text{В}^3},\quad + b = -16\;\frac{\text{мА}}{\text{В}^2},\quad + c = 17\;\frac{\text{мА}}{\text{В}} +\] +При $U_D < 0$: $I_D = 0$. + +\Goal +\begin{enumerate} + \item Составить уравнения, описывающие работу схемы + \item Реализовать численное решение этих уравнений + \item Исследовать режимы работы в зависимости от параметров $R$, $C$, $L$ + \item Найти условия автогенерации и получения чистой синусоиды + \item Определить максимальную амплитуду колебаний +\end{enumerate} + +\newpage +\section{Физическая постановка задачи} + +Рассматривается схема: источник $\mathcal{E}$, резистор $R$, параллельный колебательный контур $LC$ и туннельный диод $D$ --- все соединены последовательно. + +\textbf{Обозначения:} +\begin{itemize} + \item $U_C$ --- напряжение на конденсаторе + \item $I_L$ --- ток через катушку + \item $U_D = U_C$ --- напряжение на диоде (диод параллелен контуру) +\end{itemize} + +\textbf{Уравнения цепи} (из законов Кирхгофа): + +Для тока через конденсатор: +\[ + C\,\frac{dU_C}{dt} = I_L - I_D(U_C) +\] + +Для тока через катушку: +\[ + L\,\frac{dI_L}{dt} = \mathcal{E} - R\,I_L - U_C +\] + +Это система двух обыкновенных дифференциальных уравнений первого порядка. Её и нужно решать численно. + +\textbf{Параметры системы:} +\begin{itemize} + \item $\mathcal{E}$ --- напряжение питания (выбирается в рабочей точке на участке с отрицательным сопротивлением) + \item $R$, $C$, $L$ --- параметры, от которых зависит режим работы +\end{itemize} + +\newpage diff --git a/M5_tunnel_report/main.tex b/M5_tunnel_report/main.tex new file mode 100644 index 0000000..dcaa826 --- /dev/null +++ b/M5_tunnel_report/main.tex @@ -0,0 +1,17 @@ +\input{preamble} + +\begin{document} + +\input{title} +\pagestyle{main} +{ +\centering +\tableofcontents +} +\input{intro} +\input{models/models} +\input{numerical_methods} +\input{results} +\input{conclusion} + +\end{document} diff --git a/M5_tunnel_report/models/autogeneration.tex b/M5_tunnel_report/models/autogeneration.tex new file mode 100644 index 0000000..1fff878 --- /dev/null +++ b/M5_tunnel_report/models/autogeneration.tex @@ -0,0 +1,28 @@ +\subsection{Режимы работы генератора} + +\subsubsection{Условие автогенерации} + +Автогенерация возникает, когда отрицательное сопротивление диода компенсирует потери в резисторе. Из условия $\gamma_{\text{eff}} < 0$: +\[ +\frac{R}{L} + \frac{I_D'(U^*)}{C} < 0 +\quad \Longleftrightarrow \quad +I_D'(U^*) < -\frac{RC}{L} +\] + +Поскольку $I_D'(U^*) < 0$ на рабочем участке, это выполняется при достаточно малом $R$ и достаточно большом $L/C$. + +\subsubsection{Чистая синусоида} + +При малой амплитуде колебания близки к синусоидальным. Это возможно, если: +\begin{itemize} + \item Рабочая точка лежит вблизи минимума $|r_D|$ (наиболее «мягкое» возбуждение) + \item Добротность контура $Q = \dfrac{1}{R}\sqrt{\dfrac{L}{C}}$ достаточно велика --- контур хорошо фильтрует гармоники +\end{itemize} + +При большой добротности нелинейные члены ВАХ слабо влияют на форму сигнала --- генератор работает как почти линейный. + +\subsubsection{Максимальная амплитуда} + +Амплитуда колебаний ограничена нелинейностью ВАХ: при нарастании колебаний отрицательное сопротивление уменьшается по модулю, пока не установится баланс между генерацией и потерями. Максимальная амплитуда определяется из условия, что средняя мощность, отдаваемая диодом, равна мощности, рассеиваемой в $R$. + +Численно это соответствует установившемуся циклу на фазовой плоскости $(U_C, I_L)$ --- \textbf{предельному циклу}. diff --git a/M5_tunnel_report/models/models.tex b/M5_tunnel_report/models/models.tex new file mode 100644 index 0000000..9939bdb --- /dev/null +++ b/M5_tunnel_report/models/models.tex @@ -0,0 +1,8 @@ +\section{Физические модели} + +\input{models/vac} +\newpage +\input{models/oscillations} +\newpage +\input{models/autogeneration} +\newpage diff --git a/M5_tunnel_report/models/oscillations.tex b/M5_tunnel_report/models/oscillations.tex new file mode 100644 index 0000000..b1d22d6 --- /dev/null +++ b/M5_tunnel_report/models/oscillations.tex @@ -0,0 +1,44 @@ +\subsection{Колебания в контуре} + +\subsubsection{Система уравнений} + +Запишем систему в виде, удобном для численного решения. Введём вектор состояния $(U_C, I_L)$: + +\[ +\begin{cases} +\dfrac{dU_C}{dt} = \dfrac{1}{C}\bigl[I_L - I_D(U_C)\bigr] \\[1em] +\dfrac{dI_L}{dt} = \dfrac{1}{L}\bigl[\mathcal{E} - R\,I_L - U_C\bigr] +\end{cases} +\] + +Это нелинейная система из-за члена $I_D(U_C)$. + +\subsubsection{Линеаризация вблизи рабочей точки} + +Вблизи рабочей точки $(U^*, I^*)$ введём малые отклонения $u = U_C - U^*$, $\,i = I_L - I^*$. Тогда: +\[ +\begin{cases} +C\,\dot{u} = i - I_D'(U^*)\cdot u \\[0.5em] +L\,\dot{i} = -R\,i - u +\end{cases} +\] + +Это линейная система. Её решение --- затухающие или нарастающие колебания в зависимости от знака \textbf{эффективного затухания}: +\[ +\gamma_{\text{eff}} = \frac{R}{2L} + \frac{I_D'(U^*)}{2C} +\] + +\begin{itemize} + \item $\gamma_{\text{eff}} > 0$ --- колебания затухают, генерации нет + \item $\gamma_{\text{eff}} < 0$ --- колебания нарастают, возникает автогенерация + \item $\gamma_{\text{eff}} = 0$ --- граница генерации +\end{itemize} + +\subsubsection{Собственная частота контура} + +В линейном приближении контур колеблется на частоте: +\[ +\omega_0 = \frac{1}{\sqrt{LC}} +\] + +При нелинейном режиме (большие амплитуды) форма колебаний отличается от синусоиды. diff --git a/M5_tunnel_report/models/vac.tex b/M5_tunnel_report/models/vac.tex new file mode 100644 index 0000000..4de04e0 --- /dev/null +++ b/M5_tunnel_report/models/vac.tex @@ -0,0 +1,27 @@ +\subsection{Вольт-амперная характеристика туннельного диода} + +ВАХ туннельного диода имеет нестандартную форму: после начального роста тока он \textit{убывает} при увеличении напряжения --- это участок \textbf{отрицательного дифференциального сопротивления}. + +Аппроксимация кубической параболой при $U_D > 0$: +\[ +I_D(U_D) = aU_D^3 + bU_D^2 + cU_D +\] +\[ +a = 4\;\frac{\text{мА}}{\text{В}^3},\quad b = -16\;\frac{\text{мА}}{\text{В}^2},\quad c = 17\;\frac{\text{мА}}{\text{В}} +\] + +При $U_D \Le 0$: $\;I_D = 0$. + +\textbf{Дифференциальное сопротивление диода} --- это производная ВАХ: +\[ +r_D(U_D) = \frac{1}{I_D'(U_D)} = \frac{1}{3aU_D^2 + 2bU_D + c} +\] + +На участке, где $I_D'(U_D) < 0$, дифференциальное сопротивление отрицательно. Именно этот участок обеспечивает генерацию. + +\textbf{Рабочая точка} --- равновесие системы, при котором напряжение $U_C = U^*$ и ток $I_L = I^*$ не меняются. Она определяется пересечением нагрузочной прямой с ВАХ: +\[ +I^* = \frac{\mathcal{E} - U^*}{R}, \qquad I^* = I_D(U^*) +\] + +Для автогенерации рабочую точку выбирают на участке с отрицательным сопротивлением. diff --git a/M5_tunnel_report/numerical_methods.tex b/M5_tunnel_report/numerical_methods.tex new file mode 100644 index 0000000..54157e5 --- /dev/null +++ b/M5_tunnel_report/numerical_methods.tex @@ -0,0 +1,43 @@ +\section{Численные методы решения} + +\subsection{Метод Рунге--Кутты 4-го порядка (РК4)} + +Система двух ОДУ решается методом РК4. Обозначим вектор состояния: +\[ +\vec{u} = (U_C,\; I_L)^\top, \qquad \dot{\vec{u}} = \vec{f}(\vec{u}) +\] + +Один шаг метода: +\begin{align*} + \vec{k}_1 &= h\,\vec{f}(\vec{u}_n) \\ + \vec{k}_2 &= h\,\vec{f}(\vec{u}_n + \vec{k}_1/2) \\ + \vec{k}_3 &= h\,\vec{f}(\vec{u}_n + \vec{k}_2/2) \\ + \vec{k}_4 &= h\,\vec{f}(\vec{u}_n + \vec{k}_3) \\[0.5em] + \vec{u}_{n+1} &= \vec{u}_n + \frac{1}{6}(\vec{k}_1 + 2\vec{k}_2 + 2\vec{k}_3 + \vec{k}_4) +\end{align*} + +\subsection{Выбор шага интегрирования} + +Шаг выбирается значительно меньше периода колебаний: +\[ +h \ll T_0 = 2\pi\sqrt{LC} +\] +Обычно достаточно $h = T_0 / 1000$. + +\subsection{Нахождение рабочей точки} + +Рабочая точка $(U^*, I^*)$ находится численно как пересечение нагрузочной прямой с ВАХ: +\[ +f(U) = I_D(U) - \frac{\mathcal{E} - U}{R} = 0 +\] +решается методом бисекции или \texttt{scipy.optimize.brentq} на интервале $[0,\, \mathcal{E}]$. + +\subsection{Анализ результатов} + +После численного решения проводится: +\begin{enumerate} + \item \textbf{Фазовый портрет} --- график $(U_C(t),\, I_L(t))$. Предельный цикл указывает на установившиеся колебания. + \item \textbf{Временна́я зависимость} $U_C(t)$ --- оценка формы (синусоида или нет) и амплитуды. + \item \textbf{Спектр} --- быстрое преобразование Фурье (\texttt{numpy.fft}) для оценки чистоты синусоиды: чем меньше высшие гармоники, тем чище сигнал. + \item \textbf{Параметрическое исследование} --- сканирование по $R$, $C$, $L$ для построения карты режимов. +\end{enumerate} diff --git a/M5_tunnel_report/preamble.tex b/M5_tunnel_report/preamble.tex new file mode 100644 index 0000000..75947ca --- /dev/null +++ b/M5_tunnel_report/preamble.tex @@ -0,0 +1,100 @@ +\documentclass[a4paper,12pt]{article} + +\usepackage{geometry} +% задание полей текста +\geometry{left=10mm,right=10mm,top=20mm,bottom=20mm} +\usepackage{listings} +\usepackage{cmap} % поиск в PDF +\usepackage[T2A]{fontenc} % кодировка +\usepackage[utf8]{inputenc} +\usepackage[english, russian]{babel} % локализация и переносы +\usepackage{natbib} + +\usepackage{graphicx} +\graphicspath{{pix/}} +\usepackage{pgfplots} +\usetikzlibrary{arrows,patterns} +\pgfplotsset{width=10cm,compat=newest} +\pgfkeys{/pgf/trig format=rad} + + +\usepackage{xcolor} +\usepackage{hyperref} +\hypersetup{ + colorlinks=true, + linkcolor=blue, + filecolor=blue, + urlcolor=blue, + citecolor=blue +} +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{bm} +\usepackage{enumerate} +\usepackage[normalem]{ulem} +\usepackage{titlesec} + +% Настройка форматирования section для центрирования без растягивания +\titleformat{\section} + {\normalfont\Large\bfseries\centering}{\thesection}{1em}{} +\titleformat{\subsection} + {\normalfont\large\bfseries\centering}{\thesubsection}{1em}{} +\titleformat{\subsubsection} + {\normalfont\normalsize\bfseries\centering}{\thesubsubsection}{1em}{} + +\setlength\parindent{0pt} + +\sloppy % строго соблюдать границы текста +\linespread{1.2} % коэффициент межстрочного интервала +\setlength{\parskip}{0.5em} % вертик. интервал между абзацами + +\setcounter{secnumdepth}{0} % отключение нумерации разделов +\binoppenalty=1000 % уменьшение переносов в формулах + +% объявление новых макрокоманд + +\newcommand{\Def}{\textbf{Def.} } +\newcommand{\Th}{\textbf{Теорема.} } +\newcommand{\Thbd}{\textbf{Теорема (б/д).} } +\newcommand{\Theor}[1]{\textbf{Теорема ({#1}).} } +\newcommand{\Theorbd}[1]{\textbf{Теорема ({#1}) (б/д).} } +\newcommand{\Consequence}{\textbf{Следствие.} } +\newcommand{\Remind}{\textbf{Remind.} } +\newcommand{\Note}{\textbf{Note.} } +\newcommand{\Statement}{\textbf{Утверждение.} } +\newcommand{\Prop}{\textbf{Свойство:} } +\newcommand{\Props}{\textbf{Свойства:} } +\newcommand{\Proof}{\textbf{Доказательство:} } +\newcommand{\Prooff}{\textbf{Доказать:} } +\newcommand{\Solution}{\textbf{Решение:} } +\newcommand{\Alg}{\textbf{Algorithm.} } +\newcommand{\Lemma}{\textbf{Лемма.} } +\newcommand{\Example}{\textbf{Пример:} } +\newcommand{\Task}{\textbf{Задача.} } +\newcommand{\Goal}{\textbf{Цель работы:} } +\newcommand{\Solve}{\textbf{Решение:} } +\newcommand{\Examples}{\textbf{Примеры.} } + +\allowdisplaybreaks[4] + +\newcommand{\Endproof}{$\blacksquare$ } + +\newcommand{\tr}{\text{tr}} +\newcommand{\Le}{\leqslant} +\newcommand{\Ge}{\geqslant} +\newcommand{\A}{\mathcal{A}} +\newcommand{\M}{\mathcal{M}} +\newcommand{\F}{\mathcal{F}} +\newcommand{\Gs}{\mathcal{G}} +\newcommand{\R}{\mathbb{R}} +\newcommand{\N}{\mathbb{N}} +\newcommand{\Q}{\mathbb{Q}} +\newcommand{\Norm}{\mathcal{N}} +\newcommand{\ind}{\perp \!\!\! \perp} + +\newpagestyle{main}{ + \setheadrule{.4pt} + \sethead{}{\bullet \; \textit{Цифровизация физических процессов} \; \bullet}{} + \setfootrule{.4pt} + \setfoot{ВШПИ МФТИ}{М5. Автогенератор на туннельном диоде}{Б13-402} +} diff --git a/M5_tunnel_report/results.tex b/M5_tunnel_report/results.tex new file mode 100644 index 0000000..560a6b9 --- /dev/null +++ b/M5_tunnel_report/results.tex @@ -0,0 +1,43 @@ +\section{Результаты} + +\subsection{ВАХ туннельного диода} + +Кубическая аппроксимация воспроизводит характерную форму ВАХ: +\begin{itemize} + \item Участок роста при малых $U_D$ + \item Пик тока вблизи $U_D \approx 1\;\text{В}$ + \item Участок отрицательного дифференциального сопротивления + \item Повторный рост при больших $U_D$ +\end{itemize} + +\subsection{Условие автогенерации} + +Автогенерация возникает при: +\[ +R < R_{\text{кр}} = -I_D'(U^*) \cdot \frac{L}{C} +\] + +При фиксированных $L$ и $C$ с уменьшением $R$ система переходит от затухающих колебаний к нарастающим, а затем к установившимся. + +\subsection{Чистая синусоида} + +Наиболее близкие к синусоидальным колебания наблюдаются при: +\begin{itemize} + \item Большой добротности контура $Q \gg 1$ + \item Малой надкритичности ($R$ чуть меньше $R_{\text{кр}}$) +\end{itemize} +В этом режиме высшие гармоники в спектре пренебрежимо малы. + +\subsection{Максимальная амплитуда} + +При глубоком вхождении в режим генерации (малое $R$) амплитуда растёт, форма сигнала становится нелинейной (релаксационные колебания). Максимальная амплитуда ограничена шириной участка с отрицательным сопротивлением на ВАХ и не превышает примерно половины этого диапазона напряжений. + +\subsection{Карта режимов} + +На плоскости параметров $(R, L/C)$: +\begin{itemize} + \item Область выше кривой $R_{\text{кр}}(L/C)$ --- затухающие колебания + \item Область ниже --- автогенерация + \item При малом $R$ и малом $L/C$ --- релаксационные колебания (не синусоида) + \item При малом $R$ и большом $L/C$ --- близкие к синусоидальным колебания большой амплитуды +\end{itemize} diff --git a/M5_tunnel_report/title.tex b/M5_tunnel_report/title.tex new file mode 100644 index 0000000..85dc95d --- /dev/null +++ b/M5_tunnel_report/title.tex @@ -0,0 +1,18 @@ +\begin{titlepage} + \centering + {\scshape\Large Московский физико-технический институт \par} + \vspace{1cm} + {\scshape\Large Высшая школа программной инженерии \par} + \vspace{2cm} + {\huge \textbf{М5. Автогенератор на туннельном диоде} \par} + \vspace{1.5cm} + {\Large \textit{Численное моделирование колебаний в цепи \\ + с нелинейным элементом} \par} + \vspace{2cm} + {\Large Выполнили студенты Б13-402: \par} + {\Large Жердев Егор \par} + {\Large Савельев Данил \par} + \vfill + {\large Долгопрудный \par} + {\large \today \par} +\end{titlepage} diff --git a/src/main.py b/src/main.py index d0cfb45..c3ba00a 100644 --- a/src/main.py +++ b/src/main.py @@ -1,6 +1,7 @@ import streamlit from models.billiard import page as billiard_page +from models.diffraction import page as diffraction_page from models.ideal_gas import page as ideal_gas_page from models.pendulum import page as pendulum_page from models.rolling_ball import page as rolling_ball_page @@ -24,6 +25,11 @@ def main() -> None: title="M7Б Статистика идеального газа со столкновениями частиц", url_path="ideal_gas", ), + streamlit.Page( + diffraction_page, + title="M11 Дифракция", + url_path="diffraction", + ), ]) pages.run() diff --git a/src/models/diffraction/__init__.py b/src/models/diffraction/__init__.py new file mode 100644 index 0000000..36dd42b --- /dev/null +++ b/src/models/diffraction/__init__.py @@ -0,0 +1,110 @@ +from typing import Final + +import streamlit as st + +from models.diffraction.charts import ( + create_intensity_1d_chart, + create_intensity_2d_chart, +) +from models.diffraction.config import ApertureType, DiffractionConfig +from models.diffraction.objects import DiffractionSimulation +from models.diffraction.sidebar import create_diffraction_sidebar +from models.diffraction.utils import analytical_intensity, fresnel_number + +_FRAUNHOFER_THRESHOLD: Final[float] = 0.1 +_GEOMETRIC_THRESHOLD: Final[float] = 10.0 + + +def page() -> None: + st.set_page_config(page_title="М11. Дифракция", layout="wide") + st.title("М11. Дифракция") + st.write( + "Численное моделирование дифракции Фраунгофера и Френеля" + " методом БПФ с сравнением с аналитическими формулами." + ) + + config = create_diffraction_sidebar() + sim = DiffractionSimulation() + + try: + _display_page(sim, config) + except Exception as e: + st.error(f"Ошибка при моделировании: {e}") + + +def _display_page( # noqa: PLR0914 + sim: DiffractionSimulation, + config: DiffractionConfig, +) -> None: + result = sim.simulate(config) + + wavelength = config.wavelength_nm * 1e-9 + a = config.aperture_size_mm * 1e-3 + f_num = fresnel_number(a, wavelength, config.distance_m) + + col_info1, col_info2, col_info3 = st.columns(3) + with col_info1: + st.metric("Число Френеля F", f"{f_num:.3f}") + with col_info2: + regime = ( + "Фраунгофер" + if f_num < _FRAUNHOFER_THRESHOLD + else "Френель" + if f_num < _GEOMETRIC_THRESHOLD + else "Геом. оптика" + ) + st.metric("Режим", regime) + with col_info3: + x_first_min_mm = wavelength * config.distance_m / a * 1e3 + st.metric("Первый минимум (мм)", f"{x_first_min_mm:.3f}") + + x_mm = result.x_obs_m * 1e3 + analytical = analytical_intensity( + result.x_obs_m, config, config.aperture_type + ) + + col1, col2 = st.columns([2, 1]) + with col1: + st.subheader("1D картина дифракции") + chart_1d = create_intensity_1d_chart( + x_mm, + result.intensity_1d, + analytical, + config.aperture_type.value, + ) + st.altair_chart(chart_1d, use_container_width=True) + + with col2: + st.subheader("2D картина дифракции") + chart_2d = create_intensity_2d_chart(result.intensity_2d, x_mm) + st.altair_chart(chart_2d, use_container_width=True) + + st.subheader("Параметры апертуры") + col3, col4 = st.columns(2) + with col3: + st.write(f"• Длина волны λ = {config.wavelength_nm:.0f} нм") + st.write(f"• Расстояние z = {config.distance_m:.2f} м") + st.write(f"• Апертура a = {config.aperture_size_mm:.3f} мм") + with col4: + if config.aperture_type == ApertureType.DOUBLE_SLIT: + d_mm = config.separation_mm + st.write(f"• Расстояние между щелями d = {d_mm:.2f} мм") + fringe_mm = ( + wavelength + * config.distance_m + / (config.separation_mm * 1e-3) + * 1e3 + ) + st.write(f"• Период полос Δx = {fringe_mm:.3f} мм") + if config.aperture_type == ApertureType.GRATING: + st.write(f"• Период решётки p = {config.period_mm:.3f} мм") + st.write(f"• Число щелей N = {config.num_slits}") + order1_mm = ( + wavelength * config.distance_m / (config.period_mm * 1e-3) * 1e3 + ) + st.write(f"• Положение 1-го максимума = {order1_mm:.3f} мм") + if config.aperture_type == ApertureType.CIRCULAR: + airy_radius_mm = ( + 1.22 * wavelength * config.distance_m / (2 * a) * 1e3 + ) + st.write(f"• Радиус диска Эйри = {airy_radius_mm:.3f} мм") diff --git a/src/models/diffraction/charts.py b/src/models/diffraction/charts.py new file mode 100644 index 0000000..ec63dd2 --- /dev/null +++ b/src/models/diffraction/charts.py @@ -0,0 +1,83 @@ +import altair as alt +import numpy as np +import pandas as pd + + +def create_intensity_1d_chart( + x_mm: np.ndarray, + intensity_numerical: np.ndarray, + intensity_analytical: np.ndarray, + aperture_label: str, +) -> alt.Chart: + i_max = intensity_numerical.max() + norm = i_max if i_max > 0 else 1.0 + + df = pd.DataFrame({ + "x (мм)": np.concatenate([x_mm, x_mm]), + "Интенсивность": np.concatenate([ + intensity_numerical / norm, + intensity_analytical, + ]), + "Метод": ( + ["Численный (ФПФ)"] * len(x_mm) + + [f"Аналитика ({aperture_label})"] * len(x_mm) + ), + }) + + return ( + alt.Chart(df) + .mark_line(strokeWidth=2) + .add_params(alt.selection_interval(bind="scales")) + .encode( + x=alt.X("x (мм):Q", title="x (мм)"), + y=alt.Y("Интенсивность:Q", title="I / I₀"), + color=alt.Color("Метод:N"), + tooltip=["x (мм):Q", "Интенсивность:Q", "Метод:N"], + ) + .properties( + title="Картина дифракции (1D срез)", + width=600, + height=350, + ) + ) + + +def create_intensity_2d_chart( + intensity_2d: np.ndarray, + x_obs_mm: np.ndarray, +) -> alt.Chart: + n_display = 64 + step = max(1, len(x_obs_mm) // n_display) + i_disp = intensity_2d[::step, ::step] + x_disp = x_obs_mm[::step] + + i_max = i_disp.max() + norm = i_max if i_max > 0 else 1.0 + + rows, cols = np.meshgrid( + range(len(x_disp)), range(len(x_disp)), indexing="ij" + ) + df = pd.DataFrame({ + "x (мм)": x_disp[cols.ravel()], + "y (мм)": x_disp[rows.ravel()], + "I": i_disp.ravel() / norm, + }) + + return ( + alt.Chart(df) + .mark_rect() + .encode( + x=alt.X("x (мм):Q", bin=alt.Bin(maxbins=64), title="x (мм)"), + y=alt.Y("y (мм):Q", bin=alt.Bin(maxbins=64), title="y (мм)"), + color=alt.Color( + "mean(I):Q", + scale=alt.Scale(scheme="viridis"), + title="I / I₀", + ), + ) + .properties( + title="Картина дифракции 2D", + width=350, + height=350, + ) + ) diff --git a/src/models/diffraction/config.py b/src/models/diffraction/config.py new file mode 100644 index 0000000..337018f --- /dev/null +++ b/src/models/diffraction/config.py @@ -0,0 +1,30 @@ +from dataclasses import dataclass +from enum import Enum +from typing import Final + + +class ApertureType(Enum): + SINGLE_SLIT = "Одиночная щель" + CIRCULAR = "Круглое отверстие" + DOUBLE_SLIT = "Двойная щель" + GRATING = "Дифракционная решётка" + + +DEFAULT_WAVELENGTH_NM: Final[float] = 633.0 +DEFAULT_APERTURE_MM: Final[float] = 0.5 +DEFAULT_DISTANCE_M: Final[float] = 1.0 +DEFAULT_GRID_SIZE: Final[int] = 512 +DEFAULT_APERTURE_GRID_MM: Final[float] = 5.0 + + +@dataclass +class DiffractionConfig: + aperture_type: ApertureType = ApertureType.SINGLE_SLIT + wavelength_nm: float = DEFAULT_WAVELENGTH_NM + aperture_size_mm: float = DEFAULT_APERTURE_MM + distance_m: float = DEFAULT_DISTANCE_M + separation_mm: float = 1.0 + num_slits: int = 5 + period_mm: float = 0.2 + grid_size: int = DEFAULT_GRID_SIZE + aperture_grid_mm: float = DEFAULT_APERTURE_GRID_MM diff --git a/src/models/diffraction/objects.py b/src/models/diffraction/objects.py new file mode 100644 index 0000000..74a1630 --- /dev/null +++ b/src/models/diffraction/objects.py @@ -0,0 +1,77 @@ +from dataclasses import dataclass + +import numpy as np +import streamlit as st + +from models.diffraction.config import ApertureType, DiffractionConfig + + +@dataclass +class DiffractionResult: + intensity_2d: np.ndarray + intensity_1d: np.ndarray + x_obs_m: np.ndarray + + +class DiffractionSimulation: + @staticmethod + def create_aperture( + xi: np.ndarray, + config: DiffractionConfig, + aperture_type: ApertureType, + ) -> np.ndarray: + a = config.aperture_size_mm * 1e-3 + xi_grid, eta_grid = np.meshgrid(xi, xi) + + if aperture_type == ApertureType.SINGLE_SLIT: + return (np.abs(xi_grid) <= a).astype(float) + + if aperture_type == ApertureType.CIRCULAR: + return (xi_grid**2 + eta_grid**2 <= a**2).astype(float) + + if aperture_type == ApertureType.DOUBLE_SLIT: + d = config.separation_mm * 1e-3 + return ( + (np.abs(xi_grid - d / 2) <= a) | (np.abs(xi_grid + d / 2) <= a) + ).astype(float) + + # GRATING + p = config.period_mm * 1e-3 + mask = np.zeros_like(xi_grid, dtype=float) + for m in range(config.num_slits): + center = (m - (config.num_slits - 1) / 2) * p + mask[np.abs(xi_grid - center) <= a] = 1.0 + return mask + + @st.cache_data(ttl=300) + def simulate( # noqa: PLR0914 + _self, # noqa: N805 + config: DiffractionConfig, + ) -> DiffractionResult: + n = config.grid_size + wavelength = config.wavelength_nm * 1e-9 + z = config.distance_m + k = 2 * np.pi / wavelength + + l_ap = config.aperture_grid_mm * 1e-3 + dx = l_ap / n + xi = np.linspace(-l_ap / 2, l_ap / 2 - dx, n) + xi_grid, eta_grid = np.meshgrid(xi, xi) + + aperture = _self.create_aperture(xi, config, config.aperture_type) + q1 = np.exp(1j * k / (2 * z) * (xi_grid**2 + eta_grid**2)) + + fx = np.fft.fftfreq(n, d=dx) + + spectrum = np.fft.fft2(aperture * q1) + intensity_2d = np.fft.fftshift(np.abs(spectrum) ** 2) + x_obs = np.fft.fftshift(fx) * wavelength * z + + center = n // 2 + intensity_1d = intensity_2d[center] + + return DiffractionResult( + intensity_2d=intensity_2d, + intensity_1d=intensity_1d, + x_obs_m=x_obs, + ) diff --git a/src/models/diffraction/sidebar.py b/src/models/diffraction/sidebar.py new file mode 100644 index 0000000..e7c28db --- /dev/null +++ b/src/models/diffraction/sidebar.py @@ -0,0 +1,83 @@ +import streamlit as st + +from models.diffraction.config import ApertureType, DiffractionConfig + + +def create_diffraction_sidebar() -> DiffractionConfig: + st.sidebar.header("Параметры волны") + + wavelength_nm = st.sidebar.number_input( + "Длина волны λ (нм)", + min_value=200.0, + max_value=1000.0, + value=633.0, + step=10.0, + help="Красный He-Ne лазер: 633 нм", + ) + + st.sidebar.header("Апертура") + + aperture_label = st.sidebar.selectbox( + "Тип апертуры", + options=[a.value for a in ApertureType], + ) + aperture_type = ApertureType(aperture_label) + + aperture_size_mm = st.sidebar.number_input( + "Полуширина щели / радиус a (мм)", + min_value=0.01, + max_value=5.0, + value=0.5, + step=0.05, + ) + + separation_mm = 1.0 + num_slits = 5 + period_mm = 0.2 + + if aperture_type == ApertureType.DOUBLE_SLIT: + separation_mm = st.sidebar.number_input( + "Расстояние между центрами щелей d (мм)", + min_value=0.05, + max_value=10.0, + value=1.0, + step=0.1, + ) + + if aperture_type == ApertureType.GRATING: + num_slits = st.sidebar.number_input( + "Число щелей N", + min_value=2, + max_value=50, + value=5, + step=1, + ) + period_mm = st.sidebar.number_input( + "Период решётки p (мм)", + min_value=0.05, + max_value=5.0, + value=0.2, + step=0.05, + ) + + st.sidebar.header("Геометрия") + + distance_m = st.sidebar.number_input( + "Расстояние до экрана z (м)", + min_value=0.01, + max_value=100.0, + value=1.0, + step=0.1, + ) + + return DiffractionConfig( + aperture_type=aperture_type, + wavelength_nm=wavelength_nm, + aperture_size_mm=aperture_size_mm, + distance_m=distance_m, + separation_mm=separation_mm, + num_slits=int(num_slits), + period_mm=period_mm, + grid_size=512, + aperture_grid_mm=5.0, + ) diff --git a/src/models/diffraction/utils.py b/src/models/diffraction/utils.py new file mode 100644 index 0000000..28559ab --- /dev/null +++ b/src/models/diffraction/utils.py @@ -0,0 +1,90 @@ +from typing import Final + +import numpy as np +from scipy.special import j1 + +from models.diffraction.config import ApertureType, DiffractionConfig + +_ZERO_THRESHOLD: Final[float] = 1e-12 + + +def fresnel_number( + aperture_size_m: float, + wavelength_m: float, + distance_m: float, +) -> float: + return aperture_size_m**2 / (wavelength_m * distance_m) + + +def _sinc_squared( + x_m: np.ndarray, + config: DiffractionConfig, +) -> np.ndarray: + wavelength = config.wavelength_nm * 1e-9 + a = config.aperture_size_mm * 1e-3 + u = np.pi * a * x_m / (wavelength * config.distance_m) + with np.errstate(divide="ignore", invalid="ignore"): + return np.where(np.abs(u) < _ZERO_THRESHOLD, 1.0, (np.sin(u) / u) ** 2) + + +def _airy( + x_m: np.ndarray, + config: DiffractionConfig, +) -> np.ndarray: + wavelength = config.wavelength_nm * 1e-9 + r = config.aperture_size_mm * 1e-3 + k = 2 * np.pi / wavelength + v = k * r * x_m / config.distance_m + with np.errstate(divide="ignore", invalid="ignore"): + return np.where(np.abs(v) < _ZERO_THRESHOLD, 1.0, (2 * j1(v) / v) ** 2) + + +def _double_slit( + x_m: np.ndarray, + config: DiffractionConfig, +) -> np.ndarray: + wavelength = config.wavelength_nm * 1e-9 + a = config.aperture_size_mm * 1e-3 + d = config.separation_mm * 1e-3 + u = np.pi * a * x_m / (wavelength * config.distance_m) + v = np.pi * d * x_m / (wavelength * config.distance_m) + with np.errstate(divide="ignore", invalid="ignore"): + sinc = np.where(np.abs(u) < _ZERO_THRESHOLD, 1.0, (np.sin(u) / u) ** 2) + return sinc * np.cos(v) ** 2 + + +def _grating( + x_m: np.ndarray, + config: DiffractionConfig, +) -> np.ndarray: + wavelength = config.wavelength_nm * 1e-9 + a = config.aperture_size_mm * 1e-3 + p = config.period_mm * 1e-3 + n = config.num_slits + u = np.pi * a * x_m / (wavelength * config.distance_m) + v = np.pi * p * x_m / (wavelength * config.distance_m) + with np.errstate(divide="ignore", invalid="ignore"): + sinc = np.where(np.abs(u) < _ZERO_THRESHOLD, 1.0, (np.sin(u) / u) ** 2) + sin_v = np.sin(v) + grating_factor = np.where( + np.abs(sin_v) < _ZERO_THRESHOLD, + 1.0, + (np.sin(n * v) / (n * sin_v)) ** 2, + ) + return sinc * grating_factor + + +_APERTURE_FORMULAS = { + ApertureType.SINGLE_SLIT: _sinc_squared, + ApertureType.CIRCULAR: _airy, + ApertureType.DOUBLE_SLIT: _double_slit, + ApertureType.GRATING: _grating, +} + + +def analytical_intensity( + x_m: np.ndarray, + config: DiffractionConfig, + aperture_type: ApertureType, +) -> np.ndarray: + return _APERTURE_FORMULAS[aperture_type](x_m, config) diff --git a/tests/m11/__init__.py b/tests/m11/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/m11/test_diffraction.py b/tests/m11/test_diffraction.py new file mode 100644 index 0000000..ca26b0a --- /dev/null +++ b/tests/m11/test_diffraction.py @@ -0,0 +1,170 @@ +import numpy as np +import pytest + +from models.diffraction.config import ApertureType, DiffractionConfig +from models.diffraction.objects import DiffractionSimulation +from models.diffraction.utils import ( + analytical_intensity, + fresnel_number, +) + + +@pytest.fixture +def default_config() -> DiffractionConfig: + return DiffractionConfig(grid_size=64, aperture_grid_mm=2.0) + + +@pytest.fixture +def sim() -> DiffractionSimulation: + return DiffractionSimulation() + + +# --- aperture geometry --- + + +def test_single_slit_aperture_shape(default_config, sim): + n = default_config.grid_size + l_ap = default_config.aperture_grid_mm * 1e-3 + xi = np.linspace(-l_ap / 2, l_ap / 2 - l_ap / n, n) + aperture = sim.create_aperture(xi, default_config, ApertureType.SINGLE_SLIT) + assert aperture.shape == (n, n) + assert aperture[n // 2, n // 2] == 1.0 + assert aperture[n // 2, 0] == 0.0 + + +def test_circular_aperture_center_open(default_config, sim): + n = default_config.grid_size + l_ap = default_config.aperture_grid_mm * 1e-3 + xi = np.linspace(-l_ap / 2, l_ap / 2 - l_ap / n, n) + aperture = sim.create_aperture(xi, default_config, ApertureType.CIRCULAR) + assert aperture[n // 2, n // 2] == 1.0 + assert aperture[0, 0] == 0.0 + + +def test_double_slit_has_two_open_regions(sim): + config = DiffractionConfig( + aperture_type=ApertureType.DOUBLE_SLIT, + aperture_size_mm=0.05, + separation_mm=0.5, + grid_size=64, + aperture_grid_mm=2.0, + ) + n = config.grid_size + l_ap = config.aperture_grid_mm * 1e-3 + xi = np.linspace(-l_ap / 2, l_ap / 2 - l_ap / n, n) + aperture = sim.create_aperture(xi, config, ApertureType.DOUBLE_SLIT) + row = aperture[n // 2] + open_pixels = int(np.sum(row)) + assert open_pixels >= 2 + + +def test_grating_has_n_slit_groups(sim): + config = DiffractionConfig( + aperture_type=ApertureType.GRATING, + aperture_size_mm=0.04, + period_mm=0.2, + num_slits=5, + grid_size=128, + aperture_grid_mm=2.0, + ) + n = config.grid_size + l_ap = config.aperture_grid_mm * 1e-3 + xi = np.linspace(-l_ap / 2, l_ap / 2 - l_ap / n, n) + aperture = sim.create_aperture(xi, config, ApertureType.GRATING) + row = aperture[n // 2] + transitions = int(np.sum(np.abs(np.diff(row)))) + assert transitions == 2 * config.num_slits + + +# --- Fresnel number --- + + +def test_fresnel_number(): + # F = a² / (λz) = (0.5e-3)² / (633e-9 * 1.0) ≈ 0.395 + f = fresnel_number( + aperture_size_m=0.5e-3, + wavelength_m=633e-9, + distance_m=1.0, + ) + assert abs(f - 0.395) < 0.001 + + +# --- analytical formulas --- + + +def test_single_slit_peak_is_one(default_config): + x = np.array([0.0]) + i = analytical_intensity(x, default_config, ApertureType.SINGLE_SLIT) + assert abs(i[0] - 1.0) < 1e-10 + + +def test_single_slit_zeros_at_lambda_z_over_a(default_config): + wavelength = default_config.wavelength_nm * 1e-9 + a = default_config.aperture_size_mm * 1e-3 + z = default_config.distance_m + x_zero = np.array([wavelength * z / a]) + i = analytical_intensity(x_zero, default_config, ApertureType.SINGLE_SLIT) + assert i[0] < 1e-6 + + +def test_circular_peak_is_one(default_config): + x = np.array([0.0]) + i = analytical_intensity(x, default_config, ApertureType.CIRCULAR) + assert abs(i[0] - 1.0) < 1e-10 + + +def test_grating_main_maximum_at_order_1(sim): + config = DiffractionConfig( + aperture_type=ApertureType.GRATING, + wavelength_nm=633.0, + aperture_size_mm=0.05, + distance_m=1.0, + period_mm=0.2, + num_slits=5, + grid_size=64, + aperture_grid_mm=2.0, + ) + wavelength = config.wavelength_nm * 1e-9 + p = config.period_mm * 1e-3 + z = config.distance_m + x_max = np.array([wavelength * z / p]) + i = analytical_intensity(x_max, config, ApertureType.GRATING) + assert i[0] > 0.1 + + +# --- simulation output --- + + +def test_simulation_output_shapes(default_config, sim): + result = sim.simulate(default_config) + n = default_config.grid_size + assert result.intensity_2d.shape == (n, n) + assert result.intensity_1d.shape == (n,) + assert result.x_obs_m.shape == (n,) + + +def test_simulation_no_nan_or_inf(default_config, sim): + result = sim.simulate(default_config) + assert np.all(np.isfinite(result.intensity_2d)) + assert np.all(np.isfinite(result.intensity_1d)) + + +def test_simulation_intensity_nonnegative(default_config, sim): + result = sim.simulate(default_config) + assert np.all(result.intensity_2d >= 0) + assert np.all(result.intensity_1d >= 0) + + +def test_single_slit_peak_at_center(sim): + config = DiffractionConfig( + aperture_type=ApertureType.SINGLE_SLIT, + wavelength_nm=633.0, + aperture_size_mm=0.1, + distance_m=20.0, + grid_size=64, + aperture_grid_mm=2.0, + ) + result = sim.simulate(config) + n = config.grid_size + center = n // 2 + assert result.intensity_1d[center] == result.intensity_1d.max()