Input: datain - visibilities of (nint, nbl, nchan, npol)
        f_min,f_max - are the base-band begin and end frequencies.
            The frequencies can be entered in both MHz and GHz, units are factored out in all uses.
        maxDT - the maximal delay (in time bins) of the maximal dispersion.
            Appears in the paper as N_{\Delta}
            A typical input is maxDT = N_f
        dataType - To naively use FFT, one must use floating point types.
            Due to casting, use either complex64 or complex128.
    Output: dataout, 3d array, with dimensions [nint, N_d0, nbl, nchan, npol]
            where N_d0 is the maximal number of bins the dispersion curve travels at one frequency bin
    
    For details, see algorithm 1 in Zackay & Ofek (2014)

def FDMT_initialization(datain, f_min, f_max, maxDT, dataType):
    """
    Input: datain - visibilities of (nint, nbl, nchan, npol)
        f_min,f_max - are the base-band begin and end frequencies.
            The frequencies can be entered in both MHz and GHz, units are factored out in all uses.
        maxDT - the maximal delay (in time bins) of the maximal dispersion.
            Appears in the paper as N_{\Delta}
            A typical input is maxDT = N_f
        dataType - To naively use FFT, one must use floating point types.
            Due to casting, use either complex64 or complex128.
    Output: dataout, 3d array, with dimensions [nint, N_d0, nbl, nchan, npol]
            where N_d0 is the maximal number of bins the dispersion curve travels at one frequency bin
    
    For details, see algorithm 1 in Zackay & Ofek (2014)
    """
    # Data initialization is done prior to the first FDMT iteration
    # See Equations 17 and 19 in Zackay & Ofek (2014)

    [nint, nbl, nchan, npol] = datain.shape

    deltaF = (f_max - f_min)/float(nchan)
    deltaT = int(np.ceil((maxDT-1) *(1./f_min**2 - 1./(f_min + deltaF)**2) / (1./f_min**2 - 1./f_max**2)))

    dataout = np.zeros([nint, deltaT+1, nbl, nchan, npol], dataType)
    dataout[:,0,:,:,:] = datain
    
    for i_dT in xrange(1, deltaT+1):
        dataout[i_dT:, i_dT, :, :, :] = dataout[i_dT:, i_dT-1, :, :, :] + datain[:-i_dT]
    return dataout