module wem_histos_complex_mod
  use  sep
  use  base_struct_mod
  use  base_ioutil_mod
  use  wem_struct_mod
  implicit none

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

contains
!================================================================
! read parameters from histories
!----------------------------------------------------------------    
  subroutine wem_histo_data(adj, dgrid, vgrid, igrid, pgrid, ftag_)
    logical             :: adj
    type(dat_grid)      :: dgrid
    type(vel_grid)      :: vgrid
    type(img_grid)      :: igrid
    type(pro_grid)      :: pgrid
    type(wem_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

    call auxinout(ftag_%crec_)
    call auxinout(ftag_%csou_)
    call auxinout(ftag_%vmod_)

!   padding for the wavefield propagation
!   -------------------------------------
    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)

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

!   migration
!   ---------
    if (adj) then

!      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_)
!         total shots
!         -----------
          call from_aux(ftag_%hsou_, "n2", dgrid%nshg)
          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)
!         total shots
!         -----------
          dgrid%nshg = dgrid%asx%n * dgrid%asy%n
       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

    else

       call from_param("ngx", dgrid%agx%n, 1 )
       call from_param("dgx", dgrid%agx%d, 1.)
       call from_param("ogx", dgrid%agx%o, 0.)

       call from_param("ngy", dgrid%agy%n, 1 )
       call from_param("dgy", dgrid%agy%d, 1.)
       call from_param("ogy", dgrid%agy%o, 0.)

       call hinaxis(ftag_%csou_, dgrid%aw, 3 )

       call from_param("nsx", dgrid%asx%n, 1 )
       call from_param("dsx", dgrid%asx%d, 1.)
       call from_param("osx", dgrid%asx%o, 0.)

       call from_param("nsy", dgrid%asy%n, 1 )
       call from_param("dsy", dgrid%asy%d, 1.) 
       call from_param("osy", dgrid%asy%o, 0.)

       dgrid%ntrg = dgrid%agx%n*dgrid%agy%n
       dgrid%nshg = dgrid%asx%n*dgrid%asy%n

       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

       call houaxis(ftag_%crec_, dgrid%agx, 1)
       call houaxis(ftag_%crec_, dgrid%agy, 2)
       call houaxis(ftag_%crec_, dgrid%aw , 3)
       call houaxis(ftag_%crec_, dgrid%asx, 4)
       call houaxis(ftag_%crec_, dgrid%asy, 5)
       call auxputch("esize", "i", 8, ftag_%crec_)

       if (dgrid%offset) then
          call auxputch("label1", "s", "offx", ftag_%crec_)
          call auxputch("label2", "s", "offy", ftag_%crec_)
       else
          call auxputch("label1", "s", "gx"  , ftag_%crec_)
          call auxputch("label2", "s", "gy"  , ftag_%crec_)
       end if

       call auxputch("label3", "s", "freq", ftag_%crec_)
       call auxputch("label4", "s", "sx"  , ftag_%crec_)
       call auxputch("label5", "s", "sy"  , ftag_%crec_)

    end if
    return
  end subroutine

!===========================================================================
! create image history
!  adj=F, input history file for reading data
!  adj=T, output history file for writing data
!---------------------------------------------------------------------------
  subroutine wem_histo_imag(adj, igrid, vgrid, ftag_)
    logical        :: adj
    type(img_grid) :: igrid
    type(vel_grid) :: vgrid
    type(wem_ftag) :: ftag_

    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

!   migration, write the image coordinates
!   -------------------------------------- 
    if (adj) then
!      the default image size is set to be the same as the velocity model
!      ------------------------------------------------------------------
       image_xmin_default=vgrid%ax%o; image_dx_default=vgrid%ax%d; image_xmax_default=vgrid%ax%o+(vgrid%ax%n-1)*vgrid%ax%d
       image_ymin_default=vgrid%ay%o; image_dy_default=vgrid%ay%d; image_ymax_default=vgrid%ay%o+(vgrid%ay%n-1)*vgrid%ay%d
       image_zmin_default=vgrid%az%o; image_dz_default=vgrid%az%d; image_zmax_default=vgrid%az%o+(vgrid%az%n-1)*vgrid%az%d

!      read the imaging parameters
!      ---------------------------
       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  )

!      safe guarding
!      -------------
       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")

!      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

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

!      write out image history
!      -----------------------
       call houaxis(ftag_%imag_, igrid%ahx, 1)
       call houaxis(ftag_%imag_, igrid%ahy, 2)
       call houaxis(ftag_%imag_, igrid%ax , 3)
       call houaxis(ftag_%imag_, igrid%ay , 4)
       call houaxis(ftag_%imag_, igrid%az , 5)

       call auxputch("esize" , "i", 8   , ftag_%imag_)
       call auxputch("label1", "s", "hx", ftag_%imag_)
       call auxputch("label2", "s", "hy", ftag_%imag_)
       call auxputch("label3", "s", "mx", ftag_%imag_)
       call auxputch("label4", "s", "my", ftag_%imag_)
       call auxputch("label5", "s", "mz", ftag_%imag_)
    else
       call hinaxis(ftag_%imag_, igrid%ahx, 1)
       call hinaxis(ftag_%imag_, igrid%ahy, 2)
       call hinaxis(ftag_%imag_, igrid%ax , 3)
       call hinaxis(ftag_%imag_, igrid%ay , 4)
       call hinaxis(ftag_%imag_, igrid%az , 5)
    end if
    return
  end subroutine 
!--------------------------------------------------------------------------------------
  subroutine wem_histo_hess(hgrid, igrid, hess_)
    type(hes_grid)   :: hgrid
    type(img_grid)   :: igrid
    character(len=*) :: hess_
    hgrid%ax_modl = igrid%ax
    hgrid%ay_modl = igrid%ay
    hgrid%az_modl = igrid%az

    call from_param("nx_offd", hgrid%ax_offd%n, 1 )
    call from_param("ox_offd", hgrid%ax_offd%o, 0.)
    call from_param("dx_offd", hgrid%ax_offd%d, 1.)

    call from_param("ny_offd", hgrid%ay_offd%n, 1 )
    call from_param("oy_offd", hgrid%ay_offd%o, 0.)
    call from_param("dy_offd", hgrid%ay_offd%d, 1.)

    call from_param("nz_offd", hgrid%az_offd%n, 1 )
    call from_param("oz_offd", hgrid%az_offd%o, 0.)
    call from_param("dz_offd", hgrid%az_offd%d, 1.)
 
    call auxinout(hess_)

    call houaxis(hess_, hgrid%ax_modl, 1)
    call houaxis(hess_, hgrid%ay_modl, 2)
    call houaxis(hess_, hgrid%az_modl, 3)
    call houaxis(hess_, hgrid%ax_offd, 4) 
    call houaxis(hess_, hgrid%ay_offd, 5)
    call houaxis(hess_, hgrid%az_offd, 6)
    call auxputch("label1", "s", "x_modl", hess_)
    call auxputch("label2", "s", "y_modl", hess_)
    call auxputch("label3", "s", "z_modl", hess_)
    call auxputch("label4", "s", "x_offd", hess_)
    call auxputch("label5", "s", "y_offd", hess_)
    call auxputch("label6", "s", "z_offd", hess_)

    return
  end subroutine

end module
