Revert "Remove term from alpha that is equal to 0 anyways"

This reverts commit bc05a28b305f6689b86c64c8cdb794d9d8598744.

2 files changed, 30 insertions, 5 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 69b86f3..ab56a7e 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -132,7 +132,7 @@ int eea(
     ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
 
     energy=compute_energy(u, K1, K2);
-    alpha=compute_alpha(u, K1, K2, g, L);
+    alpha=compute_alpha(u, K1, K2, N1, N2, g, L, ifft.fft);
     enstrophy=compute_enstrophy(u, K1, K2, L);
     
     // running average
@@ -431,7 +431,7 @@ int ns_rhs(
   if (irreversible) {
     alpha=nu;
   } else {
-    alpha=compute_alpha(u,K1,K2,g,L);
+    alpha=compute_alpha(u,K1,K2,N1,N2,g,L,ifft.fft);
   }
 
   /*
@@ -586,8 +586,11 @@ double compute_alpha(
   _Complex double* u,
   int K1,
   int K2,
+  int N1,
+  int N2,
   _Complex double* g,
-  double L
+  double L,
+  _Complex double* T
 ){
   _Complex double num=0;
   double denom=0;
@@ -598,7 +601,7 @@ double compute_alpha(
 
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      num+=L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*getval_sym(g, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2));
+      num+=(L*L/4/M_PI/M_PI*(kx*kx+ky*ky)*getval_sym(g, kx,ky,K2)+sqrt(kx*kx+ky*ky)*T[klookup(kx,ky,N1,N2)])*conj(getval_sym(u, kx,ky,K2));
       denom+=__real__ (kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2));
     }
   }
@@ -607,6 +610,28 @@ double compute_alpha(
 }
 
 
+/*
+// compute alpha
+double compute_alpha(
+  _Complex double* u,
+  int K1,
+  int K2,
+  _Complex double* g
+){
+  _Complex double num=0;
+  _Complex double denom=0;
+  int kx,ky;
+
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(getval_sym(u, kx,ky,K2))*(1+(ky!=0?kx*kx/ky/ky:0));
+      num+=(kx*kx+ky*ky)*getval_sym(u, kx,ky,K2)*conj(g[getval(kx,ky,K1,K2)])*(1+(ky!=0?kx*kx/ky/ky:0));
+    }
+  }
+
+  return __real__ (num/denom);
+}*/
+
 // compute energy
 double compute_energy(
   _Complex double* u,
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 6a5f32b..da3894e 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -48,7 +48,7 @@ int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft
 int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
 
 // compute alpha
-double compute_alpha( _Complex double* u, int K1, int K2, _Complex double* g, double L);
+double compute_alpha( _Complex double* u, int K1, int K2, int N1, int N2, _Complex double* g, double L, _Complex double* T);
 // compute energy
 double compute_energy( _Complex double* u, int K1, int K2);
 // compute enstrophy