module grid_mod

  use type_mod

  type(cgrid), public :: grid
  !complex, dimension(:,:,:), pointer, private :: wfld_d, wfld_u
  !complex, dimension(:,:,:), pointer, private :: cdata , csfnc

contains

  !!.. calculate parameters for the computation grid 
  subroutine get_computation_grid(asx_in, asy_in, agx_in, agy_in, az_in)
    type(axis) :: asx_in, asy_in
    type(axis) :: agx_in, agy_in
    type(axis) :: az_in
    real       :: grid_xmin, grid_ymin
    real       :: grid_xmax, grid_ymax

    write(0,*)"asx_in", asx_in
    write(0,*)"asy_in", asy_in
    write(0,*)"agx_in", agx_in
    write(0,*)"agy_in", agy_in


    !!.. the spacing of the computation grid is assumed to 
    !!.. be the same as the receiver spacing
    grid%x%d = agx_in%d
    grid%y%d = agy_in%d

    !!.. the source location for a shot gather
    grid%sx = asx_in%o
    grid%sy = asy_in%o

    !!.. the minimum and maximum coordinates for the recorded receiver wavefield
    grid%gxmin = agx_in%o
    grid%gymin = agy_in%o
    grid%gxmax = agx_in%o + (agx_in%n-1)*agx_in%d
    grid%gymax = agy_in%o + (agy_in%n-1)*agy_in%d
  
    !!.. the coordinate for the depth direction
    grid%z  = az_in

    !!.. compute the migration aperture
    !!.. for simplicity, currently I am just using the maximum distance
    !!.. between the source and the receiver (better make it adaptive)

    grid_xmin = grid%gxmin; grid_xmax = grid%gxmax
    grid_ymin = grid%gymin; grid_ymax = grid%gymax

    if (grid%sx .ge. grid%gxmax) grid_xmax = grid%sx
    if (grid%sy .ge. grid%gymax) grid_ymax = grid%sy

    if (grid%sx .le. grid%gxmin) grid_xmin = grid%sx
    if (grid%sy .le. grid%gymin) grid_ymin = grid%sy

    write(0,*)"grid_xmin", grid_xmin
    write(0,*)"grid_xmax", grid_xmax
    write(0,*)"grid_ymin", grid_ymin
    write(0,*)"grid_ymax", grid_ymax

    !!.. no padding now, I'll do it later!!
    grid%x%o = grid_xmin
    grid%y%o = grid_ymin
    grid%x%n = 1.5 + (grid_xmax - grid_xmin)/grid%x%d
    grid%y%n = 1.5 + (grid_ymax - grid_ymin)/grid%y%d
 
    !!.. determine the source location on the computation grid
    grid%iosx = 1.5 + (grid%sx - grid%x%o) / grid%x%d
    grid%iosy = 1.5 + (grid%sy - grid%y%o) / grid%y%d

    !!.. determine the first receiver location on the computation grid
    grid%iogx = 1.5 + (agx_in%o - grid%x%o) / grid%x%d
    grid%iogy = 1.5 + (agy_in%o - grid%y%o) / grid%y%d

  
    write(0,*)"grid%x", grid%x%n, grid%x%o, grid%x%d
    write(0,*)"grid%y", grid%y%n, grid%y%o, grid%y%d
    write(0,*)"grid%iosx,iosy", grid%iosx, grid%iosy
    write(0,*)"grid%iogx,iogy", grid%iogx, grid%iogy

  end subroutine


  !!.. initialize down and up wavefields according to the computation grid
  !subroutine wfld_init(wfld_d_in, wfld_u_in, cdata_in, csfnc_in)
  !  complex, dimension(:,:,:), target :: wfld_d_in, wfld_u_in  !!.. wfld_d(sx, sy, w), wfld_u(gx, gy, w) , after gridding
  !  complex, dimension(:,:,:), target :: cdata_in , csfnc_in   !!.. csfnc(sx, sy, w) , cdata(gx, gy, w), before gridding

  !  wfld_d => wfld_d_in; wfld_u => wfld_u_in
  !  cdata  => cdata_in ; csfnc  => csfnc_in

  !  wfld_d = 0. ; wfld_u = 0.
    
  !  wfld_d(grid%iosx, grid%iosy, :) = csfnc(1, 1, :)

  !  wfld_u(grid%iogx:grid%iogx+size(wfld_u, 1)-1, grid%iogy:grid%iogy+size(wfld_u, 2)-1, :) = cdata

  !  nullify(wfld_d, wfld_u)
  !  nullify(cdata, csfnc)

  !end subroutine

end module













