module velsimpx  
! simple velocity transform
  use adj_mod

  implicit none
  integer, private :: nt,nx,ns,it,iz,ix,is 
  real, private :: t0,dt,x0,dx,s0,ds,x,s,sx,t,z,wt 
  logical, private :: lint 
  contains

  subroutine velsimpx_init ( t0_in,dt_in,x0_in,dx_in,s0_in,ds_in,nt_in&
    &,nx_in,ns_in,lint_in )
    real       :: t0_in,dt_in,x0_in,dx_in,s0_in,ds_in
    integer    :: nt_in,nx_in,ns_in 
    logical    :: lint_in 
    t0 = t0_in;dt = dt_in;x0 = x0_in;dx = dx_in;s0 = s0_in;ds = ds_in;
    nt = nt_in;nx = nx_in;ns = ns_in 
    lint = lint_in
  end subroutine

  function velsimpx_lop ( adj, add, modl, data) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: modl,data 
    call adjnull (adj,add,modl,data )
    call velsimpx_lop2(adj,add,modl,data )
    stat=0
  end function 

  subroutine velsimpx_lop2(adj,add,modl,data)
    logical,intent(in) :: adj,add 
    real, dimension (nt,ns)  :: modl 
    real, dimension (nt,nx)  :: data 
    real		     :: f,fx,gx
    integer		     :: i
    do is = 1,ns  
! slowness loop
      s = s0 + (is-1) * ds  
      do ix = 1,nx  
! offset loop
        x = x0 + (ix-1) * dx 
        do iz = 2,nt  
! zero-offset traveltime loop
          z = t0 + (iz-1) * dt 
          sx = abs( s * x)                     ! normalized offset
          t  = sqrt( z * z + sx * sx)          ! NMO
          it = 1.5 + (t-t0)/dt               
          if ( lint) then
! Do linear interpolation...
 		f=(t-t0)/dt
		i=f
		it=i+1
		if((it.ge.1).and.(it.lt.nt)) then
			fx=f-i
			gx=1.-fx
		end if
          end if 
          if ( it+1 .gt. nt) then
            exit                  ! out of the range
          end if
          wt = 1.0                             
! (t,x)-dependent weighting
          if ( adj) then
            if ( lint) then
             if (iz+1 < nt) then
              modl(iz,is) =           modl(iz,is) + gx*data(it,ix)*wt +fx*data(it+1,ix)*wt
             end if
            else
              modl(iz,is) =            modl(iz,is) + data(it,ix)*wt
            end if
          else
            if ( lint) then
             if (iz+1 < nt) then
              data(it,ix) =            data(it,ix) + gx*modl(iz,is)*wt
              data(it+1,ix) =            data(it+1,ix) + fx*modl(iz,is)*wt
             end if
            else
              data(it,ix) =            data(it,ix) + modl(iz,is)*wt
            end if
          end if
        end do 
! End time loop
      end do 
! End offset loop
    end do 
! End zero-offset traveltime loop
  end subroutine 
  subroutine velsimpx_close()
  end subroutine
end module
  
