Заданная система уравнений для скалярных комплексных амплитуд должна быть выполнена в любой точке границы раздела сред:
$$\forall y \quad E_y(x = -0) = E_y(x = +0);\\
\begin{cases}
E_y(x = -0) =f \cdot e^{-ik_{1y}y}, \; f\neq f(y);\\
E_y(x = +0) =g \cdot e^{-ik_{2y}y}, \; g\neq g(y).
\end{cases}$$Отсюда с необходимостью следует $k_{1y} = k_{2y}.$ Аналогично, рассматривая вторую границу раздела сред, получаем $k_{2y} = k_{3y}.$

Используя выражение для $k$ и равенство $k_y=k_{1y}=k_{2y}=k_{3y}$:
$$k=\sqrt{k^2_{1x}+k^2_{y}}\Rightarrow k_{1x}=\sqrt{\cfrac{\varepsilon_1 \omega^2}{c^2}-k_y^2}=\cfrac{\sqrt{\varepsilon_1}\omega}{c}\cos{\theta_1}.$$Аналогично,
$$k_{2x}=\cfrac{\sqrt{\varepsilon_2}\omega}{c}\cos{\theta_2}=\cfrac{\sqrt{\varepsilon_2}\omega}{c}\sqrt{1-\cfrac{\varepsilon_1}{\varepsilon_2}\sin^2{\theta_1}},$$$$k_{3x}= \cfrac{\sqrt{\varepsilon_3}\omega}{c}\cos{\theta_3}=\cfrac{\sqrt{\varepsilon_3}\omega}{c}\sqrt{1-\cfrac{\varepsilon_1}{\varepsilon_3}\sin^2{\theta_1}}.$$

Запишем граничные условия на границе 1-2 для тангенциальных компонент $E$ и $B$ соответственно, считая $x=0$:
$$\mathcal{E}+r\mathcal{E}=a\mathcal{E}+b\mathcal{E},$$$$k_{1x}\mathcal{E}-k_{1x}\mathcal{E}r=k_{2x}a\mathcal{E}-k_{2x}b\mathcal{E} .$$Запишем граничные условия на границе 2-3 для тангенциальных компонент $E$ и $B$ соответственно:
$$a\mathcal{E}e^{ik_{2x}d}+b\mathcal{E}e^{-ik_{2x}d}=t\mathcal{E}\ ,$$$$k_{3x}t\mathcal{E}=k_{2x}\mathcal{E} ae^{ik_{2x}d}-k_{2x}\mathcal{E}be^{-ik_{2x}d} .$$Таким образом получаем систему уравнений:
\[
\begin{cases}
1+r=a+b,\\
k_{1x}(1-r)=k_{2x}(a-b), \\
ae^{i\phi}+be^{-i\phi}=t,\\
k_{3x}t=k_{2x}(ae^{i\phi}-be^{-i\phi}),
\end{cases}
\]
где введено обозначение $\phi = k_{2x}d.$

Из записанной в предыдущем пункте системы получаем:
$$\begin{cases}
t = (a e^{i\phi} - b e^{-i\phi})\dfrac{k_{2x}}{k_{3x}},\\
a = b e^{-2i\phi} \cdot \dfrac{k_{2x} + k_{3x}}{k_{2x} - k_{3x}},\\
a(k_{1x} + k_{2x}) + b(k_{1x} - k_{2x}) = 2 k_{1x},\\
r = \dfrac{a(k_{1x} - k_{2x}) + b(k_{1x} + k_{2x})}{2k_{1x}}.
\end{cases}$$Решая, приходим к результату:
$$r_s = \dfrac{(k_{1x} - k_{2x}) (k_{2x}+ k_{3x}) +(k_{1x} + k_{2x}) (k_{2x}- k_{3x})e^{2i\phi} }{(k_{1x} - k_{2x}) (k_{2x}- k_{3x})e^{2i\phi} + (k_{1x} +k_{2x}) (k_{2x}+ k_{3x})}.$$Вводя обозначения
$$\begin{cases}
A = \dfrac{k_{2x} + k_{3x}}{k_{2x} - k_{3x}},\\
B = \dfrac{k_{1x} + k_{2x}}{k_{1x} - k_{2x}}.
\end{cases}$$получаем выражение, приведенное к требуемому виду.

Аналогично А3 записываем условия на двух границах раздела с учетом связи электрического и магнитного полей в ЭМ-волне: $\sqrt{\varepsilon\varepsilon_0\mu_0}\mathcal{E} = \mathcal{B}$:

\[
\begin{cases}
1+r=a+b,\\
\dfrac{k_{1x}}{\varepsilon_1}(1-r)=\dfrac{k_{2x}}{\varepsilon_2}(a-b), \\
ae^{ik_{2x}d}+be^{-ik_{2x}d}=t,\\
\dfrac{k_{3x}}{\varepsilon_3}t=\dfrac{k_{2x}}{\varepsilon_2}(ae^{ik_{2x}d}-be^{-ik_{2x}d}).
\end{cases}
\]
Или, с учетом предложенных в условии переобозначений:
\[
\begin{cases}
1+r=a+b,\\
\kappa_{1x}(1-r)=\kappa_{2x}(a-b), \\
ae^{i\phi}+be^{-i\phi}=t,\\
\kappa_{3x}t=\kappa_{2x}(ae^{i\phi}-be^{-i\phi}).
\end{cases}
\]

Заметим, что система в пункте A5 совпадает с системой из пункта А3 с точностью до замены $\kappa_{ix} \leftrightarrow k_{ix}$. Соответственно, сразу записываем ответ:
$$r_p = \frac{A +B e^{2i \phi}}{AB + e^{2\phi}},$$ $$\begin{cases}
A = \dfrac{\kappa_{2x} + \kappa_{3x}}{\kappa_{2x} - \kappa_{3x}},\\
B = \dfrac{\kappa_{1x} + \kappa_{2x}}{\kappa_{1x} - \kappa_{2x}}.
\end{cases}$$

Будем получать графики для трех толщин: 10, 20 и 30 нм. При малых углах падения видно, что диапазон изменения $\Psi$ составляет сотые градуса. Максимальная чувствительность достигается при \[\theta_\text{opt}=77^\circ.\] Соответствующие графики для сопоставления приведены ниже

Опишем алгоритм поиска оптимальных $\Phi$ и $d$. Для начала подадим на вход программе произвольные $\Phi$ и $d$ in B1.py. Далее, на каждом шаге будем изменять немного $\Phi$ или $d$ и следить за "error", выдаваемой программой B1.py. Если "error" уменьшается, продолжим изменять $\Phi$ или $d$ в том же направлении (например, если мы на очередном шаге уменьшали $\Phi$ на $1^\circ$, и при этом "error" уменьшилась, то на следующем шаге тоже уменьшим $\Phi$ на $1^\circ$). В конце концов мы достигнем $\Phi$ и $d$ таких, что любое их изменение приводит к увеличению "error", то есть мы достигли минимально возможной "error". Это означает, что мы нашли оптимальные $\Phi$ и $d$.

На графике можно увидеть, как подбираемые данные соотносятся с данными коллеги, но численным показателем качества аппроксимации является параметр "error", выдаваемый программой B1.py.

Алгоритм поиска оптимальных $\Phi$ и $d$ полностью совпадает с описанным в пункте $\textbf{B1}$.

Зададим $\text{name} = \text{b3}$, $\text{Phi} = 77$, $T = 1000$, $t_\max = 120$ в программе B_in.txt, затем запустим gen_data.py и B.py последовательно. В папке "res" мы получим файл с требуемым графиком. Заметим, что мы подставляли $\Phi = \theta_\text{opt}$ из пункта $\textbf{A7}$, чтобы увеличить качество эксперимента.

Рассмотрим малый промежуток времени $\Delta t$ и участок границы $\rm SiO_2$-$\rm Si$ площади $S$. Тогда масса окисленного $\rm Si$ за промежуток $\Delta t$ равна $\Delta m = F_3 \, \Delta t \, S \, \mu_{\rm SiO_2}$. С другой стороны, граница $\rm SiO_2$-$\rm Si$ сдвинется на $v \, \Delta t$, так что возникшая масса $\rm SiO_2$ равна $\Delta m = \rho_{\rm SiO_2} \, S \, v \, \Delta t$. Поэтому $\rho_{\rm SiO_2} \, v \, S\, \Delta t = F_3 \, \Delta t \, S \, \mu_{\rm SiO_2}$, так что $\rho_{\rm SiO_2} \, v = F_3 \, \mu_{\rm SiO_2}$.

Стационарные условия дают $F_1 = F_2 = F_3$, i.e. $h (C^\ast - C_0) = \frac{D}{d} (C_0 - C_i) = k C_i$. Из второго равенства мы получаем $C_0 = (1 + \frac{kd}{D}) C_i$. Подставим это в равенство $h (C^\ast - C_0) = k C_i$ и получим $h (C^\ast - (1 + \frac{kd}{D}) C_i) = k C_i$, что влечет $C^\ast = (\frac{k}{h} + \frac{kd}{D} + 1) C_i$, так что
\[C_i = \frac{C^\ast}{1 + \frac{k}{h} + \frac{kd}{D}}, \quad C_0 = \left (1 + \frac{kd}{D} \right ) C_i = C^\ast \frac{1 + \frac{kd}{D}}{1 + \frac{k}{h} + \frac{kd}{D}}\]

Из пункта $\textbf{B4}$ мы имеем $v = \dot d = \frac{F_3 \mu_{\rm SiO_2}}{\rho_{\rm SiO_2}} = \frac{k C_i \mu_{\rm SiO_2}}{\rho_{\rm SiO_2}}$. Учитывая выражение для $C_i$, полученное в пункте $\textbf{B5}$, мы получаем
\[ \dot{d} = C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}} \frac{1}{\frac{1}{k} + \frac{1}{h} + \frac{d}{D}} \]или, эквивалентно,
\[\frac{\dot d \cdot d }{D} + \dot d \cdot \left ( \frac{1}{k} + \frac{1}{h} \right) - C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}} = 0\]Заметим, что левая часть равенства равна
\[ \frac{d}{dt} \left ( \frac{d^2}{2D} + d \cdot \left ( \frac{1}{k} + \frac{1}{h} \right) - t \cdot C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}} \right ) = 0\]Поэтому
\[\frac{d^2}{2D} + d \cdot \left ( \frac{1}{k} + \frac{1}{h} \right) - t \cdot C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}} = \text{const}\]где $\text{const}$ может быть найдена из начальных условий $d(t_0) = d_0$. Именно, \[\text{const} = \frac{d_0^2}{2D} + d_0 \cdot \left ( \frac{1}{k} + \frac{1}{h} \right) - t_0 \cdot C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}}\]Тем самым мы получаем
\[d^2 + 2D \cdot \left ( \frac{1}{k} + \frac{1}{h} \right) \cdot d = 2D C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}} ( t + \tau)\]где $\tau = \frac{\text{const} \cdot \rho_{\rm Si O_2}}{2D C^* \mu_{\rm SiO_2}}$.

Для начала создадим необходимые для программы B7.py данные. Запишем $T$, $\Phi = \theta_\text{opt} = 77^\circ$, $t_\max = 120$ в файл B_in.txt, затем запустим gen_data.py и B.py, и наконец запустим B7.py пять раз, на каждом шаге вводя "name", $T$ и $\Phi$ из файла B_in.txt. Запишем каждый вывод программы B7.py в файл B8_in.txt и получим таблицу, приведенную ниже.

Отметим, что даже для $T = 850^\circ C$ выводится ошибка "Final layer is too thin to approximate", а для $T = 900^\circ C$ получается слишком большая погрешность, так что разумно начать с $T = 950^\circ C$ с шагом $50^\circ C$.

Для начала построим графики $\ln B \left (\frac{1}{T} \right)$ и $\ln \frac{B}{A} \left (\frac{1}{T} \right)$, используя программу B8.py. В консоли мы можем увидеть коэффициенты наклона предложенных линейных зависимостей. Именно, коэффициент наклона графика $\ln B \left (\frac{1}{T} \right)$ равен $-(14671 \pm 625) \: {}^\circ C$, а коэффициент наклона графика $\ln \frac{B}{A} \left (\frac{1}{T} \right)$ равен $-(23418 \pm 1730) \: {}^\circ C$.

Из части $\text{B6}$ мы получаем
\[B = 2D C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}} \propto e^{-\frac{\Delta E_a}{k_B T}}\]Поэтому коэффициент наклона графика $\ln B \left (\frac{1}{T} \right)$ равен $-\frac{\Delta E_a}{k_B}$. Значит, $\Delta E_a = (2.208 \pm 0.086) \cdot 10^{-19} \: \text{J}$, или, эквивалентно, \[\Delta E_a = \frac{(2.208 \pm 0.086) \cdot 10^{-19}}{1.60 \cdot 10^{-19}} \: eV = (1.38 \pm 0.05) \: eV\]Далее, неравенство $h \gg k$ влечет $A \approx \frac{2D}{k}$, так что
\[\frac{B}{A} \approx k C^* \frac{\mu_{\rm SiO_2}}{\rho_{\rm Si O_2}} \propto e^{-\frac{\Delta E_b}{k_B T}}\]коэффициент наклона графика $\ln \frac{B}{A} \left (\frac{1}{T} \right)$ равен $-\frac{\Delta E_b}{k_B}$. Поэтому $\Delta E_b = (3.23 \pm 0.23) \cdot 10^{-19} \: \text{J}$, или, эквивалентно, \[\Delta E_b = \frac{(3.23 \pm 0.23) \cdot 10^{-19}}{1.60 \cdot 10^{-19}} \: eV = (2.02 \pm 0.14) \: eV\]