program main
  use sep
  use delay_io_mod
  implicit none

  integer                            :: ngx, ngy, nws, nsx, nsy
  real                               :: dgx, dgy, dws, dsx, dsy
  real                               :: ogx, ogy, ows, osx, osy
  complex, dimension(:)    , pointer :: soufnc
  complex, dimension(:,:,:), pointer :: data
  real   , dimension(:,:)  , pointer :: phase_random
  integer                            :: isx, isy, iws, igx, igy
  integer                            :: index_loc1, index_loc2
  real                               :: gx, gy, sx, sy, ws, refx, refy
  real                               :: phase_delay, tdelay, loc1, loc2
  real                               :: random1, random2, linear1, linear2, wtx, angle, v0, p
  character(len=256)                 :: soufnc_, data_, phase_random_
  logical                            :: random, linear
  real, parameter                    :: PI = 3.14159265359

  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()

  soufnc_ = "soufnc"
  data_   = "output"
  phase_random_ = "phase_random"
  call auxin(soufnc_)
  call auxout(data_)
  if (random) call auxout(phase_random_)

  call from_aux(soufnc_, "n1", nws)
  call from_aux(soufnc_, "d1", dws)
  call from_aux(soufnc_, "o1", ows)

  call from_param("ngx", ngx, 1); call from_param("dgx", dgx, 1.); call from_param("ogx", ogx, 0.)
  call from_param("ngy", ngy, 1); call from_param("dgy", dgy, 1.); call from_param("ogy", ogy, 0.) 
  call from_param("nsx", nsx, 1); call from_param("dsx", dsx, 1.); call from_param("osx", osx, 0.)
  call from_param("nsy", nsy, 1); call from_param("dsy", dsy, 1.); call from_param("osy", osy, 0.)

  call from_param("loc1", loc1, 0.)
  call from_param("loc2", loc2, 0.)
  call from_param("tdelay", tdelay, 0.)
  call from_param("v0", v0, 2000.)
  call from_param("angle", angle, 0.)
  call from_param("refx", refx, 0.)
  call from_param("random", random, .false.)
  call from_param("linear", linear, .false.)

  p = sin(angle*PI/180.)/v0

  allocate(soufnc(nws))
  allocate(data(ngx,ngy,nws))
  allocate(phase_random(nws, 2))
  data = 0.
  phase_random = 0.
  call sep_read(soufnc, soufnc_)

  write(0,*)"ngx, ngy, nws, nsx, nsy", ngx, ngy, nws, nsx, nsy
  write(0,*)"ogx, ogy, ows, osx, osy", ogx, ogy, ows, osx, osy
  write(0,*)"dgx, dgy, dws, dsx, dsy", dgx, dgy, dws, dsx, dsy
  write(0,*)"loc", loc1, loc2

  index_loc1 = 1.5 + (loc1-ogx)/dgx
  index_loc2 = 1.5 + (loc2-ogx)/dgx

  write(0,*)"index_loc", index_loc1, index_loc2

  call random_seed()

  !! now only 2-D`
  do isy=1, nsy
     data = 0.
  do isx=1, nsx

     do iws=1, nws
        ws = ows + (iws-1)*dws
        phase_delay = ws*tdelay
!        call random_number(random1)
!        call random_number(random2)
!        random1 = random1 * 2*PI
!        random2 = random2 * 2*PI
!        data(index_loc1,isy,iws) = soufnc(iws)*cmplx(cos(phase_delay), sin(phase_delay))*cmplx(cos(random1),sin(random1))
!        data(index_loc2,isy,iws) = soufnc(iws)*cmplx(cos(phase_delay), sin(phase_delay))*cmplx(cos(random2),sin(random2))
         

        do igx=index_loc1-2, index_loc1+2
           gx = ogx + (igx-1)*dgx

           if (random) then
              call random_number(random1); random1 = random1*2*PI
           else
              random1 = 0.
           end if

           if (linear) then
              linear1 = ws*p*(gx-refx)
           else
              linear1 = 0. 
           end if
           
           wtx = 0.5 + 0.5 * cos(PI * real(igx-index_loc1)/3)
           data(igx,isy,iws) = wtx*soufnc(iws)*cexp(cmplx(0., phase_delay+random1+linear1))
        end do

        do igx=index_loc2-2, index_loc2+2
           gx = ogx + (igx-1)*dgx
           
           if (random) then
              call random_number(random2); random2 = random2*2*PI
           else
              random2 = 0.
           end if

           if (linear) then
              linear2 = ws*p*(gx-refx)
           else
              linear2 = 0.
           end if

           wtx = 0.5 + 0.5 * cos(PI * real(igx-index_loc2)/3)
           data(igx,isy,iws) = wtx*soufnc(iws)*cexp(cmplx(0., phase_delay+random2+linear2))
        end do


        phase_random(iws, 1) = random1
        phase_random(iws, 2) = random2
     end do

     if (random) call sep_write(phase_random, phase_random_)

     call iodata(data_, data, ngx, ngy, nws, nsx, nsy, isx, isy, "rite")

  end do
  end do

  if (random) then
     call auxputch("n1", "i", nws, phase_random_)
     call auxputch("n2", "i",   2, phase_random_)
  end if

  call auxputch("n1", "i", ngx, data_)
  call auxputch("d1", "r", dgx, data_)
  call auxputch("o1", "r", ogx, data_)

  call auxputch("n2", "i", ngy, data_)
  call auxputch("d2", "r", dgy, data_)
  call auxputch("o2", "r", ogy, data_)

  call auxputch("n3", "i", nws, data_)
  call auxputch("d3", "r", dws, data_)
  call auxputch("o3", "r", ows, 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", 8, data_)

  call sep_end_prog()
  call exit()
end program
