def SecondOrderCorrFuncfCoherentXrayScattering(self,state,gridShape):
        """
        Returns the intensity of scattered light as a function
		of the reciprocal lattice vector q

		.. math::
			I(q) = \\frac{1}{V} \\tilde{f}(q) \\tilde{f}(-q)  
        
		where 
		
		.. math::
			f(x) = e^{2 \pi i g*u(x)}

        Assume that there is a cubic symmetry, so that g = (h/a, k/a, l/a)  (set a=1)
        """
        V = float(array(gridShape).prod())
        u = state#.CalculateDisplacementField()
        ug = self.g[0]*u[x]+self.g[1]*u[y]+self.g[2]*u[z] 
        Kug = fft.fftn(ug) 
        u2g = ug*ug
        Ku2g = fft.fftn(u2g)
        I_0 = 1.
        I_1 = (1.j)*2.*pi*(Kug-Kug.conj())
        I_2 = 2.*pi**2*(2*Kug*Kug.conj()-Ku2g-Ku2g.conj())
        I_S = (I_0+I_2)/V
        I_A = I_1/V
        return I_S.real, I_A.real

def garfield(B, P, T = Identity(), seed = None):
	if seed != None:
		random.seed(seed)
	wn = random.normal(0, 1, B.shape)
	f  = fft.ifftn(fft.fftn(wn) * np.sqrt(P(B.K))).real
	#f /= f.std()
	return fft.ifftn(fft.fftn(f) * T(B.K)).real

    def potential_energy(self, potential, summed=False):
        r"""Calculate the potential energy :math:`E_{\text{pot}} := \langle\Psi|V|\Psi\rangle`
        of the different components :math:`\psi_i`.

        :param potential: The potential energy operator :math:`V(x)`.
        :param summed: Whether to sum up the potential energies :math:`E_i` of the individual
                       components :math:`\psi_i`. Default is ``False``.
        :return: A list with the potential energies of the individual components
                 or the overall potential energy of the wavefunction. (Depending on the optional arguments.)
        """
        # Compute the prefactor
        T = self._grid.get_extensions()
        N = self._grid.get_number_nodes()
        prefactor = product(array(T) / (1.0 * array(N) ** 2))

        # Reshape from (1, prod_d^D N_d) to (N_1, ..., N_D) shape
        potential = [pot.reshape(N) for pot in potential]

        # Apply the matrix potential to the ket
        tmp = [zeros(component.shape, dtype=complexfloating) for component in self._values]
        for row in xrange(0, self._number_components):
            for col in xrange(0, self._number_components):
                tmp[row] = tmp[row] + potential[row * self._number_components + col] * self._values[col]

        # Fourier transform the components
        ftcbra = [fftn(component) for component in self._values]
        ftcket = [fftn(component) for component in tmp]

        # Compute the braket in Fourier space
        epot = [prefactor * sum(conjugate(cbra) * cket) for cbra, cket in zip(ftcbra, ftcket)]

        if summed is True:
            epot = sum(epot)

        return epot

    def solve_neutral(self, phi_g, rho_g, eps=None):
        """Solve Poissons equation for a neutral and periodic charge density.

        Parameters
        ----------
        phi_g: ndarray
            Potential (output array).
        rho_g: ndarray
            Charge distribution (in units of -e).

        """

        assert phi_g.dtype == self.dtype
        assert rho_g.dtype == self.dtype
        
        if self.gd.comm.size == 1:
            # Note, implicit downcast from complex to float when the dtype of
            # phi_g is float
            phi_g[:] = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q)
        else:
            rho_g = self.gd.collect(rho_g)
            if self.gd.comm.rank == 0:
                globalphi_g = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q)
            else:
                globalphi_g = None
            # What happens here if globalphi is complex and phi is real ??????
            self.gd.distribute(globalphi_g, phi_g)
            
        return 1

def gaussian_convolution(data, ijk_linewidths):

  from numpy import float32, zeros, add, divide, outer, reshape
  if data.dtype.type != float32:
    data = data.astype(float32)

  from math import exp
  gaussians = []
  for a in range(3):
    size = data.shape[a]
    gaussian = zeros((size,), float32)
    hw = ijk_linewidths[2-a] / 2.0
    for i in range(size):
      u = min(i,size-i) / hw
      p = min(u*u/2, 100)               # avoid OverflowError with exp()
      gaussian[i] = exp(-p)
    area = add.reduce(gaussian)
    divide(gaussian, area, gaussian)
    gaussians.append(gaussian)

  g01 = outer(gaussians[0], gaussians[1])
  g012 = outer(g01, gaussians[2])
  g012 = reshape(g012, data.shape)
  
  cdata = zeros(data.shape, float32)

  from numpy.fft import fftn, ifftn
  # TODO: Fourier transform Gaussian analytically to reduce computation time
  #       about 30% (one of three fft calculations).
  ftg = fftn(g012)
  ftd = fftn(data)
  gd = ifftn(ftg * ftd)
  gd = gd.astype(float32)
  return gd

def test_plan_call():
    for shape in tested_shapes:
        plan = Plan(
            input_array=ranf_unit_complex(shape),
            output_array=numpy.empty(shape, dtype=numpy.complex),
            direction=Direction.forward,
        )
        testing.assert_allclose(
            plan(),
            fft.fftn(plan.input_array)
        )
        testing.assert_allclose(
            plan(normalize=True),
            fft.fftn(plan.input_array) / plan.input_array.size
        )
        plan = Plan(
            input_array=ranf_unit_complex(shape),
            output_array=numpy.empty(shape, dtype=numpy.complex),
            direction=Direction.backward
        )
        testing.assert_allclose(
            plan(),
            fft.ifftn(plan.input_array) * plan.input_array.size
        )
        testing.assert_allclose(
            plan(normalize=True),
            fft.ifftn(plan.input_array)
        )

    def fft_r2g(self, fr, shift_fg=False):
        """
        FFT of array ``fr`` given in real space.
        """
        ndim, shape = fr.ndim, fr.shape

        if ndim == 1:
            fr = np.reshape(fr, self.shape)
            return self.fft_r2g(fr, shift_fg=shift_fg).flatten()

        elif ndim == 3:
            assert self.size == np.prod(shape[-3:])
            fg = fftn(fr)
            if shift_fg: fg = fftshift(fg)

        elif ndim > 3:
            assert self.size == np.prod(shape[-3:])
            axes = np.arange(ndim)[-3:]
            fg = fftn(fr, axes=axes)
            if shift_fg: fg = fftshift(fg, axes=axes)

        else:
            raise NotImplementedError("ndim < 3 are not supported")

        return fg / self.size

def fftconvolve(in1, in2, mode='same'):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1 + s2 - 1

    # Always use 2**n-sized FFT
    fsize = (2 ** np.ceil(np.log2(size))).astype('int')
    IN1 = fftn(in1, fsize)
    IN1 *= fftn(in2, fsize)
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = ifftn(IN1)[fslice].copy()
    del IN1 
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if np.product(s1, axis=0) > np.product(s2, axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret, osize)
    elif mode == "valid":
        return _centered(ret, abs(s2 - s1) + 1)

    return conv[:s[0], :s[1], :s[2]]

def nccfft(s, p, room, fact=1):
    """Used for all Patterns that do not fall other categories.

    Cross correlates normalized Source and Pattern images while
    taking advantage of FFTs for the convolution step. 
    
    ----------
    s, p : Image
       Pattern and Source images for comparison.
       
    fact : int, optional
       Factor by which both Source and Pattern are
       scaled down.

    ----------
    out1 : ndarray[float]
       Confidence matrix for matches.

    out2 : float
       Threshold for deciding if a match has been found.

    out3 : float
       Mean of the confidence matrix.
    """
    
    # subtract mean from Pattern
    pmm = p - p.mean()
    pstd = p.std()
    n = p.size
    
    # make matrix of ones the same size as pattern
    u = np.ones(p.shape)

    # pad matrices (necessary for convolution)
    s = pad_by(s, room)

    upad = pad_to_size_of(u, s)
    pmmpad = pad_to_size_of(pmm, s)
    
    # compute neccessary ffts
    fftppad = fftn(pmmpad)
    ffts = fftn(s)
    fftss = fftn(s**2)
    fftu = fftn(upad)

    # compute conjugates
    cfppad = np.conj(fftppad)
    cfu = np.conj(fftu)

    # do multiplications and ifft's
    top = ifftn(cfppad * ffts)
    bot1 = n * ifftn(cfu * fftss)
    bot2 = ifftn(cfu * ffts) ** 2

    # finish it off!
    bottom = pstd * np.sqrt(bot1 - bot2)
    full = top / bottom

    return np.where(full.real.max() == full.real)

def FieldFFT(a, a_hat):
	"""Calculate the component-wise 2D or 3D FFT of a vector field a, and stores it in a_hat."""
	
	if DimOfVectorFieldDomain(a) == 2:
		a_hat[:,:,0], a_hat[:,:,1] = fft.fftn(a[:,:,0]), fft.fftn(a[:,:,1])
	else:
		a_hat[...,0], a_hat[...,1], a_hat[...,2] = fft.fftn(a[...,0]), fft.fftn(a[...,1]), fft.fftn(a[...,2])

	return a_hat

def fftd(I, dims=None):

    # Compute fft
    if dims is None:
        X = fftn(I)
    elif dims == 2:
        X = fft2(I, axes=(0, 1))
    else:
        X = fftn(I, axes=tuple(range(dims)))

    return X

def wiener_deconvolve(x, kernel, l):
    s1 = array(x.shape)
    s2 = array(kernel.shape)
    size = s1 + s2 - 1

    X = fftn(x, size)
    H = fftn(kernel, size)
    Hc = np.conj(H)

    ret = ifftshift(ifftn( X*Hc / (H*Hc*X + l**2))).real
    return _centered(ret, s1)

 def solve_neutral(self, phi_g, rho_g, eps=None):
     if self.gd.comm.size == 1:
         phi_g[:] = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q).real
     else:
         rho_g = self.gd.collect(rho_g)
         if self.gd.comm.rank == 0:
             globalphi_g = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q).real
         else:
             globalphi_g = None
         self.gd.distribute(globalphi_g, phi_g)
     return 1

    def peerWindow(im, bounds):
        im2 = deepcopy(im)
        
        im2.arr = im2.arr[bounds[2]:bounds[3],
                          bounds[0]:bounds[1]]
        im2.warr = im2.warr[bounds[2]:bounds[3],
                            bounds[0]:bounds[1], :]
        im2.fft = fftn(im2.arr)
        im2.fft2 = fftn(im2.arr ** 2)
        im2.stdev = im2.arr.std()
        im2.mean = im2.arr.mean()

        return im2

    def propagate(self):
        r"""Given the wavefunction values :math:`\Psi(\Gamma)` at time :math:`t`, calculate
        new values :math:`\Psi^\prime(\Gamma)` at time :math:`t + \tau`. We perform exactly
        one single timestep of size :math:`\tau` within this function.
        """
        # How many components does Psi have
        N = self._psi.get_number_components()

        # Unpack the values from the current WaveFunction
        values = self._psi.get_values()

        # First step with the potential
        tmp = [zeros(value.shape, dtype=complexfloating) for value in values]
        for row in range(0, N):
            for col in range(0, N):
                tmp[row] = tmp[row] + self._VE[row * N + col] * values[col]

        # Go to Fourier space
        tmp = [fftn(component) for component in tmp]

        # First step with the kinetic operator
        tmp = [self._TE * component for component in tmp]

        # Go back to real space
        tmp = [ifftn(component) for component in tmp]

        # Central step with V-tilde
        tmp2 = [zeros(value.shape, dtype=complexfloating) for value in values]
        for row in range(0, N):
            for col in range(0, N):
                tmp2[row] = tmp2[row] + self._VEtilde[row * N + col] * tmp[col]

        # Go to Fourier space
        tmp = [fftn(component) for component in tmp2]

        # Second step with the kinetic operator
        tmp = [self._TE * component for component in tmp]

        # Go back to real space
        tmp = [ifftn(component) for component in tmp]

        # Second step with the potential
        values = [zeros(component.shape, dtype=complexfloating) for component in tmp]
        for row in range(0, N):
            for col in range(0, N):
                values[row] = values[row] + self._VE[row * N + col] * tmp[col]

        # Pack values back to WaveFunction object
        # TODO: Consider squeeze(.) of data before repacking
        self._psi.set_values(values)

def convolve(target, kernel, center=None):
    if center is None:
        center = np.array(kernel.shape)/2

    shape = [max(ts, ks) for ts, ks in zip(target.shape, kernel.shape)]
    f_kernel = fft.fftn(kernel, shape)
    f_input = fft.fftn(target, shape)
    f_result = f_kernel*f_input
    result = fft.ifftn(f_result)

    for i in range(len(center)):
        roll = -center[i]
        result = np.roll(result, roll, axis=i)

    return result

    def emb_fftn(self, input_x, output_dim, act_axes):
        '''
        embedded fftn: abstraction of fft for future gpu computing 
        '''
        output_x=numpy.zeros(output_dim, dtype=dtype)
        #print('output_dim',input_dim,output_dim,range(0,numpy.size(input_dim)))
#         output_x[[slice(0, input_x.shape[_ss]) for _ss in range(0,len(input_x.shape))]] = input_x
        output_x[crop_slice_ind(input_x.shape)] = input_x
#         print('GPU flag',self.gpu_flag)
#         print('pyfftw flag',self.pyfftw_flag)
#         if self.gpu_flag == 1:
#             self.thr.to_device(output_x.astype(dtype), dest=self.data_dev)
#             output_x=self.gpufftn(self.data_dev).get()
             
#         except:
#         elif self.gpu_flag ==0:
#         elif self.pyfftw_flag == 1:
# #             try:
#     #                 print('using pyfftw interface')
#     #                 print('threads=',self.threads)
#             output_x=pyfftw.interfaces.scipy_fftpack.fftn(output_x, output_dim, act_axes, 
#                                                               threads=self.threads,overwrite_x=True)
# #             except: 
#         else:
    #                 print('using OLD interface')                
#                 output_x=scipy.fftpack.fftn(output_x, output_dim, act_axes,overwrite_x=True)
        output_x=fftpack.fftn(output_x, output_dim, act_axes)


        return output_x

def perdecomp (image):
  # Compute boundary image
  h,w,d = image.shape
  v = zeros (image.shape)
  v[:,0,:] = v[:,0,:] + image[:,0,:] - image[:,w-1,:]
  v[:,w-1,:] = v[:,w-1,:] + image[:,w-1,:] - image[:,0,:]
  v[0,:,:] = v[0,:,:] + image[0,:,:] - image[h-1,:,:]
  v[h-1,:,:] = v[h-1,:,:] + image[h-1,:,:] - image[0,:,:]

  # Compute multiplier
  x = arange (0., 1., 1./w)
  y = arange (0., 1., 1./h)
  xx,yy = meshgrid (x,y)
  multi = 4 - 2.*cos(2*pi*xx) - 2.*cos(2*pi*yy)
  multi[0,0] = 1.

  # Compute DFT of boundary image
  sh = fftn (v, axes=(0, 1))

  # Multiply by inverse of multiplier
  sh = sh / multi.reshape((h,w,1))
  sh[0,0,:] = zeros ((d))

  # Then, compute s as the iDFT of sh
  smooth = real (ifftn (sh, axes=(0, 1)))
  periodic = image - smooth

  return harmonize(periodic),harmonize(smooth)

    def source_terms(self, mara, retphi=False):
        from numpy.fft import fftfreq, fftn, ifftn
        ng = mara.number_guard_zones()
        G = self.G
        L = 1.0
        Nx, Ny, Nz = mara.fluid.shape
        Nx -= 2*ng
        Ny -= 2*ng
        Nz -= 2*ng
        P = mara.fluid.primitive[ng:-ng,ng:-ng,ng:-ng]
        rho = P[...,0]
        vx = P[...,2]
        vy = P[...,3]
        vz = P[...,4]

        K = [fftfreq(Nx)[:,np.newaxis,np.newaxis] * (2*np.pi*Nx/L),
             fftfreq(Ny)[np.newaxis,:,np.newaxis] * (2*np.pi*Ny/L),
             fftfreq(Nz)[np.newaxis,np.newaxis,:] * (2*np.pi*Nz/L)]
        delsq = -(K[0]**2 + K[1]**2 + K[2]**2)
        delsq[0,0,0] = 1.0 # prevent division by 0

        rhohat = fftn(rho)
        phihat = (4*np.pi*G) * rhohat / delsq
        fx = -ifftn(1.j * K[0] * phihat).real
        fy = -ifftn(1.j * K[1] * phihat).real
        fz = -ifftn(1.j * K[2] * phihat).real

        S = np.zeros(mara.fluid.shape + (5,))
        S[ng:-ng,ng:-ng,ng:-ng,0] = 0.0
        S[ng:-ng,ng:-ng,ng:-ng,1] = rho * (fx*vx + fy*vy + fz*vz)
        S[ng:-ng,ng:-ng,ng:-ng,2] = rho * fx
        S[ng:-ng,ng:-ng,ng:-ng,3] = rho * fy
        S[ng:-ng,ng:-ng,ng:-ng,4] = rho * fz
        return (S, ifftn(phihat).real) if retphi else S

    def steepest_descent_update(self, sdi, IWxp, template):
        # compute error image
        # error_img:  height  x  width  x  n_channels
        error_img = IWxp.pixels - template.pixels

        # compute FFT error image
        # fft_error_img:  height  x  width  x  n_channels
        fft_axes = range(IWxp.n_dims)
        fft_error_img = fftshift(fftn(error_img, axes=fft_axes),
                                 axes=fft_axes)

        # reshape FFT error image
        # fft_error_img:  (height x width)  x  n_channels
        fft_error_img = np.reshape(fft_error_img, (-1, IWxp.n_channels))

        # compute filtered steepest descent images
        # _filter_bank:        (height x width)  x
        # fft_error_img:       (height x width)  x  n_channels
        # filtered_error_img:  (height x width)  x  n_channels
        filtered_error_img = (self._filter_bank[..., None] * fft_error_img)

        # reshape _error_img
        # _error_img:  (height x width x n_channels)
        self._error_img = filtered_error_img.flatten()

        # compute steepest descent update
        # sdi:         (height x width x n_channels)  x  n_parameters
        # _error_img:  (height x width x n_channels)
        # sdu:             n_parameters
        return sdi.T.dot(np.conjugate(self._error_img))