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..f241509 100644 --- a/src/main.py +++ b/src/main.py @@ -5,6 +5,7 @@ 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 +from models.tunnel_diode import page as tunnel_diode_page def main() -> None: @@ -19,6 +20,11 @@ def main() -> None: streamlit.Page( pendulum_page, title="M5 Физический маятник", url_path="pendulum" ), + streamlit.Page( + tunnel_diode_page, + title="M5 Туннельный диод (генератор)", + url_path="tunnel_diode", + ), streamlit.Page( ideal_gas_page, title="M7Б Статистика идеального газа со столкновениями частиц", diff --git a/src/models/tunnel_diode/__init__.py b/src/models/tunnel_diode/__init__.py new file mode 100644 index 0000000..8503916 --- /dev/null +++ b/src/models/tunnel_diode/__init__.py @@ -0,0 +1,199 @@ +import numpy as np +import streamlit as st + +from models.tunnel_diode.charts import ( + create_il_time_chart, + create_iv_curve_chart, + create_phase_portrait_chart, + create_spectrum_chart, + create_uc_time_chart, +) +from models.tunnel_diode.objects import TunnelDiodeModel, TunnelDiodeParameters +from models.tunnel_diode.sidebar import create_tunnel_diode_sidebar +from models.tunnel_diode.utils import ( + REGIME_DAMPED, + REGIME_GROWING, + REGIME_LIMIT_CYCLE, + REGIME_RELAXATION, + classify_regime, + compute_spectrum, + is_self_oscillating, + oscillation_metrics, + total_harmonic_distortion, +) + +REGIME_LABELS = { + REGIME_DAMPED: "Затухающие колебания", + REGIME_LIMIT_CYCLE: "Предельный цикл (генерация)", + REGIME_RELAXATION: "Релаксационные колебания", + REGIME_GROWING: "Нарастающие колебания", +} + +THD_PURE = 0.05 +THD_RELAX = 0.30 +MIN_TIME_POINTS = 2 + + +def page() -> None: + st.set_page_config( + page_title="М5. Туннельный диод (генератор)", layout="wide" + ) + st.title("М5. Автогенератор на туннельном диоде") + st.write( + "Численное моделирование LC-контура с туннельным диодом." + " ВАХ диода аппроксимируется кубической параболой;" + " на участке отрицательного дифференциального сопротивления" + " возникают незатухающие колебания." + ) + + model = TunnelDiodeModel() + params = create_tunnel_diode_sidebar() + + analysis_type = st.sidebar.radio( + "Тип анализа:", + [ + "Рабочая точка и ВАХ", + "Колебания во времени", + "Фазовый портрет и спектр", + ], + ) + + try: + if analysis_type == "Рабочая точка и ВАХ": + display_iv_analysis(model, params) + elif analysis_type == "Колебания во времени": + display_time_analysis(model, params) + elif analysis_type == "Фазовый портрет и спектр": + display_phase_analysis(model, params) + except Exception as e: + st.error(f"Ошибка при моделировании: {e}") + + +def display_iv_analysis( + model: TunnelDiodeModel, params: TunnelDiodeParameters +) -> None: + st.subheader("Вольт-амперная характеристика") + + u_grid = np.linspace(0.0, max(params.voltage_source, 3.0), 400) + i_grid = np.asarray(model.iv_curve(u_grid)) + + try: + u_star, i_star = model.find_operating_point(params) + except ValueError as e: + st.error(f"Не удалось найти рабочую точку: {e}") + return + + load_u = np.linspace(0.0, params.voltage_source, 50) + load_i = (params.voltage_source - load_u) / params.resistance + + chart = create_iv_curve_chart( + u_grid, i_grid, (u_star, i_star), (load_u, load_i) + ) + st.altair_chart(chart, use_container_width=True) + + slope = float(model.iv_derivative(u_star)) + crit_r = model.critical_resistance(params) + omega_0 = model.natural_frequency(params) + q_factor = model.quality_factor(params) + damping = model.effective_damping(params) + + col1, col2, col3 = st.columns(3) + with col1: + st.metric("Рабочая точка U*", f"{u_star:.4f} В") + st.metric("Рабочая точка I*", f"{i_star * 1e3:.4f} мА") + with col2: + st.metric("I'_D(U*)", f"{slope * 1e3:.4f} мА/В") + st.metric("Критическое R", f"{crit_r:.2f} Ом") + with col3: + st.metric("ω₀", f"{omega_0:.3e} рад/с") + st.metric("Добротность Q", f"{q_factor:.3f}") + + st.write(f"**Эффективное затухание γ_eff:** {damping:.3e}") + if is_self_oscillating(model, params): + st.success( + "γ_eff < 0 — выполнено условие автогенерации." + " R меньше критического." + ) + else: + st.info( + "γ_eff ≥ 0 — система устойчива к малым возмущениям," + " автогенерации нет." + ) + + +def display_time_analysis( + model: TunnelDiodeModel, params: TunnelDiodeParameters +) -> None: + u_c, i_l, t = model.simulate(params) + if len(t) < MIN_TIME_POINTS: + st.error("Симуляция не вернула данных") + return + + st.subheader("U_C(t)") + st.altair_chart(create_uc_time_chart(u_c, t), use_container_width=True) + + st.subheader("I_L(t)") + st.altair_chart(create_il_time_chart(i_l, t), use_container_width=True) + + metrics = oscillation_metrics(u_c, t) + if "error" in metrics: + st.warning(metrics["error"]) + return + + regime = classify_regime(model, params, u_c, t) + st.subheader("Анализ режима") + col1, col2, col3 = st.columns(3) + with col1: + st.metric("Режим", REGIME_LABELS.get(regime, regime)) + with col2: + st.metric("Амплитуда колебаний", f"{metrics['amplitude']:.4f} В") + with col3: + if metrics["period"] > 0: + st.metric("Период колебаний", f"{metrics['period'] * 1e6:.3f} мкс") + else: + st.metric("Период колебаний", "—") + + st.write( + f"**Среднее U_C:** {metrics['u_mean']:.4f} В," + f" **min/max:** {metrics['u_min']:.4f} / {metrics['u_max']:.4f} В" + ) + + +def display_phase_analysis( + model: TunnelDiodeModel, params: TunnelDiodeParameters +) -> None: + u_c, i_l, t = model.simulate(params) + if len(t) < MIN_TIME_POINTS: + st.error("Симуляция не вернула данных") + return + + try: + operating = model.find_operating_point(params) + except ValueError as e: + st.error(f"Не удалось найти рабочую точку: {e}") + return + + st.subheader("Фазовый портрет") + st.altair_chart( + create_phase_portrait_chart(u_c, i_l, operating), + use_container_width=True, + ) + + st.subheader("Спектр") + freqs, spectrum = compute_spectrum(u_c, t) + st.altair_chart( + create_spectrum_chart(freqs, spectrum), + use_container_width=True, + ) + + thd = total_harmonic_distortion(freqs, spectrum) + col1, col2 = st.columns(2) + with col1: + st.metric("THD (гармоники / основная)", f"{thd:.4f}") + with col2: + if thd < THD_PURE: + st.success("Чистая синусоида") + elif thd > THD_RELAX: + st.warning("Сильно нелинейная форма (релаксационные колебания)") + else: + st.info("Умеренная нелинейность") diff --git a/src/models/tunnel_diode/charts.py b/src/models/tunnel_diode/charts.py new file mode 100644 index 0000000..f7198dd --- /dev/null +++ b/src/models/tunnel_diode/charts.py @@ -0,0 +1,182 @@ +import altair as alt +import numpy as np +import pandas as pd + + +def create_iv_curve_chart( + u_grid: np.ndarray, + i_grid: np.ndarray, + operating_point: tuple[float, float], + load_line: tuple[np.ndarray, np.ndarray], +) -> alt.LayerChart: + iv_df = pd.DataFrame({ + "U_D (В)": u_grid, + "I_D (мА)": i_grid * 1e3, + "type": "ВАХ диода", + }) + load_df = pd.DataFrame({ + "U_D (В)": load_line[0], + "I_D (мА)": load_line[1] * 1e3, + "type": "Нагрузочная прямая", + }) + point_df = pd.DataFrame({ + "U_D (В)": [operating_point[0]], + "I_D (мА)": [operating_point[1] * 1e3], + }) + + full_df = pd.concat([iv_df, load_df], ignore_index=True) + + lines = ( + alt.Chart(full_df) + .mark_line(strokeWidth=2.0) + .encode( + x=alt.X("U_D (В):Q", title="U_D (В)"), + y=alt.Y("I_D (мА):Q", title="I_D (мА)"), + color=alt.Color( + "type:N", + scale=alt.Scale(scheme="category10"), + title="", + ), + ) + ) + + point = ( + alt.Chart(point_df) + .mark_point(filled=True, size=140, color="red", shape="diamond") + .encode( + x="U_D (В):Q", + y="I_D (мА):Q", + tooltip=["U_D (В):Q", "I_D (мА):Q"], + ) + ) + + return (lines + point).properties( + title="ВАХ диода и рабочая точка", + width=560, + height=320, + ) + + +def create_uc_time_chart(u_c: np.ndarray, t: np.ndarray) -> alt.Chart: + df = pd.DataFrame({"t (мкс)": t * 1e6, "U_C (В)": u_c}) + + return ( + alt.Chart(df) + .mark_line(strokeWidth=1.5, color="navy") + .add_params(alt.selection_interval(bind="scales")) + .encode( + x=alt.X("t (мкс):Q", title="Время t (мкс)"), + y=alt.Y("U_C (В):Q", title="U_C (В)"), + tooltip=["t (мкс):Q", "U_C (В):Q"], + ) + .properties( + title="Напряжение на конденсаторе U_C(t)", + width=720, + height=280, + ) + ) + + +def create_il_time_chart(i_l: np.ndarray, t: np.ndarray) -> alt.Chart: + df = pd.DataFrame({"t (мкс)": t * 1e6, "I_L (мА)": i_l * 1e3}) + + return ( + alt.Chart(df) + .mark_line(strokeWidth=1.5, color="darkgreen") + .add_params(alt.selection_interval(bind="scales")) + .encode( + x=alt.X("t (мкс):Q", title="Время t (мкс)"), + y=alt.Y("I_L (мА):Q", title="I_L (мА)"), + tooltip=["t (мкс):Q", "I_L (мА):Q"], + ) + .properties( + title="Ток через катушку I_L(t)", + width=720, + height=280, + ) + ) + + +def create_phase_portrait_chart( + u_c: np.ndarray, + i_l: np.ndarray, + operating_point: tuple[float, float], +) -> alt.LayerChart: + df = pd.DataFrame({ + "U_C (В)": u_c, + "I_L (мА)": i_l * 1e3, + "step": np.arange(len(u_c)), + }) + op_df = pd.DataFrame({ + "U_C (В)": [operating_point[0]], + "I_L (мА)": [operating_point[1] * 1e3], + }) + + line = ( + alt.Chart(df) + .mark_line(strokeWidth=1.0, color="lightgray") + .encode( + x=alt.X("U_C (В):Q", title="U_C (В)"), + y=alt.Y("I_L (мА):Q", title="I_L (мА)"), + order="step:Q", + ) + ) + + dots = ( + alt.Chart(df) + .mark_circle(size=10, opacity=0.6) + .encode( + x="U_C (В):Q", + y="I_L (мА):Q", + color=alt.Color( + "step:Q", + title="Шаг", + scale=alt.Scale(scheme="viridis"), + ), + tooltip=["U_C (В):Q", "I_L (мА):Q", "step:Q"], + ) + ) + + op_point = ( + alt.Chart(op_df) + .mark_point(filled=True, size=180, color="red", shape="diamond") + .encode(x="U_C (В):Q", y="I_L (мА):Q") + ) + + return (line + dots + op_point).properties( + title="Фазовый портрет (U_C, I_L)", + width=420, + height=420, + ) + + +def create_spectrum_chart(freqs: np.ndarray, spectrum: np.ndarray) -> alt.Chart: + if len(freqs) == 0 or len(spectrum) == 0: + return alt.Chart(pd.DataFrame({"f (кГц)": [], "Амплитуда": []})) + + peak = float(np.max(spectrum)) + spectrum_norm = spectrum / peak if peak > 0 else spectrum + df = pd.DataFrame({ + "f (кГц)": freqs / 1e3, + "Амплитуда": spectrum_norm, + }) + + return ( + alt.Chart(df) + .mark_line(strokeWidth=1.2, color="purple") + .add_params(alt.selection_interval(bind="scales")) + .encode( + x=alt.X("f (кГц):Q", title="Частота (кГц)"), + y=alt.Y( + "Амплитуда:Q", + title="Нормированная амплитуда", + scale=alt.Scale(domain=[0, 1.05]), + ), + tooltip=["f (кГц):Q", "Амплитуда:Q"], + ) + .properties( + title="Спектр напряжения U_C", + width=720, + height=280, + ) + ) diff --git a/src/models/tunnel_diode/config.py b/src/models/tunnel_diode/config.py new file mode 100644 index 0000000..3c46fe7 --- /dev/null +++ b/src/models/tunnel_diode/config.py @@ -0,0 +1,12 @@ +DIODE_COEFFICIENTS = { + "a": 4.0e-3, + "b": -16.0e-3, + "c": 17.0e-3, +} + +DEFAULT_CIRCUIT = { + "voltage_source": 1.5, + "resistance": 50.0, + "inductance": 1.0e-3, + "capacitance": 1.0e-7, +} diff --git a/src/models/tunnel_diode/objects.py b/src/models/tunnel_diode/objects.py new file mode 100644 index 0000000..8fc5f10 --- /dev/null +++ b/src/models/tunnel_diode/objects.py @@ -0,0 +1,146 @@ +from dataclasses import dataclass + +import numpy as np +import streamlit as st +from scipy import integrate, optimize + +from models.tunnel_diode.config import DEFAULT_CIRCUIT, DIODE_COEFFICIENTS + + +@dataclass +class TunnelDiodeParameters: + voltage_source: float = DEFAULT_CIRCUIT["voltage_source"] + resistance: float = DEFAULT_CIRCUIT["resistance"] + inductance: float = DEFAULT_CIRCUIT["inductance"] + capacitance: float = DEFAULT_CIRCUIT["capacitance"] + u_c0: float = 0.0 + i_l0: float = 0.0 + t_max: float = 5.0e-4 + n_points: int = 5000 + + +class TunnelDiodeModel: + def __init__(self) -> None: + self.a = DIODE_COEFFICIENTS["a"] + self.b = DIODE_COEFFICIENTS["b"] + self.c = DIODE_COEFFICIENTS["c"] + + def iv_curve(self, u_d: np.ndarray | float) -> np.ndarray | float: + u = np.asarray(u_d, dtype=float) + positive = u > 0 + result = np.where( + positive, + self.a * u**3 + self.b * u**2 + self.c * u, + 0.0, + ) + if np.isscalar(u_d): + return float(result) + return result + + def iv_derivative(self, u_d: np.ndarray | float) -> np.ndarray | float: + u = np.asarray(u_d, dtype=float) + positive = u > 0 + result = np.where( + positive, + 3 * self.a * u**2 + 2 * self.b * u + self.c, + 0.0, + ) + if np.isscalar(u_d): + return float(result) + return result + + def find_operating_point( + self, params: TunnelDiodeParameters + ) -> tuple[float, float]: + if params.resistance <= 0: + raise ValueError("Сопротивление должно быть положительным") + + def residual(u: float) -> float: + return float( + self.iv_curve(u) + - (params.voltage_source - u) / params.resistance + ) + + u_low = 1.0e-6 + u_high = max(params.voltage_source, 0.1) + + try: + u_star = optimize.brentq(residual, u_low, u_high, xtol=1e-9) + except ValueError: + u_grid = np.linspace(u_low, u_high, 1000) + residuals = np.array([residual(u) for u in u_grid]) + sign_changes = np.where(np.diff(np.sign(residuals)) != 0)[0] + if len(sign_changes) == 0: + raise ValueError("Рабочая точка не найдена") from None + idx = int(sign_changes[0]) + u_star = optimize.brentq( + residual, u_grid[idx], u_grid[idx + 1], xtol=1e-9 + ) + + i_star = float(self.iv_curve(u_star)) + return float(u_star), i_star + + def equations_of_motion( + self, _t: float, state: np.ndarray, params: TunnelDiodeParameters + ) -> np.ndarray: + u_c, i_l = state + du_dt = (i_l - float(self.iv_curve(u_c))) / params.capacitance + di_dt = ( + params.voltage_source - params.resistance * i_l - u_c + ) / params.inductance + return np.array([du_dt, di_dt]) + + @st.cache_data(ttl=300) + def simulate( + _self, # noqa: N805 + params: TunnelDiodeParameters, + ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + initial_state = np.array([params.u_c0, params.i_l0]) + t_eval = np.linspace(0.0, params.t_max, params.n_points) + + solution = integrate.solve_ivp( + fun=lambda t, state: _self.equations_of_motion(t, state, params), + t_span=(0.0, params.t_max), + y0=initial_state, + t_eval=t_eval, + method="RK45", + rtol=1e-8, + atol=1e-10, + ) + + u_c = solution.y[0] + i_l = solution.y[1] + t = solution.t + return u_c, i_l, t + + def natural_frequency( # noqa: PLR6301 + self, params: TunnelDiodeParameters + ) -> float: + if params.inductance <= 0 or params.capacitance <= 0: + return 0.0 + return 1.0 / np.sqrt(params.inductance * params.capacitance) + + def quality_factor( # noqa: PLR6301 + self, params: TunnelDiodeParameters + ) -> float: + if params.resistance <= 0 or params.capacitance <= 0: + return float("inf") + return ( + 1.0 + / params.resistance + * np.sqrt(params.inductance / params.capacitance) + ) + + def critical_resistance(self, params: TunnelDiodeParameters) -> float: + u_star, _ = self.find_operating_point(params) + slope = float(self.iv_derivative(u_star)) + if slope >= 0: + return 0.0 + return -slope * params.inductance / params.capacitance + + def effective_damping(self, params: TunnelDiodeParameters) -> float: + u_star, _ = self.find_operating_point(params) + slope = float(self.iv_derivative(u_star)) + return params.resistance / (2 * params.inductance) + slope / ( + 2 * params.capacitance + ) diff --git a/src/models/tunnel_diode/sidebar.py b/src/models/tunnel_diode/sidebar.py new file mode 100644 index 0000000..c85d9a2 --- /dev/null +++ b/src/models/tunnel_diode/sidebar.py @@ -0,0 +1,97 @@ +import streamlit as st + +from models.tunnel_diode.objects import TunnelDiodeParameters + + +def create_tunnel_diode_sidebar() -> TunnelDiodeParameters: + st.sidebar.header("Источник и резистор") + + voltage_source = st.sidebar.number_input( + "ЭДС источника ℰ (В)", + min_value=0.1, + max_value=5.0, + value=1.5, + step=0.05, + format="%.3f", + ) + + resistance = st.sidebar.number_input( + "Сопротивление R (Ом)", + min_value=1.0, + max_value=1000.0, + value=50.0, + step=1.0, + format="%.1f", + ) + + st.sidebar.header("LC-контур") + + inductance_uh = st.sidebar.number_input( + "Индуктивность L (мкГн)", + min_value=1.0, + max_value=10000.0, + value=1000.0, + step=10.0, + format="%.1f", + help="1 мкГн = 10⁻⁶ Гн", + ) + + capacitance_nf = st.sidebar.number_input( + "Ёмкость C (нФ)", + min_value=0.1, + max_value=10000.0, + value=100.0, + step=1.0, + format="%.2f", + help="1 нФ = 10⁻⁹ Ф", + ) + + st.sidebar.header("Начальные условия") + + u_c0 = st.sidebar.number_input( + "Начальное напряжение U_C(0) (В)", + min_value=-2.0, + max_value=5.0, + value=0.0, + step=0.05, + format="%.3f", + ) + + i_l0_ma = st.sidebar.number_input( + "Начальный ток I_L(0) (мА)", + min_value=-100.0, + max_value=100.0, + value=0.0, + step=0.5, + format="%.3f", + ) + + st.sidebar.header("Параметры симуляции") + + t_max_us = st.sidebar.number_input( + "Время симуляции (мкс)", + min_value=10.0, + max_value=10000.0, + value=500.0, + step=10.0, + format="%.0f", + ) + + n_points = st.sidebar.slider( + "Точек траектории", + min_value=1000, + max_value=20000, + value=5000, + step=500, + ) + + return TunnelDiodeParameters( + voltage_source=voltage_source, + resistance=resistance, + inductance=inductance_uh * 1e-6, + capacitance=capacitance_nf * 1e-9, + u_c0=u_c0, + i_l0=i_l0_ma * 1e-3, + t_max=t_max_us * 1e-6, + n_points=n_points, + ) diff --git a/src/models/tunnel_diode/utils.py b/src/models/tunnel_diode/utils.py new file mode 100644 index 0000000..33cb351 --- /dev/null +++ b/src/models/tunnel_diode/utils.py @@ -0,0 +1,146 @@ +import numpy as np + +from models.tunnel_diode.objects import TunnelDiodeModel, TunnelDiodeParameters + +REGIME_DAMPED = "damped" +REGIME_LIMIT_CYCLE = "limit_cycle" +REGIME_RELAXATION = "relaxation" +REGIME_GROWING = "growing" + +THD_PURE_THRESHOLD = 0.05 +THD_RELAXATION_THRESHOLD = 0.30 +STEADY_TAIL_FRACTION = 0.5 +MIN_SIGNAL_LENGTH = 4 +MIN_ZERO_CROSSINGS = 2 +MIN_AMPLITUDE_THRESHOLD = 1e-3 +SMALL_AMPLITUDE_DAMPED = 0.05 +DAMPED_DECAY_RATIO = 0.5 +GROWTH_RATIO = 2.0 +HARMONIC_RANGE = 6 +MIN_FUNDAMENTAL_INDEX = 1 +MIN_SPECTRUM_LENGTH = 3 + + +def oscillation_metrics(u_c: np.ndarray, t: np.ndarray) -> dict: + if len(u_c) < MIN_SIGNAL_LENGTH: + return {"error": "Слишком короткий сигнал"} + + tail_start = int(len(u_c) * STEADY_TAIL_FRACTION) + tail = u_c[tail_start:] + t_tail = t[tail_start:] + + u_mean = float(np.mean(tail)) + u_max = float(np.max(tail)) + u_min = float(np.min(tail)) + amplitude = (u_max - u_min) / 2 + + period = 0.0 + centered = tail - u_mean + zero_crossings = [] + for i in range(1, len(centered)): + if centered[i - 1] <= 0 < centered[i]: + dt_lin = ( + -centered[i - 1] + * (t_tail[i] - t_tail[i - 1]) + / (centered[i] - centered[i - 1]) + ) + zero_crossings.append(t_tail[i - 1] + dt_lin) + + if len(zero_crossings) >= MIN_ZERO_CROSSINGS: + intervals = np.diff(zero_crossings) + period = float(np.mean(intervals)) + + head = u_c[: max(1, tail_start)] + head_amplitude = ( + float(np.max(head) - np.min(head)) / 2 if len(head) > 1 else 0.0 + ) + + return { + "amplitude": amplitude, + "u_mean": u_mean, + "u_max": u_max, + "u_min": u_min, + "period": period, + "head_amplitude": head_amplitude, + } + + +def compute_spectrum( + u_c: np.ndarray, t: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + if len(u_c) < MIN_SIGNAL_LENGTH or len(t) < MIN_SIGNAL_LENGTH: + return np.array([]), np.array([]) + + tail_start = int(len(u_c) * STEADY_TAIL_FRACTION) + signal = u_c[tail_start:] - np.mean(u_c[tail_start:]) + dt = float(t[1] - t[0]) if len(t) > 1 else 1.0 + + spectrum = np.abs(np.fft.rfft(signal)) + freqs = np.fft.rfftfreq(len(signal), d=dt) + return freqs, spectrum + + +def total_harmonic_distortion(freqs: np.ndarray, spectrum: np.ndarray) -> float: + if len(spectrum) < MIN_SPECTRUM_LENGTH: + return 0.0 + + spectrum_no_dc = spectrum.copy() + spectrum_no_dc[0] = 0.0 + fund_idx = int(np.argmax(spectrum_no_dc)) + if fund_idx < MIN_FUNDAMENTAL_INDEX: + return 0.0 + + fundamental = float(spectrum_no_dc[fund_idx]) + if fundamental <= 0: + return 0.0 + + harmonics_power = 0.0 + for k in range(2, HARMONIC_RANGE): + harmonic_idx = fund_idx * k + if harmonic_idx >= len(spectrum_no_dc): + break + harmonics_power += float(spectrum_no_dc[harmonic_idx]) ** 2 + + return float(np.sqrt(harmonics_power) / fundamental) + + +def _is_damped(tail_amp: float, head_amp: float) -> bool: + if tail_amp < MIN_AMPLITUDE_THRESHOLD: + return True + return ( + tail_amp < head_amp * DAMPED_DECAY_RATIO + and tail_amp < SMALL_AMPLITUDE_DAMPED + ) + + +def classify_regime( + model: TunnelDiodeModel, + params: TunnelDiodeParameters, + u_c: np.ndarray, + t: np.ndarray, +) -> str: + metrics = oscillation_metrics(u_c, t) + if "error" in metrics: + return REGIME_DAMPED + + head_amp = metrics["head_amplitude"] + tail_amp = metrics["amplitude"] + + if _is_damped(tail_amp, head_amp): + return REGIME_DAMPED + + if tail_amp > head_amp * GROWTH_RATIO: + return REGIME_GROWING + + freqs, spectrum = compute_spectrum(u_c, t) + thd = total_harmonic_distortion(freqs, spectrum) + + if thd > THD_RELAXATION_THRESHOLD: + return REGIME_RELAXATION + return REGIME_LIMIT_CYCLE + + +def is_self_oscillating( + model: TunnelDiodeModel, params: TunnelDiodeParameters +) -> bool: + return model.effective_damping(params) < 0 diff --git a/tests/m5_tunnel_diode/__init__.py b/tests/m5_tunnel_diode/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/m5_tunnel_diode/test_tunnel_diode_physics.py b/tests/m5_tunnel_diode/test_tunnel_diode_physics.py new file mode 100644 index 0000000..25a22f4 --- /dev/null +++ b/tests/m5_tunnel_diode/test_tunnel_diode_physics.py @@ -0,0 +1,290 @@ +import numpy as np +import pytest + +from models.tunnel_diode.config import DIODE_COEFFICIENTS +from models.tunnel_diode.objects import TunnelDiodeModel, TunnelDiodeParameters +from models.tunnel_diode.utils import ( + REGIME_DAMPED, + REGIME_LIMIT_CYCLE, + classify_regime, + compute_spectrum, + is_self_oscillating, + oscillation_metrics, + total_harmonic_distortion, +) + + +@pytest.fixture +def model() -> TunnelDiodeModel: + return TunnelDiodeModel() + + +@pytest.fixture +def damped_params() -> TunnelDiodeParameters: + return TunnelDiodeParameters( + voltage_source=1.5, + resistance=500.0, + inductance=1.0e-3, + capacitance=1.0e-7, + u_c0=1.5, + i_l0=0.0, + t_max=1.0e-3, + n_points=4000, + ) + + +@pytest.fixture +def oscillating_params() -> TunnelDiodeParameters: + return TunnelDiodeParameters( + voltage_source=1.5, + resistance=20.0, + inductance=1.0e-3, + capacitance=1.0e-7, + u_c0=1.5, + i_l0=0.0, + t_max=2.0e-3, + n_points=10000, + ) + + +def test_iv_curve_zero_at_zero(model: TunnelDiodeModel) -> None: + assert float(model.iv_curve(0.0)) == pytest.approx(0.0) + + +def test_iv_curve_zero_for_negative(model: TunnelDiodeModel) -> None: + assert float(model.iv_curve(-0.5)) == pytest.approx(0.0) + assert float(model.iv_curve(-1.0)) == pytest.approx(0.0) + + +def test_iv_curve_matches_polynomial(model: TunnelDiodeModel) -> None: + a = DIODE_COEFFICIENTS["a"] + b = DIODE_COEFFICIENTS["b"] + c = DIODE_COEFFICIENTS["c"] + u = 0.5 + expected = a * u**3 + b * u**2 + c * u + assert float(model.iv_curve(u)) == pytest.approx(expected, rel=1e-12) + + +def test_iv_curve_vectorized(model: TunnelDiodeModel) -> None: + u_vec = np.array([-1.0, 0.0, 0.5, 1.0, 2.0]) + result = np.asarray(model.iv_curve(u_vec)) + assert result.shape == u_vec.shape + assert result[0] == 0.0 + assert result[1] == 0.0 + assert result[2] > 0 + + +def test_iv_derivative_consistency(model: TunnelDiodeModel) -> None: + u = 1.5 + h = 1e-5 + numerical = ( + float(model.iv_curve(u + h)) - float(model.iv_curve(u - h)) + ) / (2 * h) + analytical = float(model.iv_derivative(u)) + assert numerical == pytest.approx(analytical, rel=1e-4) + + +def test_iv_has_negative_resistance_region(model: TunnelDiodeModel) -> None: + u_grid = np.linspace(0.5, 2.0, 50) + derivatives = np.asarray(model.iv_derivative(u_grid)) + assert np.any(derivatives < 0) + + +def test_operating_point_satisfies_kirchhoff( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + u_star, i_star = model.find_operating_point(damped_params) + expected_i = ( + damped_params.voltage_source - u_star + ) / damped_params.resistance + assert i_star == pytest.approx(expected_i, rel=1e-6) + assert i_star == pytest.approx(float(model.iv_curve(u_star)), rel=1e-6) + + +def test_natural_frequency_matches_lc( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + expected = 1.0 / np.sqrt( + damped_params.inductance * damped_params.capacitance + ) + assert model.natural_frequency(damped_params) == pytest.approx(expected) + + +def test_quality_factor_formula( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + expected = ( + np.sqrt(damped_params.inductance / damped_params.capacitance) + / damped_params.resistance + ) + assert model.quality_factor(damped_params) == pytest.approx(expected) + + +def test_critical_resistance_formula( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + u_star, _ = model.find_operating_point(damped_params) + slope = float(model.iv_derivative(u_star)) + if slope < 0: + expected = -slope * damped_params.inductance / damped_params.capacitance + assert model.critical_resistance(damped_params) == pytest.approx( + expected, rel=1e-9 + ) + + +def test_simulation_initial_conditions( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + u_c, i_l, t = model.simulate(damped_params) + assert u_c[0] == pytest.approx(damped_params.u_c0) + assert i_l[0] == pytest.approx(damped_params.i_l0) + assert t[0] == pytest.approx(0.0) + + +def test_simulation_returns_correct_length( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + u_c, i_l, t = model.simulate(damped_params) + assert len(u_c) == damped_params.n_points + assert len(i_l) == damped_params.n_points + assert len(t) == damped_params.n_points + + +def test_damped_regime_decays( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + u_c, _, _ = model.simulate(damped_params) + head = u_c[: len(u_c) // 4] + tail = u_c[3 * len(u_c) // 4 :] + head_amplitude = float(np.max(head) - np.min(head)) + tail_amplitude = float(np.max(tail) - np.min(tail)) + assert tail_amplitude < head_amplitude + + +def test_damped_classified_as_damped( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + u_c, _, t = model.simulate(damped_params) + regime = classify_regime(model, damped_params, u_c, t) + assert regime == REGIME_DAMPED + + +def test_oscillating_regime_has_finite_amplitude( + model: TunnelDiodeModel, oscillating_params: TunnelDiodeParameters +) -> None: + u_c, _, t = model.simulate(oscillating_params) + metrics = oscillation_metrics(u_c, t) + assert metrics["amplitude"] > 0.01 + + +def test_oscillating_classified_not_damped( + model: TunnelDiodeModel, oscillating_params: TunnelDiodeParameters +) -> None: + u_c, _, t = model.simulate(oscillating_params) + regime = classify_regime(model, oscillating_params, u_c, t) + assert regime != REGIME_DAMPED + + +def test_self_oscillation_consistent_with_damping( + model: TunnelDiodeModel, + damped_params: TunnelDiodeParameters, + oscillating_params: TunnelDiodeParameters, +) -> None: + assert is_self_oscillating(model, oscillating_params) + assert not is_self_oscillating(model, damped_params) + + +def test_oscillation_period_close_to_lc( + model: TunnelDiodeModel, oscillating_params: TunnelDiodeParameters +) -> None: + u_c, _, t = model.simulate(oscillating_params) + metrics = oscillation_metrics(u_c, t) + expected_period = ( + 2 + * np.pi + * np.sqrt( + oscillating_params.inductance * oscillating_params.capacitance + ) + ) + assert metrics["period"] == pytest.approx(expected_period, rel=0.30) + + +def test_compute_spectrum_returns_real_values( + model: TunnelDiodeModel, oscillating_params: TunnelDiodeParameters +) -> None: + u_c, _, t = model.simulate(oscillating_params) + freqs, spectrum = compute_spectrum(u_c, t) + assert len(freqs) == len(spectrum) + assert len(freqs) > 0 + assert np.all(spectrum >= 0) + + +def test_compute_spectrum_handles_short_signal() -> None: + freqs, spectrum = compute_spectrum( + np.array([0.1, 0.2]), np.array([0.0, 1.0]) + ) + assert len(freqs) == 0 + assert len(spectrum) == 0 + + +def test_thd_zero_for_pure_sine() -> None: + n = 4096 + t = np.linspace(0.0, 1.0, n) + signal = np.sin(2 * np.pi * 50 * t) + freqs, spectrum = compute_spectrum(signal, t) + thd = total_harmonic_distortion(freqs, spectrum) + assert thd < 0.05 + + +def test_thd_high_for_square_wave() -> None: + n = 4096 + t = np.linspace(0.0, 1.0, n) + signal = np.sign(np.sin(2 * np.pi * 50 * t)) + freqs, spectrum = compute_spectrum(signal, t) + thd = total_harmonic_distortion(freqs, spectrum) + assert thd > 0.10 + + +def test_pure_sine_classified_as_limit_cycle() -> None: + n = 8000 + t = np.linspace(0.0, 4.0e-3, n) + u_c = 1.0 + 0.1 * np.sin(2 * np.pi * 1.0e4 * t) + + model = TunnelDiodeModel() + params = TunnelDiodeParameters( + voltage_source=1.5, + resistance=20.0, + inductance=1.0e-3, + capacitance=1.0e-7, + u_c0=1.0, + t_max=4.0e-3, + n_points=n, + ) + regime = classify_regime(model, params, u_c, t) + assert regime == REGIME_LIMIT_CYCLE + + +def test_zero_voltage_zero_current(model: TunnelDiodeModel) -> None: + params = TunnelDiodeParameters( + voltage_source=0.0, + resistance=100.0, + inductance=1.0e-3, + capacitance=1.0e-7, + u_c0=0.0, + i_l0=0.0, + t_max=1.0e-4, + n_points=1000, + ) + u_c, i_l, _ = model.simulate(params) + assert np.all(np.abs(u_c) < 1e-9) + assert np.all(np.abs(i_l) < 1e-9) + + +def test_equations_of_motion_at_operating_point_zero_derivative( + model: TunnelDiodeModel, damped_params: TunnelDiodeParameters +) -> None: + u_star, i_star = model.find_operating_point(damped_params) + state = np.array([u_star, i_star]) + derivatives = model.equations_of_motion(0.0, state, damped_params) + assert abs(derivatives[0]) < 1e-9 + assert abs(derivatives[1]) < 1e-9