program main
  use sep
  use base_ioutil_mod
  use base_wavelt_mod
  implicit none

  type(axis)         :: asx, asy, agx, agy, at
  type(axis)         :: ax, ay, az, atmp
  integer            :: wvlt_type, it, igx, igx2, ix, isx, index_z_refl
  real               :: t0, t1, fm, t, wtx, PI, x
  real     , pointer :: tsou(:,:,:,:,:), vmod(:,:,:), refl(:,:,:), trec(:,:,:,:,:)
  character(len=256) :: trec_, tsou_, vmod_, refl_

  PI = 3.14159265359

  call initpar()

  trec_ = "trec"
  vmod_ = "vmod"
  tsou_ = "tsou"
  refl_ = "refl"

  asx%n = 3   ; asx%o = -600; asx%d = 200;
  asy%n = 1   ; asy%o = 0  ; asy%d = 10;
  agx%n = 401 ; agx%o = -2000  ; agx%d = 10;
  agy%n = 1   ; agy%o = 0  ; agy%d = 10;
  at%n  = 4000; at%o  = 0  ; at%d  = 0.002

  ax%n = 401; ax%o = -2000; ax%d = 10
  ay%n = 1  ; ay%o = 0; ay%d = 1
  az%n = 151; az%o = 0; az%d = 10

  !allocate(tsou(at%n, agx%n, agy%n, asx%n, asy%n))
  allocate(tsou(at%n,1,1,1,1))
  allocate(vmod(ax%n, ay%n, az%n))
  allocate(trec(at%n, agx%n, agy%n, asx%n, asy%n))
  allocate(refl(ax%n, ay%n, az%n))

  call houaxis(vmod_, ax, 1)
  call houaxis(vmod_, ay, 2)
  call houaxis(vmod_, az, 3)

  call houaxis(refl_, ax, 1)
  call houaxis(refl_, ay, 2)
  call houaxis(refl_, az, 3)

  call houaxis(trec_, at , 1)
  call houaxis(trec_, agx, 2)
  call houaxis(trec_, agy, 3)
  call houaxis(trec_, asx, 4)
  call houaxis(trec_, asy, 5)

atmp%n = 1; atmp%d = 1; atmp%o = 0.
  call houaxis(tsou_, at , 1)
  call houaxis(tsou_, atmp, 2)
  call houaxis(tsou_, atmp, 3)
  call houaxis(tsou_, asy, 4)
  call houaxis(tsou_, asy, 5)

  t0  = 1.1; t1=1.3
  igx = 1.5 + (600-agx%o)/agx%d
  !isx = 1.5 + (-600-agx%o)/agx%d

  igx2 = 1.5 + (900-agx%o)/agx%d

  wvlt_type = 1
  fm  = 20

  do isx=1, asx%n
     do it=1, at%n
        t = at%o + (it-1)*at%d
        tsou(it,1,1,1,1) = wavelet(t-0.1, fm, wvlt_type)

!     do ix=isx-2, isx+2
!        wtx = 0.5 + 0.5 * cos(PI * real(ix-isx)/3)
!        tsou(it, ix, 1, 1, 1) = wtx*wavelet(t-0.1, fm, wvlt_type)
!     end do

        do ix=igx-2, igx+2
           wtx = 0.5 + 0.5 * cos(PI * real(ix-igx)/3)
           trec(it, ix, 1, isx, 1) = wtx*wavelet(t-t0, fm, wvlt_type)
        end do

!     do ix=igx2-2, igx2+2
!        wtx = 0.5 + 0.5 * cos(PI * real(ix-igx2)/3)
!        trec(it, ix, 1, 1, 1) = trec(it,ix,1,1,1)+wtx*wavelet(t-t1, fm, wvlt_type)
!     end do

     end do
  end do

!  do ix=ax%n/2+1, ax%n
!     x = (ix-1-ax%n/2)*ax%d
!     vmod(ix,:,:) = 2000+(x)*0.77
!  end do
!  do ix=ax%n/2, 1, -1
!     x = (ax%n/2-(ix-1))*ax%d
!     vmod(ix,:,:) = 2000+(x)*0.77
!  end do

!  do ix=1, ax%n
!     x = ax%o + (ix-1)*ax%d
!     vmod(ix,:,:) = 800 + (x-ax%o)*0.77
!  end do

  vmod = 2000
!vmod = 2000
!vmod(1:ax%n/2,:,:) = 500
  refl = 0.
  index_z_refl = 1.5 + (800-az%o)/az%d
  refl(:,:,index_z_refl) = 1.

  call sep_write(tsou, tsou_)
  call sep_write(trec, trec_)
  call sep_write(vmod, vmod_)
  call sep_write(refl, refl_)

end program 
