module wei_transf_mod
  use wei_myfft_mod
  use wei_hamming_mod
  use sep

  integer, private  :: nw1, nw2, nw3, nw4, nt_default
  integer, private  :: nt, nw, ngx, ngy, nsx, nsy
  real   , private  :: ot, ow, ogx, ogy, osx, osy
  real   , private  :: dt, dw, dgx, dgy, dsx, dsy
  real   , private  :: f1, f2, f3, f4, w1, w2, w3, w4
  complex, dimension(:)    , pointer, private :: ctrace
  complex, dimension(:,:,:), pointer, private :: modl
  real   , dimension(:,:,:), pointer, private :: data
  logical                           , private :: adj

  real, parameter :: PI = 3.14159265359
  character(len=128)                , private :: data_="data"
  character(len=128)                , private :: modl_="modl" 
contains
!--------------------------------------------------------------------------------------  
  integer function wei_transf_init() result(stat)
    call wei_transf_param_init()
    stat = 0
  end function
!--------------------------------------------------------------------------------------
  integer function wei_transf_clos() result(stat)
    deallocate(modl)
    deallocate(data)
    deallocate(ctrace)
    stat = 0
  end function
!--------------------------------------------------------------------------------------
  integer function wei_transf_oper() result(stat)
    integer      :: isx, isy, igx, igy, it, iw

write(0,*)"nsx, nsy, ngx, ngy, nw, nt", nsx, nsy, ngx, ngy, nw, nt

    if (adj) then
       do isy=1, nsy
       do isx=1, nsx
          call wei_transf_io_data("reed", data, isx, isy)
          modl = 0.
          do igy=1, ngy
          do igx=1, ngx
             ctrace = ctrace + data(:, igx, igy)
             call FFT1DT(ctrace, 1)
             call hammingw(ctrace, nw1, nw2, nw3, nw4)
             do iw=1, nw
                modl(iw, igx, igy) = modl(iw, igx, igy) + ctrace(iw+nw1-1)
             end do
             ctrace = 0.
          end do
          end do
          call wei_transf_io_modl("rite", modl, isx, isy)
       end do
       end do

    else
       do isy=1, nsy
       do isx=1, nsx
          call wei_transf_io_modl("reed", modl, isx, isy)
          data = 0.
          do igy=1, ngy
          do igx=1, ngx
             do iw=1, nw
                ctrace(nw1+iw-1) = modl(iw, igx, igy)
             end do
             call hammingw(ctrace, nw1, nw2, nw3, nw4)
             call FFT1DT(ctrace, -1)
             data(:, igx, igy) = data(:, igx, igy) + real(ctrace)
             ctrace = 0.
          end do
          end do
          call wei_transf_io_data("rite", data, isx, isy) 
       end do
       end do
    end if
    stat = 0
  end function
!--------------------------------------------------------------------------------------  
  subroutine wei_transf_io_data(flag, data, isx, isy)
    character(len=4)                :: flag
    real, dimension(:,:,:), pointer :: data
    integer, dimension(5)           :: n_cube, j_wind, n_wind, f_wind
    integer                         :: n_dim_cube, err
    integer                         :: isx, isy

    n_dim_cube = 5
    n_cube = (/nt, ngx, ngy, nsx, nsy/)
    j_wind = (/1, 1, 1, 1, 1/)
    n_wind = (/nt, ngx, ngy, 1, 1/)
    f_wind = (/0, 0, 0, isx-1, isy-1/)
    if (flag == "reed") then
       data = 0.
       err = sreed_window(data_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,data)
    elseif (flag == "rite") then
       err = srite_window(data_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,data)
    else
       call seperr("error: wei_transf_io_data")
    end if
  end subroutine
!------------------------------------------------------------------------------------------
  subroutine wei_transf_io_modl(flag, modl, isx, isy)
    character(len=4)                   :: flag
    complex, dimension(:,:,:), pointer :: modl
    integer, dimension(5)              :: n_cube, j_wind, n_wind, f_wind
    integer                            :: n_dim_cube, err
    integer                            :: isx, isy

    n_dim_cube = 5
    n_cube = (/nw, ngx, ngy, nsx, nsy/)
    j_wind = (/1, 1, 1, 1, 1/)
    n_wind = (/nw, ngx, ngy, 1, 1/)
    f_wind = (/0, 0, 0, isx-1, isy-1/)
    if (flag == "reed") then
       modl = 0.
       err = sreed_window(modl_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,modl)
    elseif (flag == "rite") then
       err = srite_window(modl_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,modl)
    else
       call seperr("error: wei_transf_io_modl")
    end if
  end subroutine
!--------------------------------------------------------------------------------------------
  subroutine wei_transf_param_init()

    call from_param("f1", f1,  5.)
    call from_param("f2", f2, 10.)
    call from_param("f3", f3, 35.)
    call from_param("f4", f4, 40.)

    call from_param("adj", adj, .true.)

    w1=2*PI*f1
    w2=2*PI*f2
    w3=2*PI*f3
    w4=2*PI*f4
    !!.. in case of migration, forward FFT, from t to w
    !!.. dat => time domain
    !!.. mod => freq domain
    call auxinout(data_)
    call auxinout(modl_)
    if (adj) then
       write(0,*)"Forward FFT..."
       call from_aux(data_, "n1",  nt); call from_aux(data_, "d1",  dt); call from_aux(data_, "o1",  ot)
       call from_aux(data_, "n2", ngx); call from_aux(data_, "d2", dgx); call from_aux(data_, "o2", ogx)
       call from_aux(data_, "n3", ngy); call from_aux(data_, "d3", dgy); call from_aux(data_, "o3", ogy)
       call from_aux(data_, "n4", nsx); call from_aux(data_, "d4", dsx); call from_aux(data_, "o4", osx)
       call from_aux(data_, "n5", nsy); call from_aux(data_, "d5", dsy); call from_aux(data_, "o5", osy)

       dw = 2*PI/(nt*dt)
       nw1=w1/dw+0.5
       nw2=w2/dw+0.5
       nw3=w3/dw+0.5
       nw4=w4/dw+0.5
       ow = (nw1-1)*dw
       nw = nw4 - nw1 + 1
       allocate(data(nt, ngx, ngy))
       allocate(modl(nw, ngx, ngy))
       allocate(ctrace(nt))
       ctrace = 0.

       call auxputch("n1", "i",  nw, modl_); call auxputch("d1", "r",  dw, modl_); call auxputch("o1", "r",  ow, modl_)
       call auxputch("n2", "i", ngx, modl_); call auxputch("d2", "r", dgx, modl_); call auxputch("o2", "r", ogx, modl_)
       call auxputch("n3", "i", ngy, modl_); call auxputch("d3", "r", dgy, modl_); call auxputch("o3", "r", ogy, modl_)
       call auxputch("n4", "i", nsx, modl_); call auxputch("d4", "r", dsx, modl_); call auxputch("o4", "r", osx, modl_)
       call auxputch("n5", "i", nsy, modl_); call auxputch("d5", "r", dsy, modl_); call auxputch("o5", "r", osy, modl_)
       call auxputch("esize", "i", 8, modl_); call auxputch("label1", "s", "frequency", modl_)
    else
       write(0,*)"Inverse FFT..."
       call from_aux(modl_, "n1",  nw); call from_aux(modl_, "d1",  dw); call from_aux(modl_, "o1",  ow)
       call from_aux(modl_, "n2", ngx); call from_aux(modl_, "d2", dgx); call from_aux(modl_, "o2", ogx)
       call from_aux(modl_, "n3", ngy); call from_aux(modl_, "d3", dgy); call from_aux(modl_, "o3", ogy)
       call from_aux(modl_, "n4", nsx); call from_aux(modl_, "d4", dsx); call from_aux(modl_, "o4", osx)
       call from_aux(modl_, "n5", nsy); call from_aux(modl_, "d5", dsy); call from_aux(modl_, "o5", osy)

       nw1=w1/dw+0.5
       nw2=w2/dw+0.5
       nw3=w3/dw+0.5
       nw4=w4/dw+0.5

       call from_param("dt", dt, 0.004)
       nt_default = 2*PI / (dt*dw)
       call from_param("nt", nt, nt_default)
       ot = 0.

       allocate(data(nt, ngx, ngy))
       allocate(modl(nw, ngx, ngy))
       allocate(ctrace(nt))
       ctrace = 0.

       call auxputch("n1", "i",  nt, data_); call auxputch("d1", "r",  dt, data_); call auxputch("o1", "r",  ot, data_)
       call auxputch("n2", "i", ngx, data_); call auxputch("d2", "r", dgx, data_); call auxputch("o2", "r", ogx, data_)
       call auxputch("n3", "i", ngy, data_); call auxputch("d3", "r", dgy, data_); call auxputch("o3", "r", ogy, data_)
       call auxputch("n4", "i", nsx, data_); call auxputch("d4", "r", dsx, data_); call auxputch("o4", "r", osx, data_)
       call auxputch("n5", "i", nsy, data_); call auxputch("d5", "r", dsy, data_); call auxputch("o5", "r", osy, data_)
       call auxputch("esize", "i", 4, data_); call auxputch("label1", "s", "time", data_)
    end if
  end subroutine
end module








