This function determines the convolution of two inputs using the FFT. Contains an implementation for both CPU
    and GPU.

    INPUTS:
    in1             (no default):           Array containing one set of data, possibly an image.
    in2             (no default):           Gpuarray containing the FFT of the PSF.
    conv_device     (default = "cpu"):      Parameter which allows specification of "cpu" or "gpu".
    conv_mode       (default = "linear"):   Mode specifier for the convolution - "linear" or "circular".

def fft_convolve(in1, in2, conv_device="cpu", conv_mode="linear", store_on_gpu=False):
    """
    This function determines the convolution of two inputs using the FFT. Contains an implementation for both CPU
    and GPU.

    INPUTS:
    in1             (no default):           Array containing one set of data, possibly an image.
    in2             (no default):           Gpuarray containing the FFT of the PSF.
    conv_device     (default = "cpu"):      Parameter which allows specification of "cpu" or "gpu".
    conv_mode       (default = "linear"):   Mode specifier for the convolution - "linear" or "circular".
    """

    # NOTE: Circular convolution assumes a periodic repetition of the input. This can cause edge effects. Linear
    # convolution pads the input with zeros to avoid this problem but is consequently heavier on computation and
    # memory.

    if conv_device=='gpu':

        if conv_mode=="linear":
            fft_in1 = pad_array(in1)
            fft_in1 = gpu_r2c_fft(fft_in1, store_on_gpu=True)
            fft_in2 = in2

            conv_in1_in2 = fft_in1*fft_in2

            conv_in1_in2 = contiguous_slice(fft_shift(gpu_c2r_ifft(conv_in1_in2, is_gpuarray=True, store_on_gpu=True)))

            if store_on_gpu:
                return conv_in1_in2
            else:
                return conv_in1_in2.get()

        elif conv_mode=="circular":
            fft_in1 = gpu_r2c_fft(in1, store_on_gpu=True)
            fft_in2 = in2

            conv_in1_in2 = fft_in1*fft_in2

            conv_in1_in2 = fft_shift(gpu_c2r_ifft(conv_in1_in2, is_gpuarray=True, store_on_gpu=True))

            if store_on_gpu:
                return conv_in1_in2
            else:
                return conv_in1_in2.get()
    else:

        if conv_mode=="linear":
            fft_in1 = pad_array(in1)
            fft_in2 = in2

            out1_slice = tuple(slice(0.5*sz,1.5*sz) for sz in in1.shape)

            return np.require(np.fft.fftshift(np.fft.irfft2(fft_in2*np.fft.rfft2(fft_in1)))[out1_slice], np.float32, 'C')

        elif conv_mode=="circular":
            return np.fft.fftshift(np.fft.irfft2(in2*np.fft.rfft2(in1)))