module wemva_histos_mod
  use  wemva_struct_mod
  use  base_ioutil_mod
  use  sep
  implicit none

  integer, external, private :: sreed_window, srite_window, getch

contains
!================================================================
! read parameters from histories
!----------------------------------------------------------------    
  subroutine wemva_histo_data(dgrid, vgrid, igrid, pgrid, ftag_)
    type(dat_grid)      :: dgrid
    type(vel_grid)      :: vgrid
    type(img_grid)      :: igrid
    type(pro_grid)      :: pgrid
    type(wemva_ftag)    :: ftag_
    type(axis)          :: agxcs, agycs, asxcs, asycs 
    real                :: image_xmin_default, image_xmax_default, image_dx_default
    real                :: image_ymin_default, image_ymax_default, image_dy_default
    real                :: image_zmin_default, image_zmax_default, image_dz_default
    real                :: ohx_default, ohy_default
    integer             :: nshg

    !!.. read the receiver coordinates, which are assumed to be regular
    call hinaxis(ftag_%crec_, dgrid%agx, 1)
    call hinaxis(ftag_%crec_, dgrid%agy, 2)
    call hinaxis(ftag_%crec_, dgrid%aw , 3)
    !!.. total receiver traces for one shot
    dgrid%ntrg = dgrid%agx%n * dgrid%agy%n
    !!.. if source header is presented
    if (dgrid%hdrsou) then
       call auxin(ftag_%hsou_)
       call from_aux(ftag_%hsou_, "n2", dgrid%nshg)   !!.. total shots
       call from_aux(ftag_%crec_, "n4", nshg)
       if (dgrid%nshg /= nshg) call seperr("ERROR: hsou-nshg /= crec-nshg")
    else
       call hinaxis(ftag_%crec_, dgrid%asx, 4)
       call hinaxis(ftag_%crec_, dgrid%asy, 5)
       dgrid%nshg = dgrid%asx%n * dgrid%asy%n   !!.. total shots
    end if
    !!.. if composite sources are used 
    if (dgrid%comsou) then
       call hinaxis(ftag_%csou_, agxcs, 1)
       call hinaxis(ftag_%csou_, agycs, 2)
       call hinaxis(ftag_%csou_, asxcs, 4)
       call hinaxis(ftag_%csou_, asycs, 5)
       dgrid%ntrs = agxcs%n * agycs%n           !!.. size of composite source
       dgrid%nshs = asxcs%n * asycs%n           !!.. number of composite sources
       if (dgrid%ntrs /= dgrid%ntrg) call seperr("ERROR: dgrid%ntrs /= dgrid%ntrg")
       if (dgrid%nshs /= dgrid%nshg) call seperr("ERROR: dgrid%nshs /= dgrid%nshg")
    else
       dgrid%ntrs = 1
       dgrid%nshs = 1
    end if

    if (dgrid%offset) then
       dgrid%gxmin_save = dgrid%agx%o
       dgrid%gymin_save = dgrid%agy%o
    end if

    !!.. velocity parameter initialization
    call hinaxis(ftag_%bslw_, vgrid%ax, 1)
    call hinaxis(ftag_%bslw_, vgrid%ay, 2)
    call hinaxis(ftag_%bslw_, vgrid%az, 3)

    !!.. the default image size is set to be the same as the velocity model
    image_xmin_default = vgrid%ax%o
    image_xmax_default = vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d
    image_dx_default   = vgrid%ax%d

    image_ymin_default = vgrid%ay%o
    image_ymax_default = vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d
    image_dy_default   = vgrid%ay%d

    image_zmin_default = vgrid%az%o
    image_zmax_default = vgrid%az%o + (vgrid%az%n-1)*vgrid%az%d
    image_dz_default   = vgrid%az%d

    call from_param("image_xmin", igrid%image_xmin, image_xmin_default)
    call from_param("image_xmax", igrid%image_xmax, image_xmax_default)
    call from_param("image_ymin", igrid%image_ymin, image_ymin_default)
    call from_param("image_ymax", igrid%image_ymax, image_ymax_default)
    call from_param("image_zmin", igrid%image_zmin, image_zmin_default)
    call from_param("image_zmax", igrid%image_zmax, image_zmax_default)
    call from_param("image_dx"  , igrid%image_dx  , image_dx_default  )
    call from_param("image_dy"  , igrid%image_dy  , image_dy_default  )
    call from_param("image_dz"  , igrid%image_dz  , image_dz_default  )

    igrid%image_xmin = vgrid%ax%o + int(1.5+(igrid%image_xmin-vgrid%ax%o)/vgrid%ax%d-1)*vgrid%ax%d
    igrid%image_ymin = vgrid%ay%o + int(1.5+(igrid%image_ymin-vgrid%ay%o)/vgrid%ay%d-1)*vgrid%ay%d
    igrid%image_zmin = vgrid%az%o + int(1.5+(igrid%image_zmin-vgrid%az%o)/vgrid%az%d-1)*vgrid%az%d

    igrid%image_xmax = vgrid%ax%o + int(1.5+(igrid%image_xmax-vgrid%ax%o)/vgrid%ax%d-1)*vgrid%ax%d
    igrid%image_ymax = vgrid%ay%o + int(1.5+(igrid%image_ymax-vgrid%ay%o)/vgrid%ay%d-1)*vgrid%ay%d
    igrid%image_zmax = vgrid%az%o + int(1.5+(igrid%image_zmax-vgrid%az%o)/vgrid%az%d-1)*vgrid%az%d

    if(igrid%image_xmin < image_xmin_default) call seperr("image_xmin is out of the input velocity model")
    if(igrid%image_xmax > image_xmax_default) call seperr("image_xmax is out of the input velocity model")
    if(igrid%image_ymin < image_ymin_default) call seperr("image_ymin is out of the input velocity model")
    if(igrid%image_ymax > image_ymax_default) call seperr("image_ymax is out of the input velocity model")
    if(igrid%image_zmin < image_zmin_default) call seperr("image_zmin is out of the input velocity model")
    if(igrid%image_zmax > image_zmax_default) call seperr("image_zmax is out of the input velocity model")

    call from_param("nxpad_beg", pgrid%nxpad_beg, 0)
    call from_param("nxpad_end", pgrid%nxpad_end, 0)
    call from_param("nypad_beg", pgrid%nypad_beg, 0)
    call from_param("nypad_end", pgrid%nypad_end, 0)
    call from_param("nzpad_beg", pgrid%nzpad_beg, 0)
    call from_param("nzpad_end", pgrid%nzpad_end, 0)

    call from_param("nhx", igrid%ahx%n, 1         )
    call from_param("nhy", igrid%ahy%n, 1         )
    call from_param("nhz", igrid%ahz%n, 1         )

    call from_param("dhx", igrid%ahx%d, vgrid%ax%d)
    call from_param("dhy", igrid%ahy%d, vgrid%ay%d)
    call from_param("dhz", igrid%ahz%d, vgrid%az%d)

    call from_param("ohx", igrid%ahx%o, 0.        ) 
    call from_param("ohy", igrid%ahy%o, 0.        )
    call from_param("ohz", igrid%ahz%o, 0.        )
 
    !!.. define the image grid
    igrid%ax%o = igrid%image_xmin
    igrid%ax%d = igrid%image_dx
    igrid%ax%n = 1.5 + (igrid%image_xmax - igrid%image_xmin)/igrid%image_dx

    igrid%ay%o = igrid%image_ymin
    igrid%ay%d = igrid%image_dy
    igrid%ay%n = 1.5 + (igrid%image_ymax - igrid%image_ymin)/igrid%image_dy

    igrid%az%o = igrid%image_zmin
    igrid%az%d = igrid%image_dz
    igrid%az%n = 1.5 + (igrid%image_zmax - igrid%image_zmin)/igrid%image_dz

    return
  end subroutine

!===========================================================================
! create image history
!  adj=F, input history file for reading data
!  adj=T, output history file for writing data
!---------------------------------------------------------------------------
  subroutine wemva_histo_imag(igrid, imag_)
    type(img_grid)   :: igrid
    character(len=*) :: imag_

    call houaxis(imag_, igrid%ahx, 1)
    call houaxis(imag_, igrid%ahy, 2)
    call houaxis(imag_, igrid%ax , 3)
    call houaxis(imag_, igrid%ay , 4)
    call houaxis(imag_, igrid%az , 5)

    call auxputch("label1", "s", "hx", imag_)
    call auxputch("label2", "s", "hy", imag_)
    call auxputch("label3", "s", "mx", imag_)
    call auxputch("label4", "s", "my", imag_)
    call auxputch("label5", "s", "mz", imag_)

    return
  end subroutine 
!=====================================================================================
  subroutine wemva_histo_dslw(vgrid, dslw_)
    type(vel_grid)   :: vgrid
    character(len=*) :: dslw_

    call houaxis(dslw_, vgrid%ax , 1)
    call houaxis(dslw_, vgrid%ay , 2)
    call houaxis(dslw_, vgrid%az , 3)

    call auxputch("label1", "s", "x", dslw_)
    call auxputch("label2", "s", "y", dslw_)
    call auxputch("label3", "s", "z", dslw_)   
    return
  end subroutine
!=====================================================================================
  subroutine wemva_histo_wfld(pgrid, aws, wfld_)
    type(pro_grid)   :: pgrid
    type(axis)       :: aws
    character(len=*) :: wfld_

    call houaxis(wfld_, pgrid%ax, 1)
    call houaxis(wfld_, pgrid%ay, 2)
    call houaxis(wfld_, aws     , 3)
    call houaxis(wfld_, pgrid%az, 4)
    call auxputch("esize", "i", 8, wfld_)

    call auxputch("label1", "s", "x"    , wfld_)
    call auxputch("label2", "s", "y"    , wfld_)
    call auxputch("label3", "s", "freq.", wfld_)
    call auxputch("label4", "s", "z"    , wfld_)
    return
  end subroutine

end module



