diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex
index 71825dc..478a4c9 100644
--- a/docs/nstrophy_doc.tex
+++ b/docs/nstrophy_doc.tex
@@ -301,7 +301,7 @@ in which $\mathcal F^*$ is defined like $\mathcal F$ but with the opposite phase
   \sum_{p,q\in\mathcal K}
   \frac1{N_1N_2}
   \sum_{n\in\mathcal N}e^{-\frac{2i\pi}{N_1}n_1(p_1+q_1-k_1)-\frac{2i\pi}{N_2}n_2(p_2+q_2-k_2)}
-  (q\cdot p^\perp)\frac{|q|}{|p|}\hat u_q\hat u_p
+  (q\cdot p^\perp)\frac{|q|}{|p|}\hat u_p\hat u_q
 \end{equation}
 provided
 \begin{equation}
@@ -440,52 +440,98 @@ The enstrophy is defined as
 
 \subsection{Lyapunov exponents}
 \indent
-We will consider several methods of computing the Lyapunov exponents.
+The Lyapunov are defined from the {\it tangent flow} of the dynamics.
+Consider an equation of the form
+\begin{equation}
+  \dot u=f(t;u)
+  .
+\end{equation}
+The tangent flow is given by
+\begin{equation}
+  \dot\delta=Df(t;u(t))\delta
+\end{equation}
+where $Df$ is the Jacobian of $f$ with respect to $u$.
+The flow of this equation is denoted by
+\begin{equation}
+  \varphi_{t_0,t_1}(\delta_0)
+\end{equation}
+and defined by
+\begin{equation}
+  \frac d{dt}\varphi_{t_0,t}(\delta_0)=Df(t;\varphi_{t_0,t}(\delta_0))\varphi_{t_0,t}(\delta_0)
+  ,\quad
+  \varphi_{t_0,t_0}(\delta_0)=\delta_0
+  .
+\end{equation}
+The flow $\varphi_{t_0,t}(\delta_0)$ is linear in $\delta_0$, and so can be represented as a matrix.
+The Lyapnuov exponents are defined from the $QR$ decomposition of $\varphi_{0,t}$: if
+\begin{equation}
+  \varphi_{0,t}=Q_tR_t
+\end{equation}
+where $Q_t$ is orthogonal, and $R_t$ is triangular with positive diagonal entries $(r_t^{(j)})_j$, then the Lyapunov exponents are
+\begin{equation}
+  \lim_{t\to\infty}\frac 1t\log|r_t^{(j)}|
+  .
+\end{equation}
 \bigskip
 
-\subsubsection{Method 1: compute the Jacobian at every step}
 \indent
-To compute the Lyapunov exponents, we must first compute the Jacobian of $\hat u^{(n)}\mapsto\hat u^{(n+1)}$.
-This map is always of Runge-Kutta type, that is,
+One problem to compute the Lyapunov exponents numerically is that the spectrum depends exponentially on time, and so has a tendency to blow up or shrink very rapidly, which leads to large truncation errors.
+To avoid these, we compute the tangent flow {\it in parts}: we split time into intervals:
 \begin{equation}
-  \hat u(t_{n+1})=\mathfrak F_{t_n}(\hat u(t_n))
-  .
-\end{equation}
-Let $D\mathfrak F_{t_{n}}$ be the Jacobian of this map, in which we split the real and imaginary parts: if
-\begin{equation}
-  \hat u_k(t_n)=:\rho_k+i\iota_k
+  [0,t)=\bigcup_{i=0}^{N-1}[t_i,t_{i+1})
   ,\quad
-  \mathfrak F_{t_n}(\hat u(t_n))_k=:\phi_k+i\psi_k
-\end{equation}
-then
-\begin{equation}
-  (D\mathfrak F_{t_n})_{k,p}:=\left(\begin{array}{cc}
-    \partial_{\rho_p}\phi_k&\partial_{\iota_p}\phi_k\\
-    \partial_{\rho_p}\psi_k&\partial_{\iota_p}\psi_k
-  \end{array}\right)
-\end{equation}
-We compute this Jacobian numerically using a finite difference, by computing
-\begin{equation}
-  (D\mathfrak F_{t_n})_{k,p}:=\frac1\epsilon
-  \left(\begin{array}{cc}
-    \phi_k(\hat u+\epsilon\delta_p)-\phi_k(\hat u)&\phi_k(\hat u+i\epsilon\delta_p)-\phi_k(\hat u)\\
-    \psi_k(\hat u+\epsilon\delta_p)-\psi_k(\hat u)&\psi_k(\hat u+i\epsilon\delta_p)-\psi_k(\hat u)
-  \end{array}\right)
+  \varphi_{0,t}=\prod_{i=N-1}^0\varphi_{t_{i},t_{i+1}}
+  .
+\end{equation}
+At the thresholds between these intervals, we perform a QR decomposition:
+\begin{equation}
+  Q_0R_0:=\varphi_{0,t_0}
+  ,\quad
+  Q_iR_i:=\varphi_{t_{i-1},t_i}Q_{i-1}
+\end{equation}
+where $Q_i$ is orthogonal and $R_i$ is upper triangular with $\geqslant 0$ diagonal entries (in addition, we assume these are not $0$).
+Doing so, we find
+\begin{equation}
+  \varphi_{0,t}=Q_{N-1}\prod_{i=N-1}^0R_i
+\end{equation}
+and since the product of triangular matrices is triangular and the diagonal elements multiply, we find that the Lyapunov exponents are given by
+\begin{equation}
+  \lim_{t_{N-1}\to\infty}\frac1{t_{N-1}}\sum_{i=N-1}^0\log|r_i^{(j)}|
+  .
+\end{equation}
+To do the computation numerically, we drop the limit, and compute the logarithm of the product of the diagonal entries of $R_i$.
+\bigskip
+
+\indent
+In practice, we approximate $\varphi_{t_{i-1},t_i}$ by running a Runge-Kutta algorithm for the tangent flow equation.
+To obtain the full matrix, we consider every element of the canonical basis as an initial condition $\delta_0$.
+We then iterate the Runge-Kutta algorithm until the time $t_0$ (chosen in one of two ways, see below), at which point we perform a QR decomposition, save the diagonal entries of $R$, replace the family of initial conditions with the columns of $Q$, and continue the flow from there.
+The choice of the times $t_i$ can be done either by fixed-length intervals, specified with the option {\tt lyapunov\_reset}, or the QR decomposition can be triggered whenever $\|\delta\|_1$ exceeds a threshold, specified in {\tt lyapunov\_maxu} (after all, the intervals are used to prevent $\delta$ from becoming too large). 
+\bigskip
+
+\indent
+To compute the Lyapunov exponents, we thus need the Jacobian of $f$.
+For the irreversible equation,
+\begin{equation}
+  f(\hat u)=
+  -\frac{4\pi^2}{L^2}\nu k^2\hat u_k+\hat g_k
+  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
+\end{equation}
+and so
+\begin{equation}
+  ((D f(\hat u))\delta)_k
+  =
+  -\frac{4\pi^2}{L^2}\nu k^2\delta_k
+  +\frac{4\pi^2}{L^2|k|}DT(\hat u,k)\delta
+\end{equation}
+where
+\begin{equation}
+  DT(\hat u,k)\delta
+  =
+  \sum_{\displaystyle\mathop{\scriptstyle p,q\in\mathbb Z^2}_{p+q=k}}
+  \left(\frac{(q\cdot p^\perp)|q|}{|p|}+\frac{(p\cdot q^\perp)|p|}{|q|}\right)\hat u_p\hat \delta_q
   .
 \end{equation}
-The parameter $\epsilon$ can be set using the parameter {\tt D\_epsilon}.
-%, so
-%\begin{equation}
-%  D\hat u^{(n+1)}=\mathds 1+\delta\sum_{i=1}^q w_iD\mathfrak F(\hat u^{(n)})
-%  .
-%\end{equation}
-%We then compute
-%\begin{equation}
-%  (D\mathfrak F(\hat u))_{k,\ell}
-%  =
-%  -\frac{4\pi^2}{L^2}\nu k^2\delta_{k,\ell}
-%  +\frac{4\pi^2}{L^2|k|}\partial_{\hat u_\ell}T(\hat u,k)
-%\end{equation}
 %and, by\-~(\ref{T}),
 %\begin{equation}
 %  \partial_{\hat u_\ell}T(\hat u,k)
@@ -504,30 +550,115 @@ The parameter $\epsilon$ can be set using the parameter {\tt D\_epsilon}.
 %  \right)\hat u_{k-\ell}
 %  .
 %\end{equation}
-\bigskip
+For the reversible equation,
+\begin{equation}
+  f(\hat u)=
+  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\hat u_k
+  +\hat g_k
+  +\frac{4\pi^2}{L^2|k|}T(\hat u,k)
+\end{equation}
+so
+\begin{equation}
+  ((D f(\hat u))\delta)_k
+  =
+  -\frac{4\pi^2}{L^2}\alpha(\hat u) k^2\delta_k
+  -\frac{4\pi^2}{L^2}k^2\hat u_k D\alpha(\hat u)\delta
+  +\frac{4\pi^2}{L^2|k|}DT(\hat u,k)\delta
+\end{equation}
+where
+\begin{equation}
+  D\alpha(\hat u)\delta
+  =
+  \frac{\frac{L^2}{4\pi^2}\sum_k k^2\hat \delta_k^*\hat g_k+\sum_k|k|\hat (\delta_k^*T(\hat u,k)+\hat u_k^*DT(\hat u,k)\delta)}{\sum_kk^4|\hat u_k|^2}
+  -\alpha(\hat u)\frac{2\sum_kk^4\delta_k^*\hat u_k}{\sum_kk^4|\hat u_k|^2}
+  .
+\end{equation}
 
-\indent
-The Lyapunov exponents are then defined as
-\begin{equation}
-  \frac1{t_n}\log\mathrm{spec}(\pi_{t_n})
-  ,\quad
-  \pi_{t_n}:=\prod_{i=1}^nD\hat u(t_i)
-  .
-\end{equation}
-However, the product of $D\hat u$ may become quite large or quite small (if the exponents are not all 1).
-To avoid this, we periodically rescale the product.
-We set $\mathfrak L_r>0$ (set by adjusting the {\tt lyanpunov\_reset} parameter), and, when $t_n$ crosses a multiple of $\mathfrak L_r$, we rescale the eigenvalues of $\pi_i$ to 1.
-To do so, we perform a $QR$ decomposition:
-\begin{equation}
-  \pi_{\alpha\mathfrak L_r}
-  =R^{(\alpha)}Q^{(\alpha)}
-\end{equation}
-where $Q^{(\alpha)}$ is orthogonal and $R^{(\alpha)}$ is a diagonal matrix, and we divide by $R^{(\alpha)}$ (thus only keeping $Q^{(\alpha)}$).
-The Lyapunov exponents at time $\alpha\mathfrak L_r$ are then
-\begin{equation}
-  \frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)})
-  .
-\end{equation}
+
+%\indent
+%To compute the Lyapunov exponents, we must first compute the Jacobian of $\hat u^{(n)}\mapsto\hat u^{(n+1)}$.
+%This map is always of Runge-Kutta type, that is,
+%\begin{equation}
+%  \hat u(t_{n+1})=\mathfrak F_{t_n}(\hat u(t_n))
+%  .
+%\end{equation}
+%Let $D\mathfrak F_{t_{n}}$ be the Jacobian of this map, in which we split the real and imaginary parts: if
+%\begin{equation}
+%  \hat u_k(t_n)=:\rho_k+i\iota_k
+%  ,\quad
+%  \mathfrak F_{t_n}(\hat u(t_n))_k=:\phi_k+i\psi_k
+%\end{equation}
+%then
+%\begin{equation}
+%  (D\mathfrak F_{t_n})_{k,p}:=\left(\begin{array}{cc}
+%    \partial_{\rho_p}\phi_k&\partial_{\iota_p}\phi_k\\
+%    \partial_{\rho_p}\psi_k&\partial_{\iota_p}\psi_k
+%  \end{array}\right)
+%\end{equation}
+%We compute this Jacobian numerically using a finite difference, by computing
+%\begin{equation}
+%  (D\mathfrak F_{t_n})_{k,p}:=\frac1\epsilon
+%  \left(\begin{array}{cc}
+%    \phi_k(\hat u+\epsilon\delta_p)-\phi_k(\hat u)&\phi_k(\hat u+i\epsilon\delta_p)-\phi_k(\hat u)\\
+%    \psi_k(\hat u+\epsilon\delta_p)-\psi_k(\hat u)&\psi_k(\hat u+i\epsilon\delta_p)-\psi_k(\hat u)
+%  \end{array}\right)
+%  .
+%\end{equation}
+%The parameter $\epsilon$ can be set using the parameter {\tt D\_epsilon}.
+%%, so
+%%\begin{equation}
+%%  D\hat u^{(n+1)}=\mathds 1+\delta\sum_{i=1}^q w_iD\mathfrak F(\hat u^{(n)})
+%%  .
+%%\end{equation}
+%%We then compute
+%%\begin{equation}
+%%  (D\mathfrak F(\hat u))_{k,\ell}
+%%  =
+%%  -\frac{4\pi^2}{L^2}\nu k^2\delta_{k,\ell}
+%%  +\frac{4\pi^2}{L^2|k|}\partial_{\hat u_\ell}T(\hat u,k)
+%%\end{equation}
+%%and, by\-~(\ref{T}),
+%%\begin{equation}
+%%  \partial_{\hat u_\ell}T(\hat u,k)
+%%  =
+%%  \sum_{\displaystyle\mathop{\scriptstyle q\in\mathbb Z^2}_{\ell+q=k}}
+%%  \left(
+%%    \frac{(q\cdot \ell^\perp)|q|}{|\ell|}
+%%    +
+%%    \frac{(\ell\cdot q^\perp)|\ell|}{|q|}
+%%  \right)\hat u_q
+%%  =
+%%  (k\cdot \ell^\perp)\left(
+%%    \frac{|k-\ell|}{|\ell|}
+%%    -
+%%    \frac{|\ell|}{|k-\ell|}
+%%  \right)\hat u_{k-\ell}
+%%  .
+%%\end{equation}
+%\bigskip
+%
+%\indent
+%The Lyapunov exponents are then defined as
+%\begin{equation}
+%  \frac1{t_n}\log\mathrm{spec}(\pi_{t_n})
+%  ,\quad
+%  \pi_{t_n}:=\prod_{i=1}^nD\hat u(t_i)
+%  .
+%\end{equation}
+%However, the product of $D\hat u$ may become quite large or quite small (if the exponents are not all 1).
+%To avoid this, we periodically rescale the product.
+%We set $\mathfrak L_r>0$ (set by adjusting the {\tt lyanpunov\_reset} parameter), and, when $t_n$ crosses a multiple of $\mathfrak L_r$, we rescale the eigenvalues of $\pi_i$ to 1.
+%To do so, we perform a $QR$ decomposition:
+%\begin{equation}
+%  \pi_{\alpha\mathfrak L_r}
+%  =R^{(\alpha)}Q^{(\alpha)}
+%\end{equation}
+%where $Q^{(\alpha)}$ is orthogonal and $R^{(\alpha)}$ is a diagonal matrix, and we divide by $R^{(\alpha)}$ (thus only keeping $Q^{(\alpha)}$).
+%The Lyapunov exponents at time $\alpha\mathfrak L_r$ are then
+%\begin{equation}
+%  \frac1{\alpha\mathfrak L_r}\sum_{\beta=1}^\alpha\log\mathrm{spec}(Q^{(\beta)})
+%  .
+%\end{equation}