program main
  use wei_global_defi_mod
  use wei_global_calc_mod
  use wei_io_mod
  use wei_shotmig_kernel_mod
  use wei_vel_mod
  use wei_ssf_mod
  use wei_cig_mod
  use myfft_mod
  use sep
  use filter_mod
  implicit none

  integer  :: iws, izb, it, nt, stat
  integer  :: nw1, nw2, nw3, nw4, igx, igy
  real     :: f1, f2, f3, f4, w1, w2, w3, w4, ot, dt
  real, allocatable, dimension(:)        :: sfun
  real, allocatable, dimension(:,:,:)    :: data
  complex, allocatable, dimension(:)     :: ctrace
  

  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", dgrid%agx%n)
  call from_aux("rec", "d2", dgrid%agx%d)
  call from_aux("rec", "o2", dgrid%agx%o)

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

  call from_aux("rec", "n4", dgrid%asx%n)
  call from_aux("rec", "d4", dgrid%asx%d)
  call from_aux("rec", "o4", dgrid%asx%o)

  call from_aux("rec", "n5", dgrid%asy%n)
  call from_aux("rec", "d5", dgrid%asy%d)
  call from_aux("rec", "o5", dgrid%asy%o)

  call from_aux("vel", "n1", vgrid%ax%n)
  call from_aux("vel", "d1", vgrid%ax%d)
  call from_aux("vel", "o1", vgrid%ax%o)

  call from_aux("vel", "n2", vgrid%ay%n)
  call from_aux("vel", "d2", vgrid%ay%d)
  call from_aux("vel", "o2", vgrid%ay%o)

  call from_aux("vel", "n3", vgrid%az%n)
  call from_aux("vel", "d3", vgrid%az%d)
  call from_aux("vel", "o3", vgrid%az%o)

  nref = 1
  lloyd = F

  !igrid%ax = vgrid%ax
  !igrid%ay = vgrid%ay
  !igrid%az = vgrid%az
 
  igrid%image_xmin = -1000
  igrid%image_xmax =  1000
  igrid%image_ymin = -1000
  igrid%image_ymax =  1000
  igrid%image_dx   = vgrid%ax%d
  igrid%image_dy   = vgrid%ay%d

  igrid%ahx%n = 1
  igrid%ahy%n = 1

  azs%n = vgrid%az%n

  f1 = 3; f2 = 5; f3=8; f4=10

  aws%d = 2*PI/(nt*dt)

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

  aws%n=nw4-nw1+1
  aws%o=(nw1-1)*aws%d
  
  aw = aws

  call wei_global_calc()

  write(0,*)"igridx", igrid%ax%n, igrid%ax%o, igrid%ax%d
  write(0,*)"igridy", igrid%ay%n, igrid%ay%o, igrid%ay%d
  write(0,*)"ogridx", ogrid%ax%n, ogrid%ax%o, ogrid%ax%d
  write(0,*)"ogridy", ogrid%ay%n, ogrid%ay%o, ogrid%ay%d
  write(0,*)"pgridx", pgrid%ax%n, pgrid%ax%o, pgrid%ax%d
  write(0,*)"pgridy", pgrid%ay%n, pgrid%ay%o, pgrid%ay%d

  if(igrid%ax%n>1.and.ogrid%ax%n == 0) STOP
  if(igrid%ay%n>1.and.ogrid%ay%n == 0) STOP 

  call wei_image_history()
  call wei_io_init()
  call wei_array_alloc()

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



  write(0,*)"1D FFT for the source"
  ctrace = sfun
  call FFT1DT(ctrace, 1)
  call hammingw(ctrace, nw1, nw2, nw3, nw4)
  do iws=1, aws%n
     csou(iws) = ctrace(iws+nw1-1)
  end do

  write(0,*)"1D FFT for the receiver"
  ctrace = 0.
  do igy=1 , dgrid%agy%n
     do igx=1, dgrid%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 iws=1, aws%n
           crec(igx,igy,iws) = ctrace(iws+nw1-1)
        end do
        ctrace = 0.
     end do
  end do
  write(0,*)"deallocate ..."
  deallocate(data, sfun)


  call wei_wfld_init()

  call shotmig_kernel_init(VLin = wei_vel_init, &
                           FXin = wei_ssf_init, &
                           IGin = wei_cig_init)
  do izb=1, nzb

     write(0,*)"izb=",izb
     call wei_io_image(F, izb)
     stat = shotmig_kernel_lop(VLop = wei_vel_lop, &
                               FXop = wei_ssf_lop, &
                               IGop = wei_cig_lop)
 
     call wei_io_image(T, izb)

  end do

  call wei_array_dealloc()

end program









