diff --git a/Interpolation/inter.tex b/Interpolation/inter.tex
new file mode 100644
index 0000000..9e7a0b5
--- /dev/null
+++ b/Interpolation/inter.tex
@@ -0,0 +1,217 @@
+\documentclass[12pt]{article}
+\usepackage{graphicx}
+\usepackage{color}
+\usepackage[utf8]{inputenc}
+\usepackage[T2A]{fontenc}
+\usepackage[russian]{babel}
+\usepackage{amsmath}
+\usepackage{amssymb}
+
+\topmargin=-1.8cm
+\oddsidemargin=-15mm
+\evensidemargin=-15mm
+\textheight=24.5cm
+\textwidth=19cm
+\tolerance=1000
+
+\parskip=5pt plus 4pt minus 2pt
+\tolerance=9000
+%
+\title{Интерполяция полей скорости}
+\begin{document}
+\maketitle
+\section{Координаты и уравнения движения}
+Пусть у нас задано поле скорости $(u_{i,j},v_{i,j})$ на некоторой двумерной сетке $(\lambda_{i,j},\varphi_{i,j})$, $i=0,1,\dots,N_x$, $j=0,1,\dots,N_y$, $\lambda_{i,j}$~--- долгота, $\varphi_{i,j}$~--- широта. Мы предполагаем, что эта сетка может быть преобразована к прямоугольному виду, то есть существуют такие функции $X(\lambda,\varphi)$ и $Y(\lambda,\varphi)$, что $x_{i,j}=X(\lambda_{i,j},\varphi_{i,j})=i\Delta x$ и $y_{i,j}=Y(\lambda_{i,j},\varphi_{i,j})=j\Delta y$, где $\Delta x$, $\Delta y$~--- некоторые постоянные числа. Проще говоря, существуют координаты, в которых наша сетка имеет простой прямоугольный вид. В самом простом случае сетка $(\lambda_{i,j},\varphi_{i,j})$ уже имеет прямоугольный вид. Но, для глобальных полей такую сетку не используют в силу сингулярности на Северном полюсе и зависимости разрешения от широты. Сами скорости задаются обычно в см/сек. Однако, для описания крупномасштабных процессов это неудобно и лучше использовать скорости, выраженные в км/сут
+%
+\begin{equation}
+1\,\frac{\text{см}}{\text{сек}}=\frac{10^{-5}\,\text{км}}{\frac{1}{24*3600}\,\text{сут}}=0{,}864\frac{\text{км}}{\text{сут}}.
+\end{equation}
+%
+
+Есть несколько возможных подходов к написанию лагранжевых уравнений адвекции.
+%
+\begin{enumerate}
+\item Не преобразовывать ни скорости, ни координаты. Уравнения адвекции выглядят следующим образом:
+%
+\begin{equation}
+\dot\lambda=\frac{10800\cos{\varphi}}{\pi R}u,\quad \dot\varphi=\frac{10800}{\pi R}v,
+\end{equation}
+%
+где $R$~--- радиус Земли. Широта и долготы выражены в минутах, время~--- в сутках, $u$ и $v$~--- в км/сут. Коэффициент $\frac{\pi R}{10800}$ примерно равен морской миле~--- $1{,}852$ км.\\
+{\em Плюсы:} нет необходимости в каких-либо преобразованиях координат, входные и выходные данные задаются непосредственно в географических координатах.\\
+{\em Минусы:} вычисление косинуса на каждом шаге~--- затратная операция, в случае непрямоугольной сетки определение ячейки, в которой находятся текущие координаты (необходимо для нахождения интерполированых скоростей), является непростой операцией.
+\item Перейдя от линейных скоростей в правой части к угловым, мы получаем максимально простые уравнения адвекции
+%
+\begin{equation}
+\dot\lambda=u,\quad \dot\varphi=v.
+\label{maineq}
+\end{equation}
+%
+Здесь $u$ и $v$ выражены в мин/сут.\\
+{\em Плюсы:} те же, что и у первого способа, быстрое численное интегрирование уравнений.\\
+{\em Минусы:} всё те же проблемы с непрямоугольными сетками.
+\item Преобразовать координаты и скорости в прямоугольную сетку.
+Уравнения движения будут иметь тот же вид, что и в предыдущем случае, но для некоторых абстрактных координат и скоростей.
+%
+\begin{equation}
+\dot x=u,\quad \dot y=v.
+\end{equation}
+%
+{\em Плюсы:} максимально простое и быстрое численное интегрирование.\\
+{\em Минусы:} мы считаем в абстрактных координатах, следовательно, все входные и выходные данные должны преобразовываться. Кроме этого, при данном подходе неудобно рассчитывать различные интегральные функции от траектории (например, её длину) из-за необходимости преобразования метрики на каждом шаге, что может сильно увеличить время счёта.
+\end{enumerate}
+%
+Мы считаем последним способом ещё и потому, что в этом случае програмный код получается универсальным и не зависит от реальной сетки $(\lambda_{ij},\varphi_{ij})$.
+
+\section{Интерполяция}
+Поле скорости у нас задано только на сетке $(x_{i,j},y_{i,j},t_k)$, однако, для интегрирования уравнений адвекции необходимы значения скорости в произвольной точке пространства $(x,y)$ и в произвольный момент времени $t$. Для нахождения этих значений используются различные методы интерполяции. Мы используем два: линейную и бикубическую интерполяцию в пространстве и линейную и кубическими полиномами третьего порядка во времени.
+\subsection{Интерполяция в пространстве}
+Рассмотрим отдельную ячейку нашей пространственной сетки единичной ширины и высоты (это всегда можно сделать соответствующей нормировкой). В четырёх углах углах ячейки нам известны значения некоторой функции $z(x,y)$, а именно $z(0,0)=z_{00}$, $z(1,0)=z_{10}$, $z(0,1)=z_{01}$, $z(1,1)=z_{11}$ (см. Рис.~\ref{Fig-cell})
+%
+\begin{figure}[!htb]
+\begin{center}
+\includegraphics[width=0.4\textwidth]{net.eps}
+\caption{Элементарная ячейка сетки, в которой производится интерполяция.}
+\label{Fig-cell}
+\end{center}
+\end{figure}
+%
+Рассмотрим сперва линейную интерполяцию. В этом случае мы считаем, что функция $z(x,y)$ является линейной вдоль любой вертикальной или горизонтальной линии. Соответственно, мы можем найти $z(x,y)$, найдя сперва $z(x,0)$, $z(x,1)$ (или $z(0,y)$, $z(1,y)$, результат будет тот же)
+%
+\begin{equation}
+\begin{gathered}
+z(x,0)=(z(1,0)-z(0,0))x+z(0,0)=(z_{10}-z_{00})x+z_{00},\\
+z(x,1)=(z(1,1)-z(0,1))x+z(0,1)=(z_{11}-z_{01})x+z_{01},\\
+z(x,y)=(z(x,1)-z(x,0))y+z(x,0)=axy+bx+cy+d,\\
+a=z_{11}-z_{01}-z_{10}+z_{00},\quad b=z_{10}-z_{00},\quad c=z_{01}-z_{00},\quad d=z_{00}.
+\end{gathered}
+\label{linint}
+\end{equation}
+%
+Таким образом, несмотря на название метода, интерполирующая функция является нелинейной. Линейная интерполяция обеспечивает непрерывность функции $z(x,y)$ при пересечении границ ячейки, но её производные разрывны на этих границах.
+
+Интерполирующая функция для бикубической интерполяции имеет вид
+%
+\begin{multline}
+z(x,y)=a_{00}+a_{10}x+a_{20}x^2+a_{30}x^3+a_{01}y+a_{11}xy+a_{21}x^2y+a_{31}x^3y+\\
++a_{02}y^2+a_{12}xy^2+a_{22}x^2y^2+a_{32}x^3y^2+a_{03}y^3+a_{13}xy^3+a_{23}x^2y^3+a_{33}x^3y^3
+\label{bicubic}
+\end{multline}
+%
+и содержит $16$ коэффициентов. Для их нахождения необходимо знать не только значения функции на сетке, но и значения производных $\partial z/\partial x$, $\partial z/\partial y$ и $\partial^2 z/\partial x \partial y$. Так как мы имеем только сами значения функции, нам необходимо найти эти производные путём численного дифференцирования
+%
+\begin{equation}
+\begin{aligned}
+\frac{\partial z}{\partial x}&=z^x_{i,j}=\frac{z_{i+1,j}-z_{i-1,j}}{2},\\
+\frac{\partial z}{\partial y}&=z^y_{i,j}=\frac{z_{i,j+1}-z_{i,j-1}}{2},\\
+\frac{\partial^2 z}{\partial x\partial y}&=z^{xy}_{i,j}=\frac{z^y_{i+1,j}-z^y_{i-1,j}}{2}.
+\end{aligned}
+\end{equation}
+%
+В точках, где использование центральных разностей невозможно (вблизи границ сетки или около берега), используются правые или левые разности. Таким образом, для каждой ячейки сетки мы получаем $16$  значений функции, её первых производных и смешаной второй производной в углах ячейки. Коэффициенты интерполирующей функции получаются следующим образом
+%
+\begin{equation}
+\begin{gathered}
+a_{00}=z_{00},\quad
+a_{01}=z^y_{00},\quad
+a_{02}=3(z_{01}-z_{00})-2z^y_{00}-z^y_{01},\quad
+a_{03}=2(z_{00}-z_{01})+z^y_{00}+z^y_{01},\\
+a_{10}=z^x_{00},\quad
+a_{11}=z^{xy}_{00},\quad
+a_{12}=3(z^x_{01}-z^x_{00})-2z^{xy}_{00}-z^{xy}_{01},\quad
+a_{13}=2(z^x_{00}-z^x_{01})+z^{xy}_{00}+z^{xy}_{01},\\
+a_{20}=3(z_{10}-z_{00})-2z^x_{00}-z^x_{10},\quad
+a_{21}=3(z^y_{10}-z^y_{00})-2z^{xy}_{00}-z^{xy}_{10},\\
+\begin{aligned}
+a_{22}=9(z_{00}-z_{01}-z_{10}+z_{11})&+6(z^x_{00}-z^x_{01}+z^y_{00}-z^y_{10})+{}\\
+{}&+3(z^x_{10}-z^x_{11}+z^y_{01}-z^y_{11})+2(z^{xy}_{01}+z^{xy}_{10})+4z^{xy}_{00}+z^{xy}_{11},
+\end{aligned}\\
+\begin{aligned}
+a_{23}=6(z_{01}+z_{10}-z_{00}-z_{11})+4(z^x_{01}-z^x_{00})&+3(z^y_{10}+z^y_{11}-z^y_{00}-z^y_{01})+{}\\
+{}&+2(z^x_{11}-z^x_{10}-z^{xy}_{00}-z^{xy}_{01})-z^{xy}_{10}-z^{xy}_{11},
+\end{aligned}\\
+a_{30}=2(z_{00}-z_{10})+z^x_{00}+z^x_{10},\quad
+a_{31}=2(z^y_{00}-z^y_{10})+z^{xy}_{00}+z^{xy}_{10},\\
+\begin{aligned}
+a_{32}=6(z_{01}+z_{10}-z_{00}-z_{11})+4(z^y_{10}-z^y_{00})&+3(z^x_{01}+z^x_{11}-z^x_{00}-z^x_{10})+{}\\
+{}&+2(z^y_{11}-z^y_{01}-z^{xy}_{00}-z^{xy}_{10})-z^{xy}_{01}-z^{xy}_{11},\\
+\end{aligned}\\
+\begin{aligned}
+a_{33}=4(z_{00}-z_{01}-z_{10}+z_{11})+2(z^x_{00}-z^x_{01}+z^x_{10}-z^x_{11}+z^y_{00}&-z^y_{10}+z^y_{01}-z^y_{11})+{}\\
+{}&+z^{xy}_{00}+z^{xy}_{01}+z^{xy}_{10}+z^{xy}_{11}.
+\end{aligned}
+\end{gathered}
+\label{bicub_coeff}
+\end{equation}
+%
+Из-за более сложных формул, бикубическая интерполяция работает существенно медленнее линейной.
+
+\subsection{Интерполяция во времени}
+Так как нам нужно получать значения скоростей не только в произвольном месте, но и в произвольное время, необходимо делать дополнительно интерполяцию во времени. Эта задача существенно проще пространственной интерполяции, так как в данном случае мы имеем дело с одномерной функцией. Простейший вид интерполяции~--- линейный
+%
+\begin{equation}
+z(x,y,t)=\frac{z(x,y,t_{k+1})-z(x,y,t_k)}{\Delta t} (t-t_k)+z(x,y,t_k).
+\label{lint}
+\end{equation}
+%
+где $\Delta t$~--- шаг сетки по времени.
+
+Мы также используем интерполяцию кубическими полиномами Лагранжа\footnote{Полином Лагранжа степени $n$~--- это такой полином, который в $n+1$ точке имеет заданные значения.}. Для $t\in[t_k:t_{k+1}]$ мы имеем
+%
+\begin{equation}
+z(x,y,t)=a_0(x,y)+a_1(x,y)t+a_2(x,y)t^2+a_3(x,y)t^3.
+\label{cubt}
+\end{equation}
+%
+Коэффициенты $a_i(x,y)$ рассчитываются из значений $z_i=z(x,y,t_i)$, которые, в свою очередь, находятся пространственной интерполяцией
+%
+\begin{equation}
+\begin{gathered}
+a_0=z_k,\qquad
+a_1=\left\{
+ \begin{aligned}
+ -&\frac{11}{6}z_k+3z_{k+1}-\frac{3}{2}z_{k+2}+\frac{1}{3}z_{k+3}, &k&=0,\\
+ -&\frac{1}{3}z_{k-1}-\frac{1}{2}z_k+z_{k+1}-\frac{1}{6}z_{k+2},   &k&\in[1:N_t-2],\\
+  &\frac{1}{6}z_{k-2}-z_{k-1}+\frac{1}{2}z_k+\frac{1}{3}z_{k+1},   &k&=N_t-1,
+ \end{aligned}\right.\\
+a_2=\left\{
+ \begin{aligned}
+  &z_k-\frac{5}{2}z_{k+1}+2z_{k+2}-\frac{1}{2}z_{k+3}, &k&=0,\\
+  &\frac{1}{2}z_{k-1}-z_k+\frac{1}{2}z_{k+1},          &k&\geqslant 1,
+ \end{aligned}\right.\quad
+a_3=\left\{
+ \begin{aligned}
+  -&\frac{1}{6}z_k+\frac{1}{2}z_{k+1}-\frac{1}{2}z_{k+2}+\frac{1}{6}z_{k+3}, &k&=0,\\
+  -&\frac{1}{6}z_{k-1}+\frac{1}{2}z_k-\frac{1}{2}z_{k+1}+\frac{1}{6}z_{k+2}, &k&\in[1:N_t-2],\\
+  -&\frac{1}{6}z_{k-2}+\frac{1}{2}z_{k-1}-\frac{1}{2}z_k+\frac{1}{6}z_{k+1}, &k&=N_t-1.
+ \end{aligned}\right.
+\end{gathered}
+\end{equation}
+%
+
+\section{Общий алгоритм}
+\begin{enumerate}
+\item Скорости и координаты преобразуются в прямоугольную сетку. Создаётся файл для интерполяции, содержащий информацию о сетке (число шагов $N_x$, $N_y$, $N_t$ и величина шага $\Delta x$, $\Delta y$, $\Delta t$ по каждой координате и времени) и скорости. Можно представить себе сетку как стопку плоскостей, расчерченых на прямоугольники. Каждая плоскость~--- это карта скоростей в фиксированый момент времени.
+\item Для расчёта мы используем либо линейную интерполяцию в пространстве и времени, либо бикубическую в пространстве и полиномиальную во времени. Каждая компонента скорости интерполируется независимо от другой.
+ \begin{enumerate}
+  \item Определяем ячейку $(i,j,k)$, в которой находится точка $(X,Y,T)$ и нормированые координаты внутри этой ячейки
+  \begin{equation}
+   \begin{aligned}
+   i&=[X/\Delta x], &j&=[Y/\Delta y], &k&=[T/\Delta t],\\ x&=\frac{X-i\Delta x}{\Delta x}, &y&=\frac{Y-j\Delta y}{\Delta y}, &t&=\frac{T-k\Delta t}{\Delta t},
+  \end{aligned}
+  \end{equation}
+  где $[r]$~--- целая часть числа $r$.
+  \item По формуле \eqref{linint} находим значения $u$ и $v$ на верхней и нижней гранях ячейки в точке $(x,y)$.
+  \item По формуле \eqref{lint} находим значения $u$ и $v$ в точке $(x,y,t)$.
+ \end{enumerate}
+Для второго метода алгоритм чуть сложнее.
+ \begin{enumerate}
+  \item Ячейка и нормированые координаты находятся так же, как и в первом случае.
+  \item Для четырёх временных плоскостей, ближайших к $T$, находятся по формуле \eqref{bicub_coeff} коэффициенты бикубического полинома \eqref{bicubic}\footnote{На самом деле они рассчитываются только один раз для каждого квадратика на плоскости, а потом хранятся в памяти. Пересчёт их при каждом определении скоростей~--- слишком дорогое удовольствие.}.
+  \item Рассчитываются по формуле \eqref{bicubic} значения скоростей в точке $(x,y)$ на четырёх временных плоскостях.
+  \item Скорости окончательно рассчитываются по формуле \eqref{cubt}.
+ \end{enumerate}
+\item Полученные скорости подставляются в уравнение \eqref{maineq} и интегрируются.
+\item Выходные данные либо используются как есть (для анализа), либо трансформируются обратно в географические координаты (для статей).
+\end{enumerate}
+\end{document}
diff --git a/Interpolation/net.eps b/Interpolation/net.eps
new file mode 100644
index 0000000..3002cfa
--- /dev/null
+++ b/Interpolation/net.eps
@@ -0,0 +1,139 @@
+%!PS-Adobe-2.0 EPSF-2.0
+%%Title: net.eps
+%%BoundingBox: 85 85 210 215
+%%EndComments
+gsave
+100 100 translate
+10 dict begin
+/l 100 def
+/r 3 def
+/mlw 0.5 def
+/r0 2 def
+/rs 1.5 def
+/x 0.45 def
+/y 0.7 def
+/dlw 0.5 def
+/dash [ 2 2 ] def
+/fs 15 def
+/sfs fs 0.5 mul def
+/dn fs -0.1 mul def
+/shb -4 def
+/shu 6 def
+
+/strparams {20 dict begin /inp exch def gsave newpath 0 0 moveto
+ /stringtype {inp false charpath} def
+ /arraytype
+ {
+  inp {dup dup 0 get findfont exch 1 get scalefont setfont
+   dup 2 get /x exch def dup 3 get /y exch def
+   dup 5 get {x y rmoveto 4 get false charpath}{pop} ifelse
+   } forall
+ } def
+ inp type exec pathbbox
+ /ury exch def /urx exch def /lly exch def /llx exch def
+ % Return values
+ urx llx sub ury lly sub llx neg lly neg
+grestore end} def
+/GetParam { 3 dict begin
+ /default exch def /key exch def /dict exch def
+ dict key known {dict key get}{default} ifelse
+ end
+} def
+/Show {20 dict begin /Params exch def /String exch def
+ /HAlign Params /HAlign -1 GetParam def
+ /VAlign Params /VAlign -1 GetParam def
+ /UseHShift Params /UseHShift false GetParam def
+ /UseVShift Params /UseVShift false GetParam def
+ String strparams /VShift exch def /HShift exch def /Height exch def /Width exch def
+ % Correct HAlign and VAlign
+ HAlign 0 lt {/HAlign -1 def}{} ifelse HAlign 0 gt {/HAlign 1 def}{} ifelse
+ VAlign 0 lt {/VAlign -1 def}{} ifelse VAlign 0 gt {/VAlign 1 def}{} ifelse
+ % Move initial point
+ UseHShift {HShift 0 rmoveto}{} ifelse UseVShift {0 VShift rmoveto}{} ifelse
+ HAlign 0 eq {Width -0.5 mul 0 rmoveto}{} ifelse HAlign 1 eq {Width neg 0 rmoveto}{} ifelse
+ VAlign 0 eq {0 Height -0.5 mul rmoveto}{} ifelse VAlign 1 eq {0 Height neg rmoveto}{} ifelse
+ /stringtype {show} def
+ /arraytype {{
+  dup dup 0 get findfont exch 1 get scalefont setfont currentpoint 3 -1 roll dup dup 2 get exch 3 get rmoveto dup 4 get show
+  dup length 5 eq {pop pop pop}{
+  dup length 6 eq {5 get {pop pop}{moveto} ifelse}{
+   dup 5 get exch 6 get currentpoint 4 -2 roll % x0 y0 x y ux uy
+   {1}{3}ifelse index exch % x0 y0 x y ye ux
+   {2}{4}ifelse index exch %x0 y0 x y xe ye
+   moveto pop pop pop pop
+  }ifelse}ifelse
+ }forall} def
+ String dup type exec
+end} def
+
+
+0 setgray
+
+newpath
+0 0 r 0 360 arc
+fill
+newpath
+0 l r 0 360 arc
+fill
+newpath
+l l r 0 360 arc
+fill
+newpath
+l 0 r 0 360 arc
+fill
+newpath
+x l mul y l mul r0 0 360 arc
+fill
+newpath
+x l mul 0 rs 0 360 arc
+fill
+newpath
+0 y l mul rs 0 360 arc
+fill
+
+mlw setlinewidth
+newpath
+0 0 moveto
+0 l lineto
+l l lineto
+l 0 lineto
+0 0 lineto
+stroke
+
+dlw setlinewidth
+dash 0 setdash
+newpath
+x l mul 0 moveto
+x l mul y l mul lineto
+0 y l mul lineto
+stroke
+
+0 0 moveto [
+[ /Times-Italic fs 0 shb (z) true]
+[ /Times-Italic sfs 0 dn (00) true]
+] << /HAlign 0 /VAlign 1>> Show
+l 0 moveto [
+[ /Times-Italic fs 0 shb (z) true]
+[ /Times-Italic sfs 0 dn (10) true]
+] << /HAlign 0 /VAlign 1>> Show
+0 l moveto [
+[ /Times-Italic fs 0 shu (z) true]
+[ /Times-Italic sfs 0 dn (01) true]
+] << /HAlign 0 /VAlign -1>> Show
+l l moveto [
+[ /Times-Italic fs 0 shu (z) true]
+[ /Times-Italic sfs 0 dn (11) true]
+] << /HAlign 0 /VAlign -1>> Show
+l x mul 0 moveto [
+[ /Times-Italic fs 0 shb (x) true]
+] << /HAlign 0 /VAlign 1>> Show
+0 l y mul moveto [
+[ /Times-Italic fs -3 0 (y) true]
+] << /HAlign 1 /VAlign 0>> Show
+l x mul l y mul moveto [
+[ /Times-Italic fs 3 3 (z) true]
+] << /HAlign -1 /VAlign -1>> Show
+
+end
+grestore
+%%EOF