program main
  use sep
  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 auxin(datain_)

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

  call from_history("n1", ngx, 1); call from_history("d1", dgx, 1.); call from_history("o1", ogx, 0.)
  call from_history("n2", ngy, 1); call from_history("d2", dgy, 1.); call from_history("o2", ogy, 0.) 
  call from_history("n3", nws, 1); call from_history("d3", dws, 1.); call from_history("o3", ows, 0.)
  call from_history("n4", nsx, 1); call from_history("d4", dsx, 1.); call from_history("o4", osx, 0.)
  call from_history("n5", nsy, 1); call from_history("d5", dsy, 1.); call from_history("o5", 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
  !! now only 2-D
  do isy=1, nsy
  do isx=1, nsx

     do iws=1, nws
        do igy=1, ngy
        do igx=1, ngx
           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(dataou_, data, ngx, ngy, nws, nsx, nsy, isx, isy, "rite")

  end do
  end do

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

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

  call auxputch("n3", "i", nws, dataou_)
  call auxputch("d3", "r", dws, dataou_)
  call auxputch("o3", "r", ows, dataou_)

  call auxputch("n4", "i", nsx, dataou_)
  call auxputch("d4", "r", dsx, dataou_)
  call auxputch("o4", "r", osx, dataou_)

  call auxputch("n5", "i", nsy, dataou_)
  call auxputch("d5", "r", dsy, dataou_)
  call auxputch("o5", "r", osy, dataou_)

  call auxputch("esize", "i", 8, dataou_)

  call sep_end_prog()
  call exit()
end program
