diff --git a/docs/nstrophy_doc.tex b/docs/nstrophy_doc.tex
index d7ad7b3..01ecc36 100644
--- a/docs/nstrophy_doc.tex
+++ b/docs/nstrophy_doc.tex
@@ -407,12 +407,11 @@ and, by\-~(\ref{T}),
     \frac{(\ell\cdot q^\perp)|\ell|}{|q|}
   \right)\hat u_q
   =
-  \sum_{\displaystyle\mathop{\scriptstyle q\in\mathbb Z^2}_{\ell+q=k}}
-  (q\cdot \ell^\perp)\left(
-    \frac{|q|}{|\ell|}
+  (k\cdot \ell^\perp)\left(
+    \frac{|k-\ell|}{|\ell|}
     -
-    \frac{|\ell|}{|q|}
-  \right)\hat u_q
+    \frac{|\ell|}{|k-\ell|}
+  \right)\hat u_{k-\ell}
   .
 \end{equation}
 \bigskip
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 9cfac70..9799fad 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -1343,6 +1343,79 @@ int ns_T_nofft(
   return 0;
 }
 
+
+// Jacobian of rhs
+int ns_D_rhs(
+  _Complex double* out,
+  _Complex double* u,
+  int K1,
+  int K2,
+  double nu,
+  double L,
+  bool irreversible
+){
+  int kx,ky,lx,ly;
+  int i;
+  double alpha;
+
+  if (irreversible) {
+    alpha=nu;
+  } else {
+    // TODO
+    alpha=0;
+  }
+
+  #define MATSIZE (K1*(2*(K2+1)+K2))
+  for(i=0; i<MATSIZE*MATSIZE; i++){
+    out[i]=0;
+  }
+  // diagonal term
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      out[klookup_sym(kx,ky,K2)*MATSIZE+klookup_sym(kx,ky,K2)]=-4*M_PI*M_PI/L/L*alpha*(kx*kx+ky*ky);
+    }
+  }
+  // convolution term
+  for(kx=0;kx<=K1;kx++){
+    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
+      for(lx=0;lx<=K1;lx++){
+	for(ly=(lx>0 ? -K2 : 1);ly<=K2;ly++){
+	  // no contribution from k=l
+	  if (kx!=lx && ky!=ly){
+	    // use row-major ordering
+	    out[klookup_sym(kx,ky,K2)*MATSIZE+klookup_sym(lx,ly,K2)]=4*M_PI*M_PI/L/L*ns_d_T(u,kx,ky,lx,ly,K2);
+	  }
+	}
+      }
+    }
+  }
+
+  return(0);
+}
+
+// derivative of T with respect to u_k
+_Complex double ns_d_T(
+  _Complex double* u,
+  int kx,
+  int ky,
+  int lx,
+  int ly,
+  int K2
+){
+  int qx=kx-lx;
+  int qy=ky-ly;
+  // trivial case
+  if (qx*qx+qy*qy==0){
+    return 0.;
+  }
+  if (qx>0 || (qx==0 && qy>=1)){
+    return (ky*lx-kx*ly)*(qx*qx+qy*qy-lx*lx-ly*ly)/sqrt(lx*lx+ly*ly)/sqrt(qx*qx+qy*qy)*u[klookup_sym(qx,qy,K2)];
+  }else{
+    return (ky*lx-kx*ly)*(qx*qx+qy*qy-lx*lx-ly*ly)/sqrt(lx*lx+ly*ly)/sqrt(qx*qx+qy*qy)*u[klookup_sym(-qx,-qy,K2)];
+  }
+}
+
+
 // compute alpha
 double compute_alpha(
   _Complex double* u,
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 34a61d0..37fd076 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -71,6 +71,11 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
 // convolution term in right side of equation (computed without fft)
 int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
 
+// Jacobian of rhs
+int ns_D_rhs( _Complex double* out, _Complex double* u, int K1, int K2, double nu, double L, bool irreversible);
+// derivative of T with respect to u_k
+_Complex double ns_d_T( _Complex double* u, int kx, int ky, int lx, int ly, int K2);
+
 // compute alpha
 double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g, double L);
 // compute energy