program main
  use wei_myfft_mod
  use wei_hamming_mod
  use sep
  implicit none
  
  

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

  real, parameter :: PI = 3.14159265359


  call initpar()

  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
  if (adj) then
       write(0,*)"Forward FFT..."
       call from_history("n1",  nt); call from_history("d1",  dt); call from_history("o1",  ot)
       call from_history("n2", ngx); call from_history("d2", dgx); call from_history("o2", ogx)
       call from_history("n3", ngy); call from_history("d3", dgy); call from_history("o3", ogy)
       call from_history("n4", nsx); call from_history("d4", dsx); call from_history("o4", osx)
       call from_history("n5", nsy); call from_history("d5", dsy); call from_history("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))
       call to_history("n1",  nw); call to_history("d1",  dw); call to_history("o1",  ow)
       call to_history("n2", ngx); call to_history("d2", dgx); call to_history("o2", ogx)
       call to_history("n3", ngy); call to_history("d3", dgy); call to_history("o3", ogy)
       call to_history("n4", nsx); call to_history("d4", dsx); call to_history("o4", osx)
       call to_history("n5", nsy); call to_history("d5", dsy); call to_history("o5", osy)
       call to_history("esize", 8); call to_history("label1", "frequency")
       do isy=1, nsy
          do isx=1, nsx
             call sreed("in", data, 8*size(data))
             do igy=1, ngy
                do igx=1, ngx
                   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 srite("out", modl, 8*size(modl))
          end do
       end do

  !!.. in case of demigration, inverse FFT, from w to t
  else
       write(0,*)"Inverse FFT..."
       call from_history("n1",  nw); call from_history("d1",  dw); call from_history("o1",  ow)
       call from_history("n2", ngx); call from_history("d2", dgx); call from_history("o2", ogx)
       call from_history("n3", ngy); call from_history("d3", dgy); call from_history("o3", ogy)
       call from_history("n4", nsx); call from_history("d4", dsx); call from_history("o4", osx)
       call from_history("n5", nsy); call from_history("d5", dsy); call from_history("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 to_history("n1",  nt); call to_history("d1",  dt); call to_history("o1",  ot)
       call to_history("n2", ngx); call to_history("d2", dgx); call to_history("o2", ogx)
       call to_history("n3", ngy); call to_history("d3", dgy); call to_history("o3", ogy)
       call to_history("n4", nsx); call to_history("d4", dsx); call to_history("o4", osx)
       call to_history("n5", nsy); call to_history("d5", dsy); call to_history("o5", osy)
       call to_history("esize", 8); call to_history("label1", "time")
       do isy=1, nsy
          do isx=1, nsx
             call sreed("in", modl, 8*size(modl))
             do igy=1, ngy
                do igx=1, ngx
                   do iw=1, nw
                      ctrace(nw1+iw-1) = modl(iw, igx, igy)
                      !ctrace(nw1+iw-1+nt/2) = modl(iw, igx, igy)
                   end do
                   call hammingw(ctrace, nw1, nw2, nw3, nw4)
                   call FFT1DT(ctrace, -1)
                   data(:, igx, igy) = data(:, igx, igy) + ctrace
                   ctrace = 0.
                end do 
             end do
             call srite("out", data, 8*size(data))
          end do
       end do
  end if


end program









