program SR2MH_
  use sep
  use wf_util
  use wf_dump

  !! input:  receiver-source-time
  !! output: midpoint-offset-time

  !! Warning: ds==dr !!

  implicit none

  type(axis) :: at,ar,as,am,ah
  integer    :: it,ir,is,im,ih
  real       ::    rr,ss,mm,hh, dd

  real       :: rmin,rmax
  real       :: smin,smax
  real       :: mmin,mmax
  real       :: hmin,hmax

  real, pointer :: rsdat(:,:)
  real, pointer :: mhdat(:,:)

  !----------------------------------------------------------------
  call sep_init()

  call hinaxis( "in",ar,1); ar%l="r"
  call hinaxis( "in",as,2); as%l="s"
  call hinaxis( "in",at,3); at%l="t"

  call printaxis(as)
  call printaxis(ar)

  !! receiver sampling = source sampling
  if(ar%d /= as%d) call seperr("ar%d /= as%d")
  dd = as%d/2
  write(0,*) "dd=",dd

  if( mod(as%o,dd)/=0) call seperr("s0 not on the R grid")
  if( mod(ar%o,dd)/=0) call seperr("r0 not on the S grid")
  !----------------------------------------------------------------
  mmin=ar%o
  mmax=ar%o
  hmin=as%o
  hmax=as%o
  do is=1,as%n;    ss=as%o + (is-1)*as%d
     do ir=1,ar%n; rr=ar%o + (ir-1)*ar%d
        mm=(rr+ss)/2
        hh=(rr-ss)/2

        if(mm>mmax) mmax=mm
        if(mm<mmin) mmin=mm

        if(hh>hmax) hmax=hh
        if(hh<hmin) hmin=hh
     end do
  end do

  mmin=floor(mmin/dd)*dd 
  mmax=floor(mmax/dd)*dd +dd
  hmin=floor(hmin/dd)*dd -dd
  hmax=floor(hmax/dd)*dd

  write(0,*) "mmin=",mmin
  write(0,*) "mmax=",mmax
  write(0,*) "hmin=",hmin
  write(0,*) "hmax=",hmax
  !----------------------------------------------------------------
  am%d=dd
  am%o=mmin
  am%n=1+floor((mmax-mmin)/dd)
  am%l="m"

  ah%d=dd
  ah%o=hmin
  ah%n=1+floor((hmax-hmin)/dd)
  ah%l="h"
  !----------------------------------------------------------------
  call houaxis("out",am,1)
  call houaxis("out",ah,2)
  call houaxis("out",at,3)
  !----------------------------------------------------------------
  call sep_close()
  !----------------------------------------------------------------
  call printaxis(am)
  call printaxis(ah)
  !----------------------------------------------------------------
  allocate(rsdat(ar%n,as%n)); rsdat=0.
  allocate(mhdat(am%n,ah%n)); mhdat=0.

  do it=1,at%n
     write(0,*) it,at%n
     call sep_read (rsdat)
     mhdat = 0

     do is=1,as%n;    ss=as%o + (is-1)*as%d
        do ir=1,ar%n; rr=ar%o + (ir-1)*ar%d
           mm=(rr+ss)/2
           hh=(rr-ss)/2

           im=1+(mm-am%o)/dd
           ih=1+(hh-ah%o)/dd

           mhdat(im,ih) = rsdat(ir,is)
!           mhdat(im,ih) = mhdat(im,ih) + 1
        end do
     end do

     call sep_write(mhdat)
  end do

  call exit(0)
end program SR2MH_


