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
  real   , dimension(:)    , pointer :: phase_random
  complex, dimension(:)    , pointer :: soufnc
  complex, dimension(:,:,:), pointer :: data
  integer                            :: isx, isy, iws, igx, igy
  real                               :: gx, gy, sx, sy, ws, refx, refy
  real                               :: p, phase_delay, angle, v0, phase
  character(len=256)                 :: soufnc_, data_
  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"

  call auxin(soufnc_)
  call auxout(data_)

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

  allocate(phase_random(nws))
  allocate(soufnc(nws))
  allocate(data(ngx,ngy,nws))

  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

  call random_seed()

  !! now only 2-D
  do isy=1, nsy
  do isx=1, nsx

     do iws=1, nws
        ws = ows + (iws-1)*dws
        do igy=1, ngy
        do igx=1, ngx
           gx = ogx + (igx-1)*dgx
           call random_number(phase)
           phase = phase * 2*PI 
           data(igx,igy,iws) = soufnc(iws)*cmplx(cos(phase), sin(phase))
        end do
        end do
     end do

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

  end do
  end do

  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
