Fix 2 MAJOR bugs with the gaussfield module: Nyquist and numpy

1. All the Nyquist frequencies MUST be zeroed in the wave-number
   representation for the weighted Helmholtz decomposition to work.

2. The fftn and ifftn calls were including all 4 dimensions of the
   array, it should have been the last 3.

The issue with Nyquist frequencies not being zeroed is also present in
cow and needs to be corrected. It has no major influence on the
decomposed power spectra, but it does invalidate the behavior of the
rmdivB option.
Add reductions.py and toolbase.py