!!-----------------------------------------------------------
!! Author: Yaxun Tang, tang@sep.stanford.edu
!!-----------------------------------------------------------
module base_progrd_mod
  use  base_struct_mod
  implicit none

contains
  !!.. calculate the propagation grid
  !!.. now only for shot-profile migration
  !!.. composite shot migration will be included later
  subroutine base_progrd_gen3d(pgrid, dgrid, vgrid)
    type(pro_grid) :: pgrid
    type(dat_grid) :: dgrid
    type(vel_grid) :: vgrid
    real           :: rx_min, ry_min, rz_min
    real           :: rx_max, ry_max, rz_max
    real           :: pgrid_xmin, pgrid_xmax
    real           :: pgrid_ymin, pgrid_ymax
    real           :: pgrid_zmin, pgrid_zmax

    !!.. initial propagation grid set to be the same
    !!.. as the entire velocity model
    pgrid%ax = vgrid%ax
    pgrid%ay = vgrid%ay
    pgrid%az = vgrid%az

    pgrid%sx = dgrid%sx
    pgrid%sy = dgrid%sy
    pgrid%sz = dgrid%sz

    !!.. calculate the min and max of the receiver locations
    !!.. for a single shot gather (or composite shot gather)
    rx_min = dgrid%agx%o
    ry_min = dgrid%agy%o
    rz_min = dgrid%agz%o
    rx_max = dgrid%agx%o + (dgrid%agx%n-1)*dgrid%agx%d
    ry_max = dgrid%agy%o + (dgrid%agy%n-1)*dgrid%agy%d
    rz_max = dgrid%agz%o + (dgrid%agz%n-1)*dgrid%agz%d

    pgrid%gxmin = rx_min; pgrid%gxmax = rx_max
    pgrid%gymin = ry_min; pgrid%gymax = ry_max
    pgrid%gzmin = rz_min; pgrid%gzmax = rz_max

    !!.. determine the size of propagation grid
    pgrid_xmin = rx_min
    pgrid_xmax = rx_max
    pgrid_ymin = ry_min
    pgrid_ymax = ry_max
    pgrid_zmin = pgrid%az%o
    pgrid_zmax = pgrid%az%o + (pgrid%az%n-1)*pgrid%az%d


    if (.not.dgrid%comsou) then
       if (pgrid%sx <= pgrid_xmin) pgrid_xmin = pgrid%sx
       if (pgrid%sy <= pgrid_ymin) pgrid_ymin = pgrid%sy
   
       if (pgrid%sx >= pgrid_xmax) pgrid_xmax = pgrid%sx
       if (pgrid%sy >= pgrid_ymax) pgrid_ymax = pgrid%sy 
    end if

    !!.. add padding
    pgrid_xmin = pgrid_xmin - pgrid%nxpad_beg*pgrid%ax%d
    pgrid_xmax = pgrid_xmax + pgrid%nxpad_end*pgrid%ax%d

    pgrid_ymin = pgrid_ymin - pgrid%nypad_beg*pgrid%ay%d
    pgrid_ymax = pgrid_ymax + pgrid%nypad_end*pgrid%ay%d

    pgrid_zmin = pgrid_zmin - pgrid%nzpad_beg*pgrid%az%d
    pgrid_zmax = pgrid_zmax + pgrid%nzpad_end*pgrid%az%d


    !!.. make sure the propagation grid does not go beyond
    !!.. the enitire velocity grid
!    pgrid_xmin = max(pgrid_xmin, vgrid%ax%o)
!    pgrid_xmax = min(pgrid_xmax, vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d)
!
!    pgrid_ymin = max(pgrid_ymin, vgrid%ay%o)
!    pgrid_ymax = min(pgrid_ymax, vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d)
!
!    pgrid_zmin = max(pgrid_zmin, vgrid%az%o)
!    pgrid_zmax = min(pgrid_zmax, vgrid%az%o + (vgrid%az%n-1)*vgrid%az%d)

    if (pgrid_xmin<vgrid%ax%o .or. pgrid_xmax>vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d) then
       write(0,*)"ERROR: pgrid_xmin < vgrid_xmin .or. pgrid_xmax > vgrid_xmax"
       write(0,*)"pgrid_xmin=", pgrid_xmin, " < ", "vgrid_xmin=", vgrid%ax%o
       write(0,*)"pgrid_xmax=", pgrid_xmax, " > ", "vgrid_xmax=", (vgrid%ax%n-1)*vgrid%ax%d
       stop
    end if

    if (pgrid_ymin<vgrid%ay%o .or. pgrid_ymax>vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d) then
       write(0,*)"ERROR: pgrid_ymin < vgrid_ymin .or. pgrid_ymax > vgrid_ymax"
       write(0,*)"pgrid_ymin=", pgrid_ymin, " < ", "vgrid_ymin=", vgrid%ay%o
       write(0,*)"pgrid_ymax=", pgrid_ymax, " > ", "vgrid_ymax=", (vgrid%ay%n-1)*vgrid%ay%d
       stop
    end if

    if (pgrid_zmin<vgrid%az%o .or. pgrid_zmax>vgrid%az%o + (vgrid%az%n-1)*vgrid%az%d) then
       write(0,*)"ERROR: pgrid_zmin < vgrid_zmin .or. pgrid_zmax > vgrid_zmax"
       write(0,*)"pgrid_zmin=", pgrid_zmin, " < ", "vgrid_zmin=", vgrid%az%o
       write(0,*)"pgrid_zmax=", pgrid_zmax, " > ", "vgrid_zmax=", (vgrid%az%n-1)*vgrid%az%d
       stop
    end if

    pgrid%ax%o = pgrid_xmin
    pgrid%ay%o = pgrid_ymin
    pgrid%az%o = pgrid_zmin
    pgrid%ax%n = 1.5 + (pgrid_xmax - pgrid_xmin)/pgrid%ax%d
    pgrid%ay%n = 1.5 + (pgrid_ymax - pgrid_ymin)/pgrid%ay%d
    pgrid%az%n = 1.5 + (pgrid_zmax - pgrid_zmin)/pgrid%az%d

    if (.not.dgrid%comsou) then
       !!.. calculate the index of source locations on the propagation grid
       pgrid%index_datpro_sx = 1.5 + (pgrid%sx-pgrid%ax%o)/pgrid%ax%d
       pgrid%index_datpro_sy = 1.5 + (pgrid%sy-pgrid%ay%o)/pgrid%ay%d
       pgrid%index_datpro_sz = 1.5 + (pgrid%sz-pgrid%az%o)/pgrid%az%d
    end if
 
    !!.. calculate the index of receiver locations on the propagation grid
    pgrid%index_recpro_xmin = 1.5 + (rx_min-pgrid%ax%o)/pgrid%ax%d
    pgrid%index_recpro_xmax = 1.5 + (rx_max-pgrid%ax%o)/pgrid%ax%d
    pgrid%index_recpro_ymin = 1.5 + (ry_min-pgrid%ay%o)/pgrid%ay%d
    pgrid%index_recpro_ymax = 1.5 + (ry_max-pgrid%ay%o)/pgrid%ay%d
    pgrid%index_recpro_zmin = 1.5 + (rz_min-pgrid%ay%o)/pgrid%az%d
    pgrid%index_recpro_zmax = 1.5 + (rz_max-pgrid%ay%o)/pgrid%az%d

    return
  end subroutine

end module


