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 :: datain, dataou

  integer                            :: igx, igy, iws, isx, isy
  real                               :: gx, phase_delay, ws, p, angle, v0
  logical                            :: offset
  character(len=256)                 :: datain_, dataou_

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


  datain_ = "input"
  dataou_ = "output"

  call auxin(datain_)
  call auxout(dataou_)

  call from_aux(datain_, "n1", ngx); call from_aux(datain_, "d1", dgx); call from_aux(datain_, "o1", ogx)
  call from_aux(datain_, "n2", ngy); call from_aux(datain_, "d2", dgy); call from_aux(datain_, "o2", ogy)
  call from_aux(datain_, "n3", nws); call from_aux(datain_, "d3", dws); call from_aux(datain_, "o3", ows)
  call from_aux(datain_, "n4", nsx); call from_aux(datain_, "d4", dsx); call from_aux(datain_, "o4", osx)
  call from_aux(datain_, "n5", nsy); call from_aux(datain_, "d5", dsy); call from_aux(datain_, "o5", osy)

  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 from_param("v0", v0, 1500.)
  call from_param("angle", angle, 0.)
  call from_param("offset", offset, .false.)

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

  allocate(datain(ngx, ngy, nws))
  allocate(dataou(ngx, ngy, nws))

write(0,*)"ngx, ngy, nws, nsx, nsy", ngx, ngy, nws, nsx, nsy


  !.. only works for 2D now

  do isy=1, nsy
  do isx=1, nsx

     call iodata(datain_, datain, ngx, ngy, nws, nsx, nsy, isx, isy, "reed")

     do iws=1, nws
        ws = ows + (iws-1)*dws

        do igy=1, ngy
        do igx=1, ngx
           gx = ogx + (igx-1)*dgx
           phase_delay = ws*p*(gx-ogx)
           dataou(igx,igy,iws) = datain(igx,igy,iws) * cmplx(cos(phase_delay), sin(phase_delay))

        end do
        end do
     end do

     call iodata(dataou_, dataou, ngx, ngy, nws, nsx, nsy, isx, isy, "rite")

  end do
  end do

  call sep_end_prog()
  call exit()
end program







