code: def FDMT_initialization(datain, f_min, f_max, maxDT, dataType):
      
      # 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
docs: 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)