diff --git a/M3_report/conclusion.tex b/M3_report/conclusion.tex new file mode 100644 index 0000000..7df1e19 --- /dev/null +++ b/M3_report/conclusion.tex @@ -0,0 +1,33 @@ +\section{Заключение} + +В работе реализовано численное моделирование движения заряженной частицы в магнитной ловушке, образованной двумя соосными токовыми кольцами. + +\textbf{Основные результаты:} + +\begin{enumerate} + \item \textbf{Вычисление магнитного поля} + \begin{itemize} + \item Реализован расчёт поля кольца с током через эллиптические интегралы + \item Построена картина силовых линий в меридиональной плоскости + \item Подтверждена структура «магнитного зеркала»: максимумы поля у колец и минимум в центре + \end{itemize} + + \item \textbf{Моделирование траекторий} + \begin{itemize} + \item Реализован интегратор РК4 для уравнений движения в поле Лоренца + \item Сохранение кинетической энергии подтверждает корректность алгоритма + \item Визуализированы захваченные и транзитные траектории + \end{itemize} + + \item \textbf{Условия удержания} + \begin{itemize} + \item Численно подтверждено существование критического начального угла скорости + \item Доля захваченных частиц растёт при увеличении пробочного отношения $R_m$ + \end{itemize} +\end{enumerate} + +\subsection{Физический смысл} + +Магнитная ловушка удерживает частицы за счёт эффекта магнитного зеркала: у кольца поле сильнее, и частица, летящая к нему, «разворачивается» обратно. Частицы, летящие почти вдоль оси, зеркало не останавливает --- они вылетают через торцы ловушки. + +Именно по этому принципу работают реальные термоядерные установки типа «пробкотрон». diff --git a/M3_report/intro.tex b/M3_report/intro.tex new file mode 100644 index 0000000..adb45c4 --- /dev/null +++ b/M3_report/intro.tex @@ -0,0 +1,53 @@ +\newpage + +\section{Введение} + +\Task Два одинаковых кольца с током расположены в параллельных плоскостях, токи направлены в одну сторону. В пространстве между кольцами образуется магнитная ловушка --- область, способная удерживать заряженные частицы. Такие устройства используются, например, для удержания горячей термоядерной плазмы. Необходимо реализовать алгоритм расчёта траектории пробного заряда в такой ловушке. + +\Goal +\begin{enumerate} + \item Вычислить магнитное поле системы двух токовых колец в произвольной точке пространства + \item Реализовать численное интегрирование уравнений движения заряженной частицы + \item Визуализировать траектории частиц и картину магнитного поля + \item Исследовать, при каких начальных условиях частица удерживается в ловушке +\end{enumerate} + +\newpage +\section{Физическая постановка задачи} + +Рассматриваются два одинаковых круговых витка с током, расположенных симметрично: + +\[ + \text{Кольцо 1:} \quad z = +d/2, \quad \text{радиус } R, \quad \text{ток } I +\] +\[ + \text{Кольцо 2:} \quad z = -d/2, \quad \text{радиус } R, \quad \text{ток } I +\] + +Токи в обоих кольцах текут в одном направлении. Пробный заряд $q$ с массой $m$ движется в магнитном поле этой системы. + +\textbf{Уравнение движения:} + +На движущийся заряд действует сила Лоренца: +\[ + \vec{F} = q\,\vec{v} \times \vec{B}(\vec{r}) +\] + +Второй закон Ньютона даёт: +\[ + m\,\ddot{\vec{r}} = q\,\dot{\vec{r}} \times \vec{B}(\vec{r}) +\] + +\textbf{Важное свойство:} сила Лоренца всегда перпендикулярна скорости, поэтому она не совершает работы. Это означает, что модуль скорости (и кинетическая энергия) частицы не изменяется: +\[ + |\vec{v}| = \text{const} +\] + +\textbf{Параметры системы:} +\begin{itemize} + \item $\mu_0 = 4\pi \times 10^{-7}$ Гн/м --- магнитная постоянная + \item Радиус колец $R$ и расстояние между ними $d$ + \item Заряд $q$ и масса $m$ пробной частицы +\end{itemize} + +\newpage diff --git a/M3_report/main.tex b/M3_report/main.tex new file mode 100644 index 0000000..dcaa826 --- /dev/null +++ b/M3_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/M3_report/models/lorentz_motion.tex b/M3_report/models/lorentz_motion.tex new file mode 100644 index 0000000..5b81a33 --- /dev/null +++ b/M3_report/models/lorentz_motion.tex @@ -0,0 +1,25 @@ +\subsection{Движение заряженной частицы в магнитном поле} + +\subsubsection{Уравнения движения} + +Запишем уравнение $m\ddot{\vec{r}} = q\,\vec{v} \times \vec{B}$ как систему ОДУ первого порядка. Обозначим $\vec{r} = (x, y, z)$, $\vec{v} = (v_x, v_y, v_z)$: + +\[ +\begin{cases} +\dot{x} = v_x, \quad \dot{y} = v_y, \quad \dot{z} = v_z \\[0.5em] +\dot{v}_x = \dfrac{q}{m}(v_y B_z - v_z B_y) \\[0.5em] +\dot{v}_y = \dfrac{q}{m}(v_z B_x - v_x B_z) \\[0.5em] +\dot{v}_z = \dfrac{q}{m}(v_x B_y - v_y B_x) +\end{cases} +\] + +Магнитное поле $\vec{B} = (B_x, B_y, B_z)$ вычисляется в текущей точке $\vec{r}(t)$. + +\subsubsection{Сохранение энергии} + +Поскольку сила Лоренца не совершает работы: +\[ +\mathcal{E}_k = \frac{m|\vec{v}|^2}{2} = \text{const} +\] + +Это свойство используется как проверка точности: если $|\vec{v}|$ изменяется, значит шаг интегрирования слишком велик. diff --git a/M3_report/models/magnetic_field.tex b/M3_report/models/magnetic_field.tex new file mode 100644 index 0000000..3a9b1b2 --- /dev/null +++ b/M3_report/models/magnetic_field.tex @@ -0,0 +1,43 @@ +\subsection{Магнитное поле токового кольца} + +Поле прямого проводника или бесконечной катушки выражается простыми формулами, но для \textit{кольца} аналитического выражения в замкнутой форме не существует. Поле вычисляется через специальные функции --- полные эллиптические интегралы $K(k)$ и $E(k)$. + +\subsubsection{Поле одного кольца} + +Рассмотрим кольцо радиуса $R$ с током $I$ в плоскости $z = z_0$. Из-за осевой симметрии удобно использовать цилиндрические координаты $(\rho, z)$, где $\rho = \sqrt{x^2 + y^2}$. Ненулевые компоненты поля: + +\[ +B_\rho = \frac{\mu_0 I}{2\pi} \cdot \frac{z - z_0}{\rho \sqrt{(\rho + R)^2 + (z-z_0)^2}} +\left[ -K(k) + \frac{R^2 + \rho^2 + (z-z_0)^2}{(R - \rho)^2 + (z-z_0)^2} E(k) \right] +\] + +\[ +B_z = \frac{\mu_0 I}{2\pi} \cdot \frac{1}{\sqrt{(\rho + R)^2 + (z-z_0)^2}} +\left[ K(k) - \frac{R^2 - \rho^2 - (z-z_0)^2}{(R - \rho)^2 + (z-z_0)^2} E(k) \right] +\] + +где параметр $k$: +\[ +k^2 = \frac{4R\rho}{(\rho + R)^2 + (z - z_0)^2} +\] + +Прямо на оси ($\rho = 0$) формула упрощается до знакомого вида: +\[ +B_z\big|_{\rho=0} = \frac{\mu_0 I R^2}{2\left[R^2 + (z - z_0)^2\right]^{3/2}} +\] + +\subsubsection{Суммарное поле двух колец} + +Поле системы двух колец --- это просто сумма полей каждого из них (принцип суперпозиции): +\[ +\vec{B}(\vec{r}) = \vec{B}^{(1)}(\vec{r}) + \vec{B}^{(2)}(\vec{r}) +\] + +\subsubsection{Почему это ловушка?} + +При одинаковых направлениях токов поле вблизи каждого кольца усиливается, а в центре между ними --- ослабевает. Отношение максимального поля к минимальному называется \textbf{пробочным отношением}: +\[ +R_m = \frac{B_{\max}}{B_{\min}} > 1 +\] + +Именно это неравенство и обеспечивает удержание частиц. diff --git a/M3_report/models/models.tex b/M3_report/models/models.tex new file mode 100644 index 0000000..9183a1e --- /dev/null +++ b/M3_report/models/models.tex @@ -0,0 +1,8 @@ +\section{Физические модели} + +\input{models/magnetic_field} +\newpage +\input{models/lorentz_motion} +\newpage +\input{models/trapping_condition} +\newpage diff --git a/M3_report/models/trapping_condition.tex b/M3_report/models/trapping_condition.tex new file mode 100644 index 0000000..3946eb0 --- /dev/null +++ b/M3_report/models/trapping_condition.tex @@ -0,0 +1,17 @@ +\subsection{Условия удержания в магнитной ловушке} + +\subsubsection{Почему частица может отразиться?} + +Когда частица движется в сторону более сильного поля (к кольцу), её поперечная скорость растёт, а продольная убывает --- так как $|\vec{v}|$ сохраняется. Если поле достаточно сильное, продольная скорость обращается в ноль и частица отражается назад --- как мяч, брошенный вверх. + +Это и есть \textbf{эффект магнитного зеркала}. + +\subsubsection{Условие захвата} + +Частица \textbf{захватывается} ловушкой, если она отражается от обоих колец. Это зависит от соотношения продольной и поперечной составляющих начальной скорости. Частицы, летящие почти вдоль оси (преимущественно продольная скорость), пролетают сквозь зеркала --- они \textbf{теряются}. + +Ключевой параметр --- \textbf{пробочное отношение} $R_m = B_{\max}/B_{\min}$: чем оно больше, тем более широкий диапазон начальных условий приводит к захвату. + +\subsubsection{Статистика траекторий} + +Для изучения захвата проводится серия расчётов со случайными начальными условиями. Для каждой траектории определяется исход: \textbf{захвачена} или \textbf{потеряна}. По результатам строится зависимость доли захваченных частиц от начального направления скорости и пробочного отношения $R_m$. diff --git a/M3_report/numerical_methods.tex b/M3_report/numerical_methods.tex new file mode 100644 index 0000000..435c7fd --- /dev/null +++ b/M3_report/numerical_methods.tex @@ -0,0 +1,73 @@ +\section{Численные методы решения} + +\subsection{Вычисление магнитного поля} + +\subsubsection{Перевод координат} + +Для вычисления поля кольца в точке $(x, y, z)$ сначала находим цилиндрический радиус: +\[ +\rho = \sqrt{x^2 + y^2} +\] + +Затем по формулам из раздела 2 находим компоненты $(B_\rho, B_z)$ и переводим в декартовы: +\[ +B_x = B_\rho \cdot \frac{x}{\rho}, \qquad B_y = B_\rho \cdot \frac{y}{\rho} +\] + +При $\rho = 0$ (точка на оси) используется упрощённая осевая формула, чтобы не делить на ноль. + +\subsubsection{Эллиптические интегралы} + +Значения $K(k)$ и $E(k)$ вычисляются встроенными функциями \texttt{scipy.special.ellipk} и \texttt{scipy.special.ellipe}. + +\subsection{Интегрирование уравнений движения} + +\subsubsection{Метод Рунге--Кутты 4-го порядка (РК4)} + +Система ОДУ решается методом РК4. Обозначим вектор состояния: +\[ +\vec{u} = (x,\; y,\; z,\; v_x,\; v_y,\; v_z)^\top +\] + +Один шаг метода: +\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*} + +где $\vec{f}(\vec{u})$ --- правая часть системы из раздела~2.2. + +\subsubsection{Выбор шага} + +Шаг интегрирования должен быть значительно меньше характерного периода вращения частицы в поле: +\[ +T_c = \frac{2\pi m}{|q| B} +\] +\[ +h = \frac{T_c}{N_{\text{шаг}}}, \quad N_{\text{шаг}} \approx 100 +\] + +\subsubsection{Проверка точности} + +После каждого шага проверяем, что модуль скорости не изменился: +\[ +\frac{\bigl||\vec{v}_{n+1}| - |\vec{v}_0|\bigr|}{|\vec{v}_0|} < 10^{-5} +\] + +Если это условие нарушается, шаг нужно уменьшить. + +\subsection{Классификация траекторий} + +\begin{enumerate} + \item Задаём начальные условия: положение $\vec{r}_0$ и скорость $\vec{v}_0$ + \item Интегрируем уравнения движения + \item На каждом шаге проверяем: вышла ли частица за пределы ловушки ($|z| > z_{\text{bound}}$) + \begin{itemize} + \item Если да --- частица \textbf{потеряна} + \item Если время симуляции истекло, а частица внутри --- частица \textbf{захвачена} + \end{itemize} + \item Повторяем для набора начальных условий и строим статистику +\end{enumerate} diff --git a/M3_report/preamble.tex b/M3_report/preamble.tex new file mode 100644 index 0000000..428cac8 --- /dev/null +++ b/M3_report/preamble.tex @@ -0,0 +1,103 @@ +\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} +\pgfplotsset{width=10cm,compat=1.18} +\pgfkeys{/pgf/trig format=rad} + + +\usepackage{xcolor} +\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}{} + +\usepackage{hyperref} +\hypersetup{ + colorlinks=true, + linkcolor=blue, + filecolor=blue, + urlcolor=blue, + citecolor=blue +} + + +\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{Примеры.} } + +\newpagestyle{main}{ + \setheadrule{.4pt} + \sethead{}{\bullet \; \textit{Цифровизация физических процессов} \; \bullet}{} + \setfootrule{.4pt} + \setfoot{ВШПИ МФТИ}{М3. Магнитная пробка}{Б13-402} +} + +\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} diff --git a/M3_report/report.pdf b/M3_report/report.pdf new file mode 100644 index 0000000..1505cef Binary files /dev/null and b/M3_report/report.pdf differ diff --git a/M3_report/results.tex b/M3_report/results.tex new file mode 100644 index 0000000..bf459c6 --- /dev/null +++ b/M3_report/results.tex @@ -0,0 +1,33 @@ +\section{Результаты} + +\subsection{Магнитное поле системы} + +Рассчитанное поле демонстрирует ожидаемую структуру магнитной ловушки: +\begin{itemize} + \item На оси $z$ поле имеет два максимума (у каждого кольца) и минимум в центре + \item Пробочное отношение $R_m = B_{\max}/B_{\min}$ растёт при увеличении расстояния между кольцами + \item При $d \approx R$ достигается $R_m \approx 1.5 \div 2$ +\end{itemize} + +\subsection{Траектории заряженных частиц} + +Моделирование показывает три характерных типа движения: + +\begin{enumerate} + \item \textbf{Захваченная частица} (большая поперечная составляющая скорости): \\ + Частица «прыгает» между двумя зеркалами, не покидая ловушку. + + \item \textbf{Транзитная частица} (скорость направлена преимущественно вдоль оси): \\ + Частица проходит сквозь оба кольца и вылетает из системы. + + \item \textbf{Пограничный случай}: \\ + Частица несколько раз отражается, но постепенно уходит из ловушки. +\end{enumerate} + +\subsection{Статистика захвата} + +При моделировании ансамбля частиц со случайными начальными условиями: +\begin{itemize} + \item Частицы с большой поперечной составляющей скорости захватываются чаще + \item При бо́льшем пробочном отношении $R_m$ ловушка удерживает более широкий набор начальных условий +\end{itemize} diff --git a/M3_report/title.tex b/M3_report/title.tex new file mode 100644 index 0000000..ab76a9c --- /dev/null +++ b/M3_report/title.tex @@ -0,0 +1,18 @@ +\begin{titlepage} + \centering + {\scshape\Large Московский физико-технический институт \par} + \vspace{1cm} + {\scshape\Large Высшая школа программной инженерии \par} + \vspace{2cm} + {\huge \textbf{М3. Магнитная пробка} \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..181f8b6 100644 --- a/src/main.py +++ b/src/main.py @@ -2,6 +2,7 @@ from models.billiard import page as billiard_page from models.ideal_gas import page as ideal_gas_page +from models.magnetic_trap import page as magnetic_trap_page from models.pendulum import page as pendulum_page from models.rolling_ball import page as rolling_ball_page from models.stone_flight import page as stone_flight_page @@ -13,6 +14,11 @@ def main() -> None: stone_flight_page, title="M1 проект", url_path="stone_flight" ), streamlit.Page(billiard_page, title="M2 Бильярд", url_path="billiard"), + streamlit.Page( + magnetic_trap_page, + title="M3 Магнитная ловушка", + url_path="magnetic_trap", + ), streamlit.Page( rolling_ball_page, title="M4 Шар на столе", url_path="rolling_ball" ), diff --git a/src/models/magnetic_trap/__init__.py b/src/models/magnetic_trap/__init__.py new file mode 100644 index 0000000..4a4457e --- /dev/null +++ b/src/models/magnetic_trap/__init__.py @@ -0,0 +1,245 @@ +import numpy as np +import streamlit as st + +from models.magnetic_trap.charts import ( + create_axial_field_chart, + create_field_magnitude_chart, + create_speed_chart, + create_trajectory_xy_chart, + create_trajectory_xz_chart, + create_trapping_statistics_chart, + create_z_vs_time_chart, +) +from models.magnetic_trap.config import PHYSICAL_CONSTANTS +from models.magnetic_trap.objects import MagneticTrapModel, TrapParameters +from models.magnetic_trap.sidebar import create_magnetic_trap_sidebar +from models.magnetic_trap.utils import ( + analyze_speed_conservation, + classify_trajectory, + compute_axial_field_profile, + compute_field_grid, + mirror_ratio, + trapping_statistics, +) + +MAX_GYROPERIODS_PER_SIM = 450 +MAX_GYROPERIODS_PER_BATCH = 1300 +MIN_TRAJECTORY_POINTS = 2 + + +def _estimate_cyclotron_count(params: TrapParameters) -> tuple[float, float]: + b_estimate = ( + PHYSICAL_CONSTANTS["mu_0"] + * params.current + * params.ring_radius**2 + / ((params.ring_radius**2 + (params.ring_distance / 2) ** 2) ** 1.5) + ) + if b_estimate <= 0 or abs(params.charge) <= 0: + return 0.0, float("inf") + period = 2 * np.pi * params.mass / (abs(params.charge) * b_estimate) + return b_estimate, params.t_max / period + + +def _check_simulation_cost(params: TrapParameters, budget: float) -> bool: + _, n_periods = _estimate_cyclotron_count(params) + if n_periods > budget: + st.error( + f"Параметры приведут к слишком жёсткому уравнению" + f" (≈{n_periods:.0f} ларморовских периодов на интервале t_max," + f" допустимо ≤ {int(budget)}). Уменьшите ток I, t_max или массу" + f" частицы — либо увеличьте радиус кольца." + ) + return False + return True + + +def page() -> None: + st.set_page_config(page_title="М3. Магнитная ловушка", layout="wide") + st.title("М3. Магнитная ловушка из двух токовых колец") + st.write( + "Численное моделирование движения заряженной частицы в магнитной" + " ловушке, образованной двумя соосными кольцами с током." + ) + + model = MagneticTrapModel() + params = create_magnetic_trap_sidebar() + + analysis_type = st.sidebar.radio( + "Тип анализа:", + [ + "Магнитное поле", + "Траектория частицы", + "Статистика захвата", + ], + ) + + try: + if analysis_type == "Магнитное поле": + display_field_analysis(model, params) + elif analysis_type == "Траектория частицы": + display_trajectory_analysis(model, params) + elif analysis_type == "Статистика захвата": + display_trapping_statistics(model, params) + except Exception as e: + st.error(f"Ошибка при моделировании: {e}") + + +def display_field_analysis( + model: MagneticTrapModel, params: TrapParameters +) -> None: + st.subheader("Поле на оси системы") + z_values, b_z = compute_axial_field_profile(model, params) + st.altair_chart( + create_axial_field_chart(z_values, b_z), + use_container_width=True, + ) + + st.subheader("Картина поля в меридиональной плоскости") + rho_mesh, z_mesh, b_x, b_z_grid = compute_field_grid(model, params) + st.altair_chart( + create_field_magnitude_chart(rho_mesh, z_mesh, b_x, b_z_grid), + use_container_width=True, + ) + + st.subheader("Пробочное отношение") + ratio_data = mirror_ratio(model, params) + col1, col2, col3 = st.columns(3) + with col1: + st.metric("B_max (Тл)", f"{ratio_data['b_max']:.4e}") + with col2: + st.metric("B_min (Тл)", f"{ratio_data['b_min']:.4e}") + with col3: + st.metric("R_m = B_max / B_min", f"{ratio_data['mirror_ratio']:.3f}") + + +def display_trajectory_analysis( + model: MagneticTrapModel, params: TrapParameters +) -> None: + if not _check_simulation_cost(params, MAX_GYROPERIODS_PER_SIM): + return + x, y, z, vx, vy, vz, t = model.simulate_motion(params) + + if len(t) < MIN_TRAJECTORY_POINTS: + st.error( + f"Симуляция вернула {len(t)} точек." + f" Уменьшите dt или увеличьте t_max." + ) + return + + classification = classify_trajectory(z, t, params) + + col1, col2 = st.columns(2) + with col1: + st.subheader("Траектория в плоскости xz") + st.altair_chart( + create_trajectory_xz_chart(x, z, t), + use_container_width=True, + ) + with col2: + st.subheader("Траектория в плоскости xy") + st.altair_chart( + create_trajectory_xy_chart(x, y, t), + use_container_width=True, + ) + + st.subheader("Положение по z во времени") + st.altair_chart( + create_z_vs_time_chart(z, t, params.z_bound), + use_container_width=True, + ) + + st.subheader("Сохранение модуля скорости") + speed = np.sqrt(vx**2 + vy**2 + vz**2) + st.altair_chart( + create_speed_chart(speed, t), + use_container_width=True, + ) + + col3, col4 = st.columns(2) + with col3: + st.subheader("Исход") + outcome = classification["outcome"] + if outcome == "trapped": + st.success("Частица захвачена ловушкой") + elif outcome == "lost": + st.warning( + f"Частица покинула ловушку в момент" + f" t = {classification['exit_time']:.3e} с" + ) + else: + st.info("Нет данных") + + with col4: + st.subheader("Сохранение |v|") + analysis = analyze_speed_conservation(vx, vy, vz) + if "error" in analysis: + st.write(analysis["error"]) + else: + st.write( + f"**Начальный модуль скорости:**" + f" {analysis['initial_speed']:.4e} м/с" + ) + st.write( + f"**Конечный модуль скорости:**" + f" {analysis['final_speed']:.4e} м/с" + ) + st.write( + f"**Максимальное относительное отклонение:**" + f" {analysis['max_deviation']:.3e}" + ) + st.write( + f"**Качество сохранения:** {analysis['conservation_quality']}" + ) + + +def display_trapping_statistics( # noqa: PLR0914 + model: MagneticTrapModel, base_params: TrapParameters +) -> None: + st.subheader("Зависимость захвата от питч-угла") + + if not _check_simulation_cost(base_params, MAX_GYROPERIODS_PER_BATCH / 9): + return + + speed = float( + np.sqrt(base_params.vx0**2 + base_params.vy0**2 + base_params.vz0**2) + ) + if speed <= 0: + speed = 5.0e4 + + pitch_angles = np.linspace(5.0, 85.0, 9) + + progress = st.progress(0) + outcomes = [] + for i, angle_deg in enumerate(pitch_angles): + angle_rad = np.radians(angle_deg) + v_par = speed * np.cos(angle_rad) + v_perp = speed * np.sin(angle_rad) + + params = TrapParameters(**{ + **base_params.__dict__, + "vx0": v_perp, + "vy0": 0.0, + "vz0": v_par, + }) + + _x, _y, z, _vx, _vy, _vz, t = model.simulate_motion(params) + result = classify_trajectory(z, t, params) + outcomes.append((float(angle_deg), result["outcome"])) + progress.progress((i + 1) / len(pitch_angles)) + + st.altair_chart( + create_trapping_statistics_chart(outcomes), + use_container_width=True, + ) + + stats = trapping_statistics(model, base_params, pitch_angles, speed) + col1, col2, col3 = st.columns(3) + with col1: + st.metric("Захвачено", stats["trapped_count"]) + with col2: + st.metric("Потеряно", stats["lost_count"]) + with col3: + st.metric( + "Доля захваченных", + f"{stats['trapped_fraction'] * 100:.1f}%", + ) diff --git a/src/models/magnetic_trap/charts.py b/src/models/magnetic_trap/charts.py new file mode 100644 index 0000000..c4ce56b --- /dev/null +++ b/src/models/magnetic_trap/charts.py @@ -0,0 +1,196 @@ +import altair as alt +import numpy as np +import pandas as pd + + +def create_axial_field_chart( + z_values: np.ndarray, b_z: np.ndarray +) -> alt.Chart: + df = pd.DataFrame({"z": z_values, "B_z": b_z, "|B_z|": np.abs(b_z)}) + + return ( + alt.Chart(df) + .mark_line(strokeWidth=2, color="darkblue") + .add_params(alt.selection_interval(bind="scales")) + .encode( + x=alt.X("z:Q", title="z (м)"), + y=alt.Y("B_z:Q", title="B_z (Тл)"), + tooltip=["z:Q", "B_z:Q", "|B_z|:Q"], + ) + .properties( + title="Магнитное поле на оси системы", + width=600, + height=300, + ) + ) + + +def _trajectory_layer( + df: pd.DataFrame, + x_field: str, + y_field: str, + x_title: str, + y_title: str, +) -> alt.LayerChart: + base = alt.Chart(df).encode( + x=alt.X(f"{x_field}:Q", title=x_title), + y=alt.Y(f"{y_field}:Q", title=y_title), + ) + + line = base.mark_line(strokeWidth=1.0, color="lightgray").encode( + order=alt.Order("t:Q"), + ) + + dots = base.mark_circle(size=12, opacity=0.7).encode( + color=alt.Color( + "t:Q", + title="Время (с)", + scale=alt.Scale(scheme="viridis"), + ), + tooltip=[f"{x_field}:Q", f"{y_field}:Q", "t:Q"], + ) + + return (line + dots).add_params(alt.selection_interval(bind="scales")) + + +def create_trajectory_xz_chart( + x: np.ndarray, z: np.ndarray, t: np.ndarray +) -> alt.LayerChart: + df = pd.DataFrame({"x": x, "z": z, "t": t}) + return _trajectory_layer(df, "x", "z", "x (м)", "z (м)").properties( + title="Траектория в плоскости xz", width=400, height=400 + ) + + +def create_trajectory_xy_chart( + x: np.ndarray, y: np.ndarray, t: np.ndarray +) -> alt.LayerChart: + df = pd.DataFrame({"x": x, "y": y, "t": t}) + return _trajectory_layer(df, "x", "y", "x (м)", "y (м)").properties( + title="Траектория в плоскости xy", width=400, height=400 + ) + + +def create_speed_chart(speed: np.ndarray, t: np.ndarray) -> alt.Chart: + df = pd.DataFrame({"t": t, "|v|": speed}) + + return ( + alt.Chart(df) + .mark_line(strokeWidth=2, color="green") + .add_params(alt.selection_interval(bind="scales")) + .encode( + x=alt.X("t:Q", title="Время t (с)"), + y=alt.Y("|v|:Q", title="|v| (м/с)"), + tooltip=["t:Q", "|v|:Q"], + ) + .properties(title="Сохранение модуля скорости", width=600, height=250) + ) + + +def create_z_vs_time_chart( + z: np.ndarray, t: np.ndarray, z_bound: float +) -> alt.Chart: + df = pd.DataFrame({"t": t, "z": z}) + bounds_df = pd.DataFrame({ + "t": [t.min(), t.max(), t.min(), t.max()], + "z": [z_bound, z_bound, -z_bound, -z_bound], + "label": [ + "Граница (+)", + "Граница (+)", + "Граница (-)", + "Граница (-)", + ], + }) + + line = ( + alt.Chart(df) + .mark_line(strokeWidth=2, color="navy") + .encode( + x=alt.X("t:Q", title="Время t (с)"), + y=alt.Y("z:Q", title="z (м)"), + tooltip=["t:Q", "z:Q"], + ) + ) + + bound = ( + alt.Chart(bounds_df) + .mark_line(strokeDash=[5, 5], color="red", strokeWidth=1.5) + .encode( + x="t:Q", + y="z:Q", + detail="label:N", + ) + ) + + return (line + bound).properties( + title="Положение по z во времени", width=600, height=300 + ) + + +def create_field_magnitude_chart( + rho_mesh: np.ndarray, + z_mesh: np.ndarray, + b_x: np.ndarray, + b_z: np.ndarray, +) -> alt.Chart: + b_magnitude = np.sqrt(b_x**2 + b_z**2) + + df = pd.DataFrame({ + "rho": rho_mesh.flatten(), + "z": z_mesh.flatten(), + "|B|": b_magnitude.flatten(), + }) + + return ( + alt.Chart(df) + .mark_rect() + .encode( + x=alt.X("rho:Q", title="ρ (м)", bin=alt.Bin(maxbins=40)), + y=alt.Y("z:Q", title="z (м)", bin=alt.Bin(maxbins=40)), + color=alt.Color( + "|B|:Q", + title="|B| (Тл)", + scale=alt.Scale(scheme="inferno"), + ), + tooltip=["rho:Q", "z:Q", "|B|:Q"], + ) + .properties( + title="Картина модуля поля в меридиональной плоскости", + width=400, + height=400, + ) + ) + + +def create_trapping_statistics_chart(outcomes: list) -> alt.Chart: + df = pd.DataFrame(outcomes, columns=["pitch_angle_deg", "outcome"]) + df["captured"] = (df["outcome"] == "trapped").astype(int) + + return ( + alt.Chart(df) + .mark_circle(size=80, opacity=0.8) + .encode( + x=alt.X( + "pitch_angle_deg:Q", + title="Питч-угол (°)", + ), + y=alt.Y( + "captured:Q", + title="Захвачена (1) / Потеряна (0)", + scale=alt.Scale(domain=[-0.1, 1.1]), + ), + color=alt.Color( + "outcome:N", + scale=alt.Scale( + domain=["trapped", "lost"], + range=["green", "red"], + ), + ), + tooltip=["pitch_angle_deg:Q", "outcome:N"], + ) + .properties( + title="Зависимость захвата от питч-угла", + width=600, + height=300, + ) + ) diff --git a/src/models/magnetic_trap/config.py b/src/models/magnetic_trap/config.py new file mode 100644 index 0000000..86f3c25 --- /dev/null +++ b/src/models/magnetic_trap/config.py @@ -0,0 +1,26 @@ +import numpy as np + +PHYSICAL_CONSTANTS = { + "mu_0": 4 * np.pi * 1e-7, + "elementary_charge": 1.602176634e-19, + "electron_mass": 9.1093837015e-31, + "proton_mass": 1.67262192369e-27, +} + +PARTICLE_PRESETS = { + "electron": { + "description": "Электрон", + "charge": -PHYSICAL_CONSTANTS["elementary_charge"], + "mass": PHYSICAL_CONSTANTS["electron_mass"], + }, + "proton": { + "description": "Протон", + "charge": PHYSICAL_CONSTANTS["elementary_charge"], + "mass": PHYSICAL_CONSTANTS["proton_mass"], + }, + "custom": { + "description": "Произвольная частица", + "charge": PHYSICAL_CONSTANTS["elementary_charge"], + "mass": PHYSICAL_CONSTANTS["proton_mass"], + }, +} diff --git a/src/models/magnetic_trap/objects.py b/src/models/magnetic_trap/objects.py new file mode 100644 index 0000000..268faa6 --- /dev/null +++ b/src/models/magnetic_trap/objects.py @@ -0,0 +1,217 @@ +from dataclasses import dataclass + +import numpy as np +import streamlit as st +from scipy import integrate, special + +from models.magnetic_trap.config import PARTICLE_PRESETS, PHYSICAL_CONSTANTS + +EPSILON = 1e-12 + + +@dataclass +class TrapParameters: + ring_radius: float = 0.1 + ring_distance: float = 0.1 + current: float = 1.0e5 + particle_type: str = "proton" + charge: float = PARTICLE_PRESETS["proton"]["charge"] + mass: float = PARTICLE_PRESETS["proton"]["mass"] + x0: float = 0.0 + y0: float = 0.0 + z0: float = 0.0 + vx0: float = 0.0 + vy0: float = 5.0e4 + vz0: float = 1.0e4 + t_max: float = 1.0e-5 + dt: float = 1.0e-8 + z_bound: float = 0.2 + + +class MagneticTrapModel: + def __init__(self) -> None: + self.mu_0 = PHYSICAL_CONSTANTS["mu_0"] + + def _ring_field_cylindrical( # noqa: PLR0914 + self, + rho: np.ndarray | float, + z: np.ndarray | float, + ring_radius: float, + ring_z: float, + current: float, + ) -> tuple[np.ndarray, np.ndarray]: + rho = np.asarray(rho, dtype=float) + z = np.asarray(z, dtype=float) + dz = z - ring_z + + on_axis = rho < EPSILON + + b_rho = np.zeros_like(rho, dtype=float) + b_z = np.zeros_like(rho, dtype=float) + + if np.any(on_axis): + denom = (ring_radius**2 + dz**2) ** 1.5 + b_z_axis = self.mu_0 * current * ring_radius**2 / (2 * denom) + b_z = np.where(on_axis, b_z_axis, b_z) + + off_axis = ~on_axis + if np.any(off_axis): + rho_off = np.where(off_axis, rho, ring_radius) + dz_off = np.where(off_axis, dz, 0.0) + + sum_term = (rho_off + ring_radius) ** 2 + dz_off**2 + diff_term = (ring_radius - rho_off) ** 2 + dz_off**2 + + sum_term = np.maximum(sum_term, EPSILON) + diff_term = np.maximum(diff_term, EPSILON) + + k_squared = 4 * ring_radius * rho_off / sum_term + k_squared = np.clip(k_squared, 0.0, 1.0 - EPSILON) + + ellip_k = special.ellipk(k_squared) + ellip_e = special.ellipe(k_squared) + + prefactor = self.mu_0 * current / (2 * np.pi) + + b_rho_val = ( + prefactor + * dz_off + / (rho_off * np.sqrt(sum_term)) + * ( + -ellip_k + + (ring_radius**2 + rho_off**2 + dz_off**2) + / diff_term + * ellip_e + ) + ) + + b_z_val = ( + prefactor + / np.sqrt(sum_term) + * ( + ellip_k + + (ring_radius**2 - rho_off**2 - dz_off**2) + / diff_term + * ellip_e + ) + ) + + b_rho = np.where(off_axis, b_rho_val, b_rho) + b_z = np.where(off_axis, b_z_val, b_z) + + return b_rho, b_z + + def magnetic_field( # noqa: PLR0914 + self, + x: np.ndarray | float, + y: np.ndarray | float, + z: np.ndarray | float, + params: TrapParameters, + ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + x_arr = np.asarray(x, dtype=float) + y_arr = np.asarray(y, dtype=float) + z_arr = np.asarray(z, dtype=float) + + rho = np.sqrt(x_arr**2 + y_arr**2) + ring_z_top = params.ring_distance / 2 + ring_z_bottom = -params.ring_distance / 2 + + b_rho_top, b_z_top = self._ring_field_cylindrical( + rho, z_arr, params.ring_radius, ring_z_top, params.current + ) + b_rho_bot, b_z_bot = self._ring_field_cylindrical( + rho, z_arr, params.ring_radius, ring_z_bottom, params.current + ) + + b_rho_total = b_rho_top + b_rho_bot + b_z_total = b_z_top + b_z_bot + + on_axis = rho < EPSILON + rho_safe = np.where(on_axis, 1.0, rho) + cos_phi = np.where(on_axis, 0.0, x_arr / rho_safe) + sin_phi = np.where(on_axis, 0.0, y_arr / rho_safe) + + b_x = b_rho_total * cos_phi + b_y = b_rho_total * sin_phi + + return b_x, b_y, b_z_total + + def _equations_of_motion( # noqa: PLR0914 + self, + _t: float, + state: np.ndarray, + params: TrapParameters, + ) -> np.ndarray: + x, y, z, vx, vy, vz = state + + b_x, b_y, b_z = self.magnetic_field(x, y, z, params) + + b_x_f = float(b_x) + b_y_f = float(b_y) + b_z_f = float(b_z) + + q_over_m = params.charge / params.mass + + ax = q_over_m * (vy * b_z_f - vz * b_y_f) + ay = q_over_m * (vz * b_x_f - vx * b_z_f) + az = q_over_m * (vx * b_y_f - vy * b_x_f) + + return np.array([vx, vy, vz, ax, ay, az]) + + @st.cache_data(ttl=300) + def simulate_motion( + _self, # noqa: N805 + params: TrapParameters, + ) -> tuple[ + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + ]: + initial_state = np.array([ + params.x0, + params.y0, + params.z0, + params.vx0, + params.vy0, + params.vz0, + ]) + + t_eval = np.arange(0, params.t_max, params.dt) + + def hit_boundary(_t: float, state: np.ndarray) -> float: + return params.z_bound - abs(state[2]) + + hit_boundary.terminal = True # type: ignore[attr-defined] + hit_boundary.direction = -1.0 # type: ignore[attr-defined] + + solution = integrate.solve_ivp( + fun=lambda t, state: _self._equations_of_motion(t, state, params), # noqa: SLF001 + t_span=(0, params.t_max), + y0=initial_state, + t_eval=t_eval, + method="RK45", + rtol=1e-6, + atol=1e-9, + events=hit_boundary, + ) + + x = solution.y[0] + y = solution.y[1] + z = solution.y[2] + vx = solution.y[3] + vy = solution.y[4] + vz = solution.y[5] + t = solution.t + + return x, y, z, vx, vy, vz, t + + def cyclotron_period( # noqa: PLR6301 + self, params: TrapParameters, b_magnitude: float + ) -> float: + if b_magnitude == 0.0 or params.charge == 0.0: + return float("inf") + return 2 * np.pi * params.mass / (abs(params.charge) * b_magnitude) diff --git a/src/models/magnetic_trap/sidebar.py b/src/models/magnetic_trap/sidebar.py new file mode 100644 index 0000000..8ade5df --- /dev/null +++ b/src/models/magnetic_trap/sidebar.py @@ -0,0 +1,191 @@ +import streamlit as st + +from models.magnetic_trap.config import PARTICLE_PRESETS +from models.magnetic_trap.objects import TrapParameters + + +def create_magnetic_trap_sidebar() -> TrapParameters: # noqa: PLR0914 + st.sidebar.header("Параметры колец") + + ring_radius = st.sidebar.number_input( + "Радиус кольца R (м)", + min_value=0.02, + max_value=1.0, + value=0.1, + step=0.01, + format="%.3f", + ) + + ring_distance = st.sidebar.number_input( + "Расстояние между кольцами d (м)", + min_value=0.02, + max_value=1.0, + value=0.2, + step=0.01, + format="%.3f", + ) + + current = st.sidebar.number_input( + "Ток I (А)", + min_value=1.0e3, + max_value=2.0e5, + value=1.0e5, + step=1000.0, + format="%.1f", + help=( + "Большие токи увеличивают жёсткость уравнений — расчёт замедляется" + ), + ) + + st.sidebar.header("Пробная частица") + + particle_type = st.sidebar.selectbox( + "Тип частицы", + options=["proton", "custom"], + format_func=lambda x: PARTICLE_PRESETS[x]["description"], + index=0, + help=( + "Электрон не предлагается: масса в 1836 раз меньше протона" + " — расчёт занял бы минуты" + ), + ) + + preset = PARTICLE_PRESETS[particle_type] + if particle_type == "custom": + charge = st.sidebar.number_input( + "Заряд q (Кл)", + min_value=-3.0e-19, + max_value=3.0e-19, + value=float(preset["charge"]), + format="%.3e", + ) + mass = st.sidebar.number_input( + "Масса m (кг)", + min_value=1.0e-27, + max_value=1.0e-24, + value=float(preset["mass"]), + format="%.3e", + help=( + "Минимум близок к массе протона:" + " лёгкие частицы расчёт замедляют" + ), + ) + else: + charge = float(preset["charge"]) + mass = float(preset["mass"]) + st.sidebar.write(f"Заряд: {charge:.3e} Кл") + st.sidebar.write(f"Масса: {mass:.3e} кг") + + st.sidebar.header("Начальное положение") + + pos_bound = max(ring_radius, ring_distance / 2) + x0 = st.sidebar.number_input( + "x₀ (м)", + min_value=-pos_bound, + max_value=pos_bound, + value=0.0, + step=0.01, + format="%.4f", + ) + y0 = st.sidebar.number_input( + "y₀ (м)", + min_value=-pos_bound, + max_value=pos_bound, + value=0.0, + step=0.01, + format="%.4f", + ) + z0 = st.sidebar.number_input( + "z₀ (м)", + min_value=-pos_bound, + max_value=pos_bound, + value=0.0, + step=0.01, + format="%.4f", + ) + + st.sidebar.header("Начальная скорость") + + v_max = 1.0e7 + vx0 = st.sidebar.number_input( + "v_x (м/с)", + min_value=-v_max, + max_value=v_max, + value=0.0, + step=1.0e3, + format="%.3e", + ) + vy0 = st.sidebar.number_input( + "v_y (м/с)", + min_value=-v_max, + max_value=v_max, + value=5.0e4, + step=1.0e3, + format="%.3e", + ) + vz0 = st.sidebar.number_input( + "v_z (м/с)", + min_value=-v_max, + max_value=v_max, + value=1.0e4, + step=1.0e3, + format="%.3e", + ) + + st.sidebar.header("Параметры симуляции") + + t_max_ns = st.sidebar.number_input( + "Время симуляции (нс)", + min_value=100.0, + max_value=100_000.0, + value=10_000.0, + step=100.0, + format="%.0f", + help="Прямо влияет на время расчёта. 1 нс = 10⁻⁹ с", + ) + + dt_max_ns = max(1.0, t_max_ns / 100) + dt_default_ns = max(1.0, min(10.0, dt_max_ns)) + dt_ns = st.sidebar.number_input( + "Шаг времени dt (нс)", + min_value=1.0, + max_value=dt_max_ns, + value=dt_default_ns, + step=1.0, + format="%.0f", + help=( + "Шаг сохранения точек траектории (не шаг интегратора RK45)." + " Автоматически ограничивается ≤ t_max / 100," + " чтобы траектория содержала минимум 100 точек." + ), + ) + + t_max = t_max_ns * 1e-9 + dt = dt_ns * 1e-9 + + z_bound = st.sidebar.number_input( + "Граница ловушки |z| (м)", + min_value=ring_distance / 2, + max_value=2.0, + value=max(2 * ring_distance, ring_distance / 2), + step=0.01, + format="%.3f", + ) + + return TrapParameters( + ring_radius=ring_radius, + ring_distance=ring_distance, + current=current, + particle_type=particle_type, + charge=charge, + mass=mass, + x0=x0, + y0=y0, + z0=z0, + vx0=vx0, + vy0=vy0, + vz0=vz0, + t_max=t_max, + dt=dt, + z_bound=z_bound, + ) diff --git a/src/models/magnetic_trap/utils.py b/src/models/magnetic_trap/utils.py new file mode 100644 index 0000000..e3a6dd5 --- /dev/null +++ b/src/models/magnetic_trap/utils.py @@ -0,0 +1,166 @@ +import numpy as np + +from models.magnetic_trap.objects import MagneticTrapModel, TrapParameters + +EPSILON = 1e-15 +EXCELLENT_THRESHOLD = 1e-6 +GOOD_THRESHOLD = 1e-4 +ACCEPTABLE_THRESHOLD = 1e-2 + + +def compute_axial_field_profile( + model: MagneticTrapModel, + params: TrapParameters, + n_points: int = 200, + z_range: float | None = None, +) -> tuple[np.ndarray, np.ndarray]: + if z_range is None: + z_range = params.ring_distance * 1.5 + z_values = np.linspace(-z_range, z_range, n_points) + zeros = np.zeros_like(z_values) + _, _, b_z = model.magnetic_field(zeros, zeros, z_values, params) + return z_values, b_z + + +def mirror_ratio( + model: MagneticTrapModel, + params: TrapParameters, + n_points: int = 400, +) -> dict: + z_values, b_z = compute_axial_field_profile( + model, params, n_points=n_points, z_range=params.ring_distance / 2 + ) + b_abs = np.abs(b_z) + b_max = float(np.max(b_abs)) + b_min = float(np.min(b_abs)) + + ratio = b_max / b_min if b_min > EPSILON else float("inf") + + return { + "b_max": b_max, + "b_min": b_min, + "mirror_ratio": ratio, + "z_at_max": float(z_values[int(np.argmax(b_abs))]), + "z_at_min": float(z_values[int(np.argmin(b_abs))]), + } + + +def classify_trajectory( + z: np.ndarray, t: np.ndarray, params: TrapParameters +) -> dict: + if len(z) == 0: + return {"outcome": "no_data", "exit_time": 0.0} + + escape_indices = np.where(np.abs(z) >= params.z_bound)[0] + if len(escape_indices) > 0: + exit_idx = int(escape_indices[0]) + return { + "outcome": "lost", + "exit_time": float(t[exit_idx]), + "exit_z": float(z[exit_idx]), + } + + finished = bool(t[-1] >= params.t_max - params.dt) + if finished: + return { + "outcome": "trapped", + "exit_time": float(t[-1]), + "exit_z": float(z[-1]), + } + return { + "outcome": "lost", + "exit_time": float(t[-1]), + "exit_z": float(z[-1]), + } + + +def analyze_speed_conservation( + vx: np.ndarray, vy: np.ndarray, vz: np.ndarray +) -> dict: + if len(vx) == 0: + return {"error": "Нет данных для анализа"} + + speed = np.sqrt(vx**2 + vy**2 + vz**2) + v0 = float(speed[0]) + v_final = float(speed[-1]) + + relative_change = (v_final - v0) / v0 if v0 > EPSILON else 0.0 + max_deviation = ( + float(np.max(np.abs(speed - v0)) / v0) if v0 > EPSILON else 0.0 + ) + + if max_deviation < EXCELLENT_THRESHOLD: + quality = "excellent" + elif max_deviation < GOOD_THRESHOLD: + quality = "good" + elif max_deviation < ACCEPTABLE_THRESHOLD: + quality = "acceptable" + else: + quality = "poor" + + return { + "initial_speed": v0, + "final_speed": v_final, + "relative_change": relative_change, + "max_deviation": max_deviation, + "mean_speed": float(np.mean(speed)), + "speed_std": float(np.std(speed)), + "conservation_quality": quality, + } + + +def compute_field_grid( + model: MagneticTrapModel, + params: TrapParameters, + n_rho: int = 30, + n_z: int = 30, + rho_max: float | None = None, + z_max: float | None = None, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + if rho_max is None: + rho_max = params.ring_radius * 1.5 + if z_max is None: + z_max = params.ring_distance + rho_grid = np.linspace(EPSILON, rho_max, n_rho) + z_grid = np.linspace(-z_max, z_max, n_z) + + rho_mesh, z_mesh = np.meshgrid(rho_grid, z_grid) + zeros = np.zeros_like(rho_mesh) + + b_x, _, b_z = model.magnetic_field(rho_mesh, zeros, z_mesh, params) + + return rho_mesh, z_mesh, b_x, b_z + + +def trapping_statistics( + model: MagneticTrapModel, + base_params: TrapParameters, + pitch_angles_deg: np.ndarray, + speed: float, +) -> dict: + outcomes = [] + for angle_deg in pitch_angles_deg: + angle_rad = np.radians(angle_deg) + v_par = speed * np.cos(angle_rad) + v_perp = speed * np.sin(angle_rad) + + params = TrapParameters(**{ + **base_params.__dict__, + "vx0": v_perp, + "vy0": 0.0, + "vz0": v_par, + }) + + _x, _y, z, _vx, _vy, _vz, t = model.simulate_motion(params) + result = classify_trajectory(z, t, params) + outcomes.append((float(angle_deg), result["outcome"])) + + trapped = [a for a, o in outcomes if o == "trapped"] + lost = [a for a, o in outcomes if o == "lost"] + + return { + "outcomes": outcomes, + "trapped_count": len(trapped), + "lost_count": len(lost), + "trapped_fraction": len(trapped) / len(outcomes) if outcomes else 0.0, + } diff --git a/tests/m3/__init__.py b/tests/m3/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/m3/test_magnetic_trap_physics.py b/tests/m3/test_magnetic_trap_physics.py new file mode 100644 index 0000000..03fd69e --- /dev/null +++ b/tests/m3/test_magnetic_trap_physics.py @@ -0,0 +1,574 @@ +import numpy as np +import pytest + +from models.magnetic_trap.config import PARTICLE_PRESETS, PHYSICAL_CONSTANTS +from models.magnetic_trap.objects import MagneticTrapModel, TrapParameters +from models.magnetic_trap.utils import ( + analyze_speed_conservation, + classify_trajectory, + compute_axial_field_profile, + compute_field_grid, + mirror_ratio, + trapping_statistics, +) + + +@pytest.fixture +def trap_model() -> MagneticTrapModel: + return MagneticTrapModel() + + +@pytest.fixture +def basic_params() -> TrapParameters: + return TrapParameters( + ring_radius=0.1, + ring_distance=0.1, + current=1.0e5, + particle_type="proton", + charge=PARTICLE_PRESETS["proton"]["charge"], + mass=PARTICLE_PRESETS["proton"]["mass"], + x0=0.0, + y0=0.0, + z0=0.0, + vx0=0.0, + vy0=5.0e4, + vz0=1.0e4, + t_max=1.0e-5, + dt=1.0e-8, + z_bound=0.5, + ) + + +@pytest.fixture +def transit_params() -> TrapParameters: + return TrapParameters( + ring_radius=0.1, + ring_distance=0.1, + current=1.0e5, + particle_type="proton", + charge=PARTICLE_PRESETS["proton"]["charge"], + mass=PARTICLE_PRESETS["proton"]["mass"], + x0=0.0, + y0=0.0, + z0=0.0, + vx0=0.0, + vy0=1.0e3, + vz0=1.0e6, + t_max=1.0e-6, + dt=1.0e-9, + z_bound=0.15, + ) + + +def test_model_initialization(trap_model: MagneticTrapModel) -> None: + assert trap_model.mu_0 == pytest.approx(PHYSICAL_CONSTANTS["mu_0"]) + + +def test_axial_field_single_ring_formula(trap_model: MagneticTrapModel) -> None: + ring_radius = 0.1 + ring_z = 0.0 + current = 1.0e5 + + z_test = np.array([0.0, 0.05, 0.1, -0.1]) + + rho = np.zeros_like(z_test) + b_rho, b_z = trap_model._ring_field_cylindrical( + rho, z_test, ring_radius, ring_z, current + ) + + expected_b_z = ( + trap_model.mu_0 + * current + * ring_radius**2 + / (2 * (ring_radius**2 + (z_test - ring_z) ** 2) ** 1.5) + ) + + assert np.allclose(b_rho, 0.0, atol=1e-12) + assert np.allclose(b_z, expected_b_z, rtol=1e-9) + + +def test_field_at_center_of_two_rings( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + b_x, b_y, b_z = trap_model.magnetic_field(0.0, 0.0, 0.0, basic_params) + + half_d = basic_params.ring_distance / 2 + expected_single = ( + trap_model.mu_0 + * basic_params.current + * basic_params.ring_radius**2 + / (2 * (basic_params.ring_radius**2 + half_d**2) ** 1.5) + ) + expected_total = 2 * expected_single + + assert float(b_x) == pytest.approx(0.0, abs=1e-12) + assert float(b_y) == pytest.approx(0.0, abs=1e-12) + assert float(b_z) == pytest.approx(expected_total, rel=1e-9) + + +def test_field_symmetry_along_axis( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + z_positive = 0.02 + z_negative = -0.02 + + _, _, b_z_pos = trap_model.magnetic_field( + 0.0, 0.0, z_positive, basic_params + ) + _, _, b_z_neg = trap_model.magnetic_field( + 0.0, 0.0, z_negative, basic_params + ) + + assert float(b_z_pos) == pytest.approx(float(b_z_neg), rel=1e-10) + + +def test_mirror_ratio_greater_than_one( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + params = TrapParameters(**{ + **basic_params.__dict__, + "ring_distance": 3 * basic_params.ring_radius, + }) + ratio_data = mirror_ratio(trap_model, params) + + assert ratio_data["b_max"] > ratio_data["b_min"] + assert ratio_data["mirror_ratio"] > 1.0 + + +def test_field_max_near_rings( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + mirror_params = TrapParameters(**{ + **basic_params.__dict__, + "ring_distance": 3 * basic_params.ring_radius, + }) + z_axis = np.linspace( + -mirror_params.ring_distance / 2, + mirror_params.ring_distance / 2, + 201, + ) + zeros = np.zeros_like(z_axis) + _, _, b_z = trap_model.magnetic_field(zeros, zeros, z_axis, mirror_params) + b_abs = np.abs(b_z) + + center_idx = len(z_axis) // 2 + + b_at_center = b_abs[center_idx] + b_at_edges = max(b_abs[0], b_abs[-1]) + + assert b_at_edges > b_at_center + + +def test_lorentz_force_perpendicular_to_velocity( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + state = np.array([0.05, 0.0, 0.02, 1.0e4, 2.0e4, 3.0e4]) + derivatives = trap_model._equations_of_motion(0.0, state, basic_params) + + velocity = state[3:] + acceleration = derivatives[3:] + + dot_product = float(np.dot(velocity, acceleration)) + speed = float(np.linalg.norm(velocity)) + accel = float(np.linalg.norm(acceleration)) + + if accel > 0: + cos_angle = dot_product / (speed * accel) + assert abs(cos_angle) < 1e-9 + + +def test_zero_charge_means_no_force( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + params = TrapParameters(**{**basic_params.__dict__, "charge": 0.0}) + state = np.array([0.05, 0.0, 0.02, 1.0e4, 2.0e4, 3.0e4]) + derivatives = trap_model._equations_of_motion(0.0, state, params) + + assert derivatives[0] == pytest.approx(state[3]) + assert derivatives[1] == pytest.approx(state[4]) + assert derivatives[2] == pytest.approx(state[5]) + assert derivatives[3] == pytest.approx(0.0) + assert derivatives[4] == pytest.approx(0.0) + assert derivatives[5] == pytest.approx(0.0) + + +def test_speed_conservation_during_simulation( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + _x, _y, _z, vx, vy, vz, t = trap_model.simulate_motion(basic_params) + + assert len(t) > 0 + + analysis = analyze_speed_conservation(vx, vy, vz) + + assert analysis["max_deviation"] < 1e-3 + + +def test_simulation_initial_conditions( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + x, y, z, vx, vy, vz, _t = trap_model.simulate_motion(basic_params) + + assert x[0] == pytest.approx(basic_params.x0) + assert y[0] == pytest.approx(basic_params.y0) + assert z[0] == pytest.approx(basic_params.z0) + assert vx[0] == pytest.approx(basic_params.vx0) + assert vy[0] == pytest.approx(basic_params.vy0) + assert vz[0] == pytest.approx(basic_params.vz0) + + +def test_trapped_particle_classification( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + _x, _y, z, _vx, _vy, _vz, t = trap_model.simulate_motion(basic_params) + classification = classify_trajectory(z, t, basic_params) + + assert classification["outcome"] == "trapped" + assert abs(classification["exit_z"]) < basic_params.z_bound + + +def test_lost_particle_classification( + trap_model: MagneticTrapModel, transit_params: TrapParameters +) -> None: + _x, _y, z, _vx, _vy, _vz, t = trap_model.simulate_motion(transit_params) + classification = classify_trajectory(z, t, transit_params) + + assert classification["outcome"] == "lost" + + +def test_field_off_axis_has_radial_component( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + b_x_axis, b_y_axis, _ = trap_model.magnetic_field( + 0.0, 0.0, 0.02, basic_params + ) + assert float(b_x_axis) == pytest.approx(0.0, abs=1e-12) + assert float(b_y_axis) == pytest.approx(0.0, abs=1e-12) + + rho = 0.03 + b_x, b_y, _ = trap_model.magnetic_field(rho, 0.0, 0.02, basic_params) + assert abs(float(b_x)) > 1e-9 + assert float(b_y) == pytest.approx(0.0, abs=1e-12) + + +def test_cyclotron_period_decreases_with_field( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + period_weak = trap_model.cyclotron_period(basic_params, 0.001) + period_strong = trap_model.cyclotron_period(basic_params, 1.0) + + assert period_strong < period_weak + + +def test_zero_current_zero_field( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + params = TrapParameters(**{**basic_params.__dict__, "current": 0.0}) + b_x, b_y, b_z = trap_model.magnetic_field(0.05, 0.0, 0.02, params) + + assert float(b_x) == pytest.approx(0.0, abs=1e-15) + assert float(b_y) == pytest.approx(0.0, abs=1e-15) + assert float(b_z) == pytest.approx(0.0, abs=1e-15) + + +def test_field_scales_linearly_with_current( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + p1 = TrapParameters(**{**basic_params.__dict__, "current": 1.0e5}) + p2 = TrapParameters(**{**basic_params.__dict__, "current": 3.0e5}) + + _, _, b1 = trap_model.magnetic_field(0.0, 0.0, 0.0, p1) + _, _, b2 = trap_model.magnetic_field(0.0, 0.0, 0.0, p2) + + assert float(b2) == pytest.approx(3.0 * float(b1), rel=1e-9) + + +def test_field_reverses_with_negative_current( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + p_pos = TrapParameters(**{**basic_params.__dict__, "current": 1.0e5}) + p_neg = TrapParameters(**{**basic_params.__dict__, "current": -1.0e5}) + + _, _, b_pos = trap_model.magnetic_field(0.03, 0.0, 0.02, p_pos) + _, _, b_neg = trap_model.magnetic_field(0.03, 0.0, 0.02, p_neg) + + assert float(b_neg) == pytest.approx(-float(b_pos), rel=1e-9) + + +def test_field_decays_far_from_rings( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + _, _, b_close = trap_model.magnetic_field(0.0, 0.0, 0.05, basic_params) + _, _, b_far = trap_model.magnetic_field(0.0, 0.0, 5.0, basic_params) + + assert abs(float(b_far)) < abs(float(b_close)) / 1000 + + +def test_field_at_center_decreases_with_ring_radius( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + p_small = TrapParameters(**{ + **basic_params.__dict__, + "ring_radius": 0.05, + }) + p_big = TrapParameters(**{ + **basic_params.__dict__, + "ring_radius": 0.5, + }) + + _, _, b_small = trap_model.magnetic_field(0.0, 0.0, 0.0, p_small) + _, _, b_big = trap_model.magnetic_field(0.0, 0.0, 0.0, p_big) + + assert abs(float(b_small)) > abs(float(b_big)) + + +def test_field_anti_symmetric_b_rho_about_z_zero( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + rho = 0.04 + b_x_top, _, _ = trap_model.magnetic_field(rho, 0.0, 0.03, basic_params) + b_x_bot, _, _ = trap_model.magnetic_field(rho, 0.0, -0.03, basic_params) + + assert float(b_x_top) == pytest.approx(-float(b_x_bot), rel=1e-9) + + +def test_field_axial_symmetric_b_z_about_z_zero( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + rho = 0.04 + _, _, b_z_top = trap_model.magnetic_field(rho, 0.0, 0.03, basic_params) + _, _, b_z_bot = trap_model.magnetic_field(rho, 0.0, -0.03, basic_params) + + assert float(b_z_top) == pytest.approx(float(b_z_bot), rel=1e-9) + + +def test_field_rotational_symmetry_around_axis( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + rho = 0.04 + z = 0.03 + b_x_a, b_y_a, b_z_a = trap_model.magnetic_field(rho, 0.0, z, basic_params) + b_x_b, b_y_b, b_z_b = trap_model.magnetic_field(0.0, rho, z, basic_params) + + b_rho_a = float(b_x_a) + b_rho_b = float(b_y_b) + + assert b_rho_a == pytest.approx(b_rho_b, rel=1e-9) + assert float(b_z_a) == pytest.approx(float(b_z_b), rel=1e-9) + assert float(b_y_a) == pytest.approx(0.0, abs=1e-12) + assert float(b_x_b) == pytest.approx(0.0, abs=1e-12) + + +def test_field_grid_shapes( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + n_rho, n_z = 25, 30 + rho_mesh, z_mesh, b_x, b_z = compute_field_grid( + trap_model, basic_params, n_rho=n_rho, n_z=n_z + ) + + assert rho_mesh.shape == (n_z, n_rho) + assert z_mesh.shape == (n_z, n_rho) + assert b_x.shape == (n_z, n_rho) + assert b_z.shape == (n_z, n_rho) + + +def test_axial_field_profile_endpoints( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + n_points = 101 + z_range = 0.5 + z, b_z = compute_axial_field_profile( + trap_model, basic_params, n_points=n_points, z_range=z_range + ) + + assert z.shape == (n_points,) + assert b_z.shape == (n_points,) + assert z[0] == pytest.approx(-z_range) + assert z[-1] == pytest.approx(z_range) + + +def test_helmholtz_has_lower_mirror_ratio_than_real_trap( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + helmholtz = TrapParameters(**{ + **basic_params.__dict__, + "ring_distance": basic_params.ring_radius, + }) + real_trap = TrapParameters(**{ + **basic_params.__dict__, + "ring_distance": 4 * basic_params.ring_radius, + }) + + r_helmholtz = mirror_ratio(trap_model, helmholtz)["mirror_ratio"] + r_real_trap = mirror_ratio(trap_model, real_trap)["mirror_ratio"] + + assert r_real_trap > r_helmholtz + + +def test_classify_trajectory_handles_empty_arrays( + basic_params: TrapParameters, +) -> None: + result = classify_trajectory(np.array([]), np.array([]), basic_params) + assert result["outcome"] == "no_data" + + +def test_analyze_speed_conservation_perfect( + basic_params: TrapParameters, +) -> None: + n = 100 + vx = np.full(n, 3.0) + vy = np.full(n, 4.0) + vz = np.zeros(n) + + analysis = analyze_speed_conservation(vx, vy, vz) + + assert analysis["initial_speed"] == pytest.approx(5.0) + assert analysis["final_speed"] == pytest.approx(5.0) + assert analysis["max_deviation"] == pytest.approx(0.0) + assert analysis["conservation_quality"] == "excellent" + + +def test_analyze_speed_conservation_handles_empty() -> None: + analysis = analyze_speed_conservation( + np.array([]), np.array([]), np.array([]) + ) + assert "error" in analysis + + +def test_analyze_speed_conservation_quality_levels() -> None: + n = 50 + base = np.ones(n) + + perturbed = base.copy() + perturbed[25] = 1.0 + 5e-3 + analysis = analyze_speed_conservation(perturbed, np.zeros(n), np.zeros(n)) + assert analysis["conservation_quality"] in {"acceptable", "poor"} + + +def test_particle_charge_sign_reverses_force_direction( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + state = np.array([0.05, 0.0, 0.02, 1.0e4, 2.0e4, 3.0e4]) + + pos = TrapParameters(**{**basic_params.__dict__, "charge": 1.6e-19}) + neg = TrapParameters(**{**basic_params.__dict__, "charge": -1.6e-19}) + + deriv_pos = trap_model._equations_of_motion(0.0, state, pos) + deriv_neg = trap_model._equations_of_motion(0.0, state, neg) + + assert deriv_pos[3] == pytest.approx(-deriv_neg[3], rel=1e-9) + assert deriv_pos[4] == pytest.approx(-deriv_neg[4], rel=1e-9) + assert deriv_pos[5] == pytest.approx(-deriv_neg[5], rel=1e-9) + + +def test_particle_in_uniform_b_does_circular_motion( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + huge_radius_params = TrapParameters(**{ + **basic_params.__dict__, + "ring_radius": 5.0, + "ring_distance": 10.0, + "current": 1.0e5, + "vx0": 1.0e4, + "vy0": 0.0, + "vz0": 0.0, + "t_max": 5.0e-7, + "dt": 5.0e-10, + "z_bound": 1.0, + }) + + _x, _y, z, _vx, _vy, _vz, t = trap_model.simulate_motion(huge_radius_params) + + assert len(t) > 0 + assert float(np.max(np.abs(z))) < 1.0e-3 + + +def test_trapping_statistics_returns_full_distribution( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + angles = np.array([10.0, 30.0, 60.0, 80.0]) + speed = 1.0e5 + + stats = trapping_statistics(trap_model, basic_params, angles, speed) + + assert len(stats["outcomes"]) == len(angles) + assert stats["trapped_count"] + stats["lost_count"] == len(angles) + assert 0.0 <= stats["trapped_fraction"] <= 1.0 + + +def test_trapping_statistics_more_perpendicular_more_trapped( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + mirror_params = TrapParameters(**{ + **basic_params.__dict__, + "ring_distance": 3 * basic_params.ring_radius, + "t_max": 5.0e-6, + "dt": 5.0e-9, + }) + + perpendicular = trapping_statistics( + trap_model, + mirror_params, + np.array([80.0, 85.0, 88.0]), + speed=5.0e4, + ) + parallel = trapping_statistics( + trap_model, + mirror_params, + np.array([2.0, 5.0, 10.0]), + speed=5.0e4, + ) + + assert perpendicular["trapped_fraction"] >= parallel["trapped_fraction"] + + +def test_particle_presets_consistent_signs() -> None: + assert PARTICLE_PRESETS["proton"]["charge"] > 0 + assert PARTICLE_PRESETS["electron"]["charge"] < 0 + assert ( + PARTICLE_PRESETS["proton"]["mass"] + > PARTICLE_PRESETS["electron"]["mass"] + ) + + +def test_physical_constants_correct_values() -> None: + expected_mu_0 = 4 * np.pi * 1e-7 + assert PHYSICAL_CONSTANTS["mu_0"] == pytest.approx(expected_mu_0) + assert PHYSICAL_CONSTANTS["elementary_charge"] == pytest.approx( + 1.602176634e-19, rel=1e-9 + ) + + +def test_simulation_with_zero_current_is_straight_line( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + params = TrapParameters(**{ + **basic_params.__dict__, + "current": 0.0, + "vx0": 1.0e4, + "vy0": 0.0, + "vz0": 0.0, + "t_max": 1.0e-6, + "dt": 1.0e-9, + "z_bound": 1.0, + }) + + x, y, z, _vx, _vy, _vz, t = trap_model.simulate_motion(params) + + expected_x = params.vx0 * t + assert np.allclose(x, expected_x, atol=1e-9) + assert np.allclose(y, 0.0, atol=1e-9) + assert np.allclose(z, 0.0, atol=1e-9) + + +def test_field_off_axis_has_correct_radial_direction( + trap_model: MagneticTrapModel, basic_params: TrapParameters +) -> None: + rho = 0.04 + z_above = 0.03 + + b_x_pos, _, _ = trap_model.magnetic_field(rho, 0.0, z_above, basic_params) + b_x_neg, _, _ = trap_model.magnetic_field(-rho, 0.0, z_above, basic_params) + + assert float(b_x_pos) == pytest.approx(-float(b_x_neg), rel=1e-9)