program main
  use sep
  use base_struct_mod
  use base_fftwv3_mod
  use base_fftsft_mod
  use wem_wnbmap_mod
  use mva_resmig_mod
  implicit none

  integer*8                              :: plan_cdat_3dxF
  integer*8                              :: plan_cmod_3dxI
  integer                                :: stat, nkx, nky, nkz
  integer                                :: nxpad_beg, nxpad_end
  integer                                :: nypad_beg, nypad_end
  integer                                :: nzpad_beg, nzpad_end
  integer                                :: ikx, iky, irho
  type(axis)                             :: ax, ay, az, akz, arho
  real   , dimension(:)    , allocatable :: kx, ky
  real   , dimension(:,:,:), allocatable :: rdat
  complex, dimension(:,:,:), allocatable :: cdat, cmod
  logical                                :: adj, add
  real                                   :: PI, rho, eps

  PI = 3.14159265359

  call sep_init()

  call from_history("n1", az%n); call from_history("d1", az%d); call from_history("o1", az%o)
  call from_history("n2", ax%n); call from_history("d2", ax%d); call from_history("o2", ax%o)
  call from_history("n3", ay%n); call from_history("d3", ay%d); call from_history("o3", ay%o)

write(0,*)"nx, ny, nz", ax%n, ay%n, az%n

  call from_param("orho", arho%o, 0.8   )
  call from_param("drho", arho%d, 0.0001)
  call from_param("nrho", arho%n, 1     )
  call from_param("eps" , eps   , 0.1   )

  call to_history("n4", arho%n); call to_history("d4", arho%d); call to_history("o4", arho%o)

  call from_param("nxpad_beg", nxpad_beg, 0); call from_param("nxpad_end", nxpad_end, 0)
  call from_param("nypad_beg", nypad_end, 0); call from_param("nypad_end", nypad_end, 0)
  call from_param("nzpad_beg", nzpad_beg, 0); call from_param("nzpad_end", nzpad_end, 0)

  nkx = ax%n + nxpad_beg + nxpad_end
  nky = ay%n + nypad_beg + nypad_end
  nkz = az%n + nzpad_beg + nzpad_end

!  call to_history("n1", nkz); call to_history("n2", nkx); call to_history("n3", nky)

  allocate(kx(nkx), ky(nky))

  allocate(rdat(az%n,ax%n,ay%n))
  allocate(cdat( nkz, nkx, nky))
  allocate(cmod( nkz, nkx, nky))
 
  call sep_read(rdat)

  write(0,*)"generating plans..."

  call sfftw_dft_3dx_init(plan_cdat_3dxF, cdat, nkz, nkx, nky, +1)
  call sfftw_dft_3dx_init(plan_cmod_3dxI, cmod, nkz, nkx, nky, -1)

  cdat(nzpad_beg+1:nzpad_beg+az%n,nxpad_beg+1:nxpad_beg+ax%n,nypad_beg+1:nypad_beg+ay%n) = rdat

  write(0,*)"FFT..."
  call sfftw_dft_3dx_oper(plan_cdat_3dxF, cdat, nkz, nkx, nky)

  call k_map(kx, nkx, ax%d)
  call k_map(ky, nky, ay%d)

  akz%o = 0.
  akz%d = 2*PI/(nkz*az%d)
  akz%n = nkz/2+1
  cdat(akz%n+1:,:,:) = 0.

  write(0,*)"residual migration..."
  stat = mva_resmig_init(kx, ky, akz, eps)

  do irho=1, arho%n
     rho = arho%o + (irho-1)*arho%d
     write(0,*)"irho=", irho, "rho=", rho
     stat = mva_resmig_oper(cdat(1:akz%n,:,:), cmod(1:akz%n,:,:), rho, 0., 0.)
     call sfftw_dft_3dx_oper(plan_cmod_3dxI, cmod, nkz, nkx, nky) 
     call sep_write(real(cmod(nzpad_beg+1:nzpad_beg+az%n,nxpad_beg+1:nxpad_beg+ax%n,nypad_beg+1:nypad_beg+ay%n)))
  end do
  stat = mva_resmig_clos()


end program
