import numpy as np

from fft_routines import CONVOLVE_3D,COMPUTE_WINDOW_FUNCTION_REAL_SPACE,\
CONVOLVE_3D_WITH_FFT_INPUT,COMPUTE_FFT_3D_NO_NORM_NO_SHIFT

def spherical_averaging_main():

    print "Computing FFT of 3D Signal Field ...."
    L_box          = 160.0
    N_Grid         = 256
    inp_path       = "/path/to/3d/data/cub/e"
    inp_file       = "3d_data_cube_file_name.npy"
    signal_arr     = np.load(inp_path + inp_file)
    sigmal_arr_fft = COMPUTE_FFT_3D_NO_NORM_NO_SHIFT(signal_arr)

    # R_sphere is grid size over which spherical top hat averaging needs to be done.
    for R_sphere in [10,20,30]:
        inp_dict_fft   = {"L_box":L_box,"N_Grid":N_Grid,"R_sphere":R_sphere}

        print "Computing FFT of Spherical Top-hat Filter R_sphere =",R_sphere

        # Compute window function in Real Space
        WR_arr_3d,npix  = COMPUTE_WINDOW_FUNCTION_REAL_SPACE(inp_dict_fft,flag_window = "top_hat")

        # Compute FFT of the window function (Fourier Space)
        WR_arr_3d_fft   = COMPUTE_FFT_3D_NO_NORM_NO_SHIFT(WR_arr_3d)

        # Perform convolution using FFT of signal and window function
        # And inverse fourier transform it back to real space.
        # signal_arr_smth is spherically averaged values at each pixel.

        signal_arr_smth = CONVOLVE_3D_WITH_FFT_INPUT(signal_arr_fft,WR_arr_3d_fft)

        # Do you your halo calculations here ....


spherical_averaging_main()