module planemod2
use sep
use adj_mod
implicit none
  integer, private :: nt,nx,nq,ns
  real, private :: t0,q0,s0
  real, private :: dt,dq,ds

!  integer, private :: nthreads
!  integer, private :: verb
  real,allocatable,save,public:: head(:,:)
  real,save,public:: sx,sy,gz,pi,ox,oy,os,x,y,xs2,s2
contains
!----------------------------------------------------------------
  subroutine planemod_init(t0_in,dt_in,q0_in,dq_in,nt_in,nx_in,nq_in )

    real    :: t0_in,dt_in, q0_in,dq_in
    integer :: nt_in,nx_in,nq_in
    integer :: trc
    real:: lx,ly,s

    t0 = t0_in
    dt = dt_in
    q0 = q0_in
    dq = dq_in
    nt = nt_in
    nx = nx_in
    nq = nq_in

    allocate(head(2,nx));    head=0.
    call sep_read(head,"geom")
    call from_param("shottrc",trc,1)
    sx=head(2,trc)
    sy=head(1,trc)
    write(0,*)'  Location of correlated shotpoint (sx,sy):',sx,sy
    ox=2000. +sx !parmin(1)+sx 
    oy=-40000.  +sy !parmin(2)+sy
    !! distant source location
    x=ox
    y=oy
    !! distance from source to virtual shot location
    lx=(x-sx)**2
    ly=(y-sy)**2
    xs2=lx+ly    
    s=1./1456.
    s2=s*s
  end subroutine planemod_init
!----------------------------------------------------------------
!!page 144 "lessons in seismic computing",Slotnick
!  integer function planemod_op(modl,j)
  integer function planemod_op(adja,adda,modll,dat) result(stat)
    logical,intent(in) :: adja,adda
    real,dimension(:)  :: modll,dat 
    call adjnull (adja,adda,modll,dat )
    call planemod_op2(adja,modll,dat )
    stat=0
  end function
!----------------------------------------------------------------
  subroutine planemod_op2(adj,modl,data)
    logical,intent(in) :: adj
    real,dimension(nt,nq)::modl
    real,dimension(nt,nx)::data

    real:: s,lx,ly,lt,x2,z,t,tshot!,s2
    integer:: i,j,k,it,itt

    if(adj)then
       slow:do j=1,nq
          !! vel squared
          s=q0+(j-1)*dq
!          s2=s*s
          !! event at source trace (correlated gathers) always at t=0
          tshot=sqrt(s2*xs2)  !xs2 calcd in _init
          trc:do i=1,nx
             lx=(head(2,i)-x)**2
             ly=(head(1,i)-y)**2
             x2=lx+ly
             lt=(sqrt(s2*x2)-tshot)*s
             itt=1.5+(lt-t0)/dt
             time:do k=1,nt
                it=itt+k-1
                if ( (it >= nt).or.(it<1) ) then
                   cycle                           ! out of the range
                end if
                modl(k,j)= modl(k,j) + data(it,i)
             end do time
       end do trc
    end do slow
 else ! forward
    trcf:do i=1,nx
       lx=(head(2,i)-x)**2
       ly=(head(1,i)-y)**2
       x2=lx+ly
       slowf:do j=1,nq
          s=q0+(j-1)*dq
!          s2=s*s
          tshot=sqrt(s2*xs2)
          lt=(sqrt(s2*x2)-tshot)*s
          itt=1.5+ (lt-t0)/dt
          timef:do k=1, nt
             it=itt+k-1
             if ( (it >= nt) .or. (it<1) ) then
                cycle                          ! out of the range
             end if
             data(it,i)=data(it,i)+ modl(k,j)
          end do timef
       end do slowf
    end do trcf
 end if

end subroutine planemod_op2
!----------------------------------------------------------------
  subroutine velsimp6_init ( t0_in,dt_in,s0_in,ds_in,nt_in&
    &,nx_in,ns_in )

!$  integer, external  :: omp_get_num_threads,omp_get_thread_num

    real    :: t0_in 
    real    :: dt_in 
    real    :: s0_in 
    real    :: ds_in 

    integer    :: nt_in 
    integer    :: nx_in 
    integer    :: ns_in 
    character(len=128) :: desc
!write(0,*)' vel init', sx,sy,head(2,1),head(1,1)
    t0 = t0_in 
    dt = dt_in 
    s0 = s0_in 
    ds = ds_in 
    nt = nt_in 
    nx = nx_in 
    ns = ns_in

!    nthreads = 1
!!$OMP PARALLEL
!!$      nthreads = omp_get_num_threads()
!!$OMP END PARALLEL
!    write(0,*) 'nthreads = ',nthreads

  end subroutine velsimp6_init

!----------------------------------------------------------------
  integer function velsimp6_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 velsimp6_lop2(adj,modl,data )
write(0,*)' vel-op'
    stat=0

  end function velsimp6_lop
!----------------------------------------------------------------  
  subroutine velsimp6_lop2(adj,modl,data)

    logical,intent(in) :: adj 
    real, dimension (nt,ns)  :: modl 
    real, dimension (nt,nx)  :: data 
    integer it,iz,ix,is,nread,nth,ith
    real dx,s,svx,t,z,wt,lx,ly

if (adj) then

!!$OMP PARALLEL DO PRIVATE(is,s,ix,x,iz,z,svx,t,it,wt)
     do is= 1, ns

      s = s0 + (is-1) * ds               ! slowness

      do ix= 1, nx  
!        x = x0 + (ix-1) * dx             ! offset
         lx=(head(2,ix)-sx)**2
         ly=(head(1,ix)-sy)**2
         dx=sqrt(lx+ly)

        do iz= 2, nt  

          z = t0 + (iz-1) * dt
          svx = abs( s * dx)               ! normalized offset
          t = sqrt( z * z + svx * svx)
          it = 1.5 + (t-t0)/dt 

           if (it>nt) then
            exit
           end if

          wt = (z/t) / sqrt(t)

            modl(iz,is) =   modl(iz,is) + data(it,ix)*wt

        end do
      end do

    end do 
!!$OMP END PARALLEL DO 

else

!!$OMP PARALLEL DO PRIVATE(ix,x,is,s,iz,z,svx,t,it,wt)
   do ix= 1, nx  

!     x = x0 + (ix-1) * dx                ! offset
      lx=(head(2,ix)-sx)**2
      ly=(head(1,ix)-sy)**2
      dx=sqrt(lx+ly)
      do is= 1, ns

      s = s0 + (is-1) * ds               ! slowness

        do iz= 2, nt  

          z = t0 + (iz-1) * dt
          svx = abs( s * dx)               ! normalized offset
          t = sqrt( z * z + svx * svx)
          it = 1.5 + (t-t0)/dt 

           if (it>nt) then
            exit
           end if

          wt = (z/t) / sqrt(t)

            data(it,ix) =   data(it,ix) + modl(iz,is)*wt

        end do
      end do

    end do 
!!$OMP END PARALLEL DO 

end if

end subroutine velsimp6_lop2

end module planemod2
