From cf90b764f8253be64e03e3f604f8d08163d0d934 Mon Sep 17 00:00:00 2001
From: Martin	Diehl <m.diehl@mpie.de>
Date: Sat, 5 Sep 2015 11:51:36 +0000
Subject: [PATCH] calculation was for nyquist freq not fully correct. See Notes
 on FFT-based differentiation Steven G. Johnson, MIT Applied Mathematics
 Created April, 2011, updated May 4, 2011:

---
 processing/post/addCurl.py       | 18 ++++++++++++++----
 processing/post/addDivergence.py | 18 ++++++++++++++----
 2 files changed, 28 insertions(+), 8 deletions(-)

diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py
index afc060d1d..1b189bccf 100755
--- a/processing/post/addCurl.py
+++ b/processing/post/addCurl.py
@@ -29,13 +29,23 @@ def curlFFT(geomdim,field):
  TWOPIIMG = (0.0+2.0j*math.pi)
  for i in xrange(grid[0]):
    k_s[0] = i
-   if(i > grid[0]/2 ): k_s[0] = k_s[0] - grid[0]
+   if(grid[0]%2==0 and i == grid[0]//2):                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+     k_s[0]=0
+   elif (i > grid[0]//2): 
+     k_s[0] = k_s[0] - grid[0]
+
    for j in xrange(grid[1]):
      k_s[1] = j
-     if(j > grid[1]/2 ): k_s[1] = k_s[1] - grid[1]
-     for k in xrange(grid[2]/2+1):
+     if(grid[1]%2==0 and j == grid[1]//2):                                                          # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+       k_s[1]=0
+     elif (j > grid[1]//2): 
+       k_s[1] = k_s[1] - grid[1]
+
+     for k in xrange(grid[2]//2+1):
        k_s[2] = k
-       if(k > grid[2]/2 ): k_s[2] = k_s[2] - grid[2]
+       if(grid[2]%2==0 and k == grid[2]//2):                                                        # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+         k_s[2]=0
+
        xi = np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16')
        if dataType == 'tensor': 
          for l in xrange(3):
diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py
index 32e0bb120..27d616aad 100755
--- a/processing/post/addDivergence.py
+++ b/processing/post/addDivergence.py
@@ -24,13 +24,23 @@ def divFFT(geomdim,field):
  TWOPIIMG = (0.0+2.0j*math.pi)
  for i in xrange(grid[0]):
    k_s[0] = i
-   if(i > grid[0]/2 ): k_s[0] = k_s[0] - grid[0]
+   if(grid[0]%2==0 and i == grid[0]//2):                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+     k_s[0]=0
+   elif (i > grid[0]//2): 
+     k_s[0] = k_s[0] - grid[0]
+
    for j in xrange(grid[1]):
      k_s[1] = j
-     if(j > grid[1]/2 ): k_s[1] = k_s[1] - grid[1]
-     for k in xrange(grid[2]/2+1):
+     if(grid[1]%2==0 and j == grid[1]//2):                                                          # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+       k_s[1]=0
+     elif (j > grid[1]//2): 
+       k_s[1] = k_s[1] - grid[1]
+
+     for k in xrange(grid[2]//2+1):
        k_s[2] = k
-       if(k > grid[2]/2 ): k_s[2] = k_s[2] - grid[2]
+       if(grid[2]%2==0 and k == grid[2]//2):                                                        # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+         k_s[2]=0
+
        xi=np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16')
        if n == 9:                                                                                   # tensor, 3x3 -> 3
          for l in xrange(3):