program main
  use sep
  use wei_global_defi_mod
  use wei_global_calc_mod
  use wei_io_mod
  use wei_vel_mod
  use wei_ssf_mod
  use wei_cig_mod
  use myfft_mod

  implicit none

  real, allocatable, dimension(:,:,:)    :: data, wfld_d_data
  complex, allocatable, dimension(:,:,:) :: cdata

  real, allocatable, dimension(:)        :: sfun
  complex, allocatable, dimension(:)     :: csfun, ctrace

  complex, allocatable, dimension(:,:,:)   ,target :: wfld_d, wfld_u
  real   , allocatable, dimension(:,:,:)   ,target :: model
  complex, dimension(:,:,:)  , pointer :: pwfld_d, pwfld_u
  real   , dimension(:,:), pointer :: pmodel

  integer    :: nz, iz, stat, izs
  real       :: vel 
  type(axis) :: agx, agy, aw, asx, asy
  type(axis) :: az, ahx, ahy, akx, aky
  
  integer, parameter  :: opsign = 1
  real   , parameter  :: PI = 3.1415926
  logical, parameter  :: T = .true., F = .false.

  integer             :: iw, igx, igy, nt, it, nw1, nw2, nw3, nw4, nw
  real                :: oz, ow, ot, dz, dt, dw, f1, f2, f3, f4, w1, w2, w3, w4
  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()

  call from_aux("rec", "n1", nt)
  call from_aux("rec", "d1", dt)
  call from_aux("rec", "o1", ot)

  call from_aux("rec", "n2", agx%n)
  call from_aux("rec", "d2", agx%d)
  call from_aux("rec", "o2", agx%o)

  call from_aux("rec", "n3", agy%n)
  call from_aux("rec", "d3", agy%d)
  call from_aux("rec", "o3", agy%o)

  allocate(data(nt, agx%n, agy%n))
  allocate(sfun(nt))

  write(0,*)"read data"
  !call sep_read(data, "rec", esize=4)
  !call sep_read(sfun, "sou", esize=4)

  call sreed("rec", data, 4*nt*agx%n*agy%n)
  call sreed("sou", sfun, 4*nt)


  nz = 101; dz=10; oz=0
 
  f1 = 3; f2 = 8; f3=30; f4=45

  dw = 2*PI/(nt*dt)

  w1=2*PI*f1; nw1=w1/dw+1.5
  w2=2*PI*f2; nw2=w2/dw+1.5
  w3=2*PI*f3; nw3=w3/dw+1.5
  w4=2*PI*f4; nw4=w4/dw+1.5

  
  nw=nw4-nw1+1
  ow=(nw1-1)*dw

   
  write(0,*)"nw1, nw2, nw3, nw4, nw"
  write(0,*) nw1, nw2, nw3, nw4, nw
  write(0,*)"nt=", nt
  aw%n = nw; aw%d = dw; aw%o = ow
  az%n = nz; az%d = dz; az%o = oz
  allocate(cdata(agx%n,agy%n,aw%n))
  allocate(csfun(aw%n), ctrace(nt))
  
  cdata = 0.
  csfun = 0.
  ctrace = 0.

  write(0,*)"1D FFT for the source"
  ! 1D FFT for the source function
  ctrace = sfun
  call FFT1DT(ctrace, 1)
  call hammingw(ctrace, nw1, nw2, nw3, nw4)
  do iw=1, nw
     csfun(iw) = ctrace(iw+nw1-1)
  end do


  ctrace = 0.
  do iw=1, nw
     ctrace(iw+nw1-1) = csfun(iw)
  end do
  call FFT1DT(ctrace, -1)

  sfun = real(ctrace)

  call to_history("n1", nt, "fft")
  call sep_write(sfun, "fft")




  write(0,*)"1D FFT for the receiver"
  ! 1D FFT for the receiver wavefield
  ctrace = 0.
  do igy=1 , agy%n
     do igx=1, agx%n
        do it=1, nt
           ctrace(it) = data(it, igx, igy)
        end do
        call FFT1DT(ctrace, 1)
        call hammingw(ctrace, nw1, nw2, nw3, nw4)
        do iw=1, nw
           cdata(igx,igy,iw) = ctrace(iw+nw1-1)
        end do
        ctrace = 0.
     end do
  end do
  write(0,*)"deallocate ..."
  deallocate(data, sfun)

  write(0,*)"get grid"
  call get_computation_grid(asx, asy, agx, agy, az)

  write(0,*)"grid%x=", grid%x
  write(0,*)"grid%y=", grid%y
  write(0,*)"grid%z=", grid%z
   
  allocate(wfld_u(grid%x%n,grid%y%n,aw%n))
  allocate(wfld_d(grid%x%n,grid%y%n,aw%n))
  allocate(model(grid%x%n,grid%y%n,grid%z%n)) 

  wfld_d = 0.
  wfld_u = 0.
  model  = 0.

  wfld_d(grid%iosx, grid%iosy, :) = csfun
  wfld_u(grid%iogx:grid%iogx+size(wfld_u, 1)-1, grid%iogy:grid%iogy+size(wfld_u, 2)-1, :) = cdata

  deallocate(cdata, csfun)

  pwfld_d => wfld_d
  pwfld_u => wfld_u
  !pmodel  => model

  ahx%n=1; ahx%o=0; ahx%d=10;
  ahy%n=1; ahy%o=0; ahy%d=10; 
  
  akx%n=grid%x%n
  akx%d=2.0*PI/(akx%n*grid%x%d)
  !akx%o=-PI/grid%x%d
  akx%o= 0.

  aky%n=grid%y%n
  aky%d=2.0*PI/(aky%n*grid%y%d)
  !aky%o=-PI/grid%y%d
  aky%o= 0.



  vref = 1500 

  call auxputch("n1", "i", grid%x%n, "image")
  call auxputch("d1", "r", grid%x%d, "image")
  call auxputch("o1", "r", grid%x%o, "image")

  call auxputch("n2", "i", grid%y%n, "image")
  call auxputch("d2", "r", grid%y%d, "image")
  call auxputch("o2", "r", grid%y%o, "image")

  call auxputch("n3", "i", grid%z%n, "image")
  call auxputch("d3", "r", grid%z%d, "image")
  call auxputch("o3", "r", grid%z%o, "image")

  call auxputch("esize", "i", 4, "image")
  write(0,*)"ssf init"
  call wei_ssf_init(az, aw, akx, aky)

  call from_param("izs", izs, 50)

  do iz=1, az%n 
     write(0,*)"iz=", iz
     call sep_read(vmod, "vmod")
     call wei_get_vref_slice( F )
     stat = wei_ssf_down( F, pwfld_d )
     stat = wei_ssf_down( T, pwfld_u )

     if (iz.eq.izs) then
         write(0,*)"registering wfld_d"
         allocate(wfld_d_data(nt, grid%x%n, grid%y%n))
         do igy=1, grid%y%n
            do igx=1, grid%x%n
               ctrace = 0.
               do iw=1, nw
                  ctrace(iw+nw1-1) = wfld_u(igx, igy, iw)
               end do
               call FFT1DT(ctrace, -1)
               do it=1, nt
                  wfld_d_data(it, igx, igy) = real(ctrace(it))
               end do
            end do
         end do
         call auxputch("n1", "i", nt, "wfld_d")
         call auxputch("d1", "r", dt, "wfld_d")
         call auxputch("o1", "r", ot, "wfld_d")

         call auxputch("n2", "i", grid%x%n, "wfld_d")
         call auxputch("d2", "r", grid%x%d, "wfld_d")
         call auxputch("o2", "r", grid%x%o, "wfld_d")

         call auxputch("n3", "i", grid%y%n, "wfld_d")
         call auxputch("d3", "r", grid%y%d, "wfld_d")
         call auxputch("o3", "r", grid%y%o, "wfld_d")
         call srite("wfld_d",wfld_d_data,4*nt*agx%n*agy%n)
         deallocate(wfld_d_data)
         write(0,*)"registering finished"
         
     end if
     do iw=1, nw
        model(:,:,iz) = model(:,:,iz) + real(conjg(pwfld_d(:,:,iw))*pwfld_u(:,:,iw))
     end do

  end do
 
  call sep_write(model, "image")

end program









