module OW_slant
! 3D-linear Radon transform

  use OW_type
  use OW_parms_slt

  implicit none

  real,private,allocatable                   :: vv(:,:)

  contains

!===============================================================================
  function owslant_init(v_in) result(stat)
  real          :: v_in(img%z%n,img%xm%n)
  integer       :: i,j,stat
  allocate(vv(img%z%n,img%xm%n))
  do i=1,img%xm%n
     do j=1,img%z%n
        vv(j,i)=sum(v_in(1:j,i))/j
     end do
  end do
  stat=0
  end function
!===============================================================================
  function owslant(adj, add, modl, data) result(stat)
    integer                     :: stat 
    logical,intent(in)          :: adj,add 
    real                        :: modl(img%z%n,img%sh%n,img%xm%n) 
    real                        :: data(img%z%n,img%th%n,img%xm%n)

    if (adj) then
       modl=0.
    else
       data=0.
    end if
    call owslant_exe(adj,add,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine owslant_exe(adj,add,modl,data)
    logical,intent(in)        :: adj,add
    real		  	:: modl(:,:,:) 
    real		  	:: data(:,:,:)

!! internal variables
    integer		:: inode
    integer 		:: iax,ix,ixm,iah,ih,izm,izd,it
    integer		:: minx,maxx,indd,indm,idx,nodes,nah
    real			:: ax,tgx,xx,ah,tgh,hx,zm,zd,rdx,fx,gx,tah,coh,t,h

    integer, external         :: omp_get_thread_num,omp_get_num_threads

!if (node.eq.9999) then
!!$OMP PARALLEL
!    nodes = omp_get_num_threads()
!!$OMP END PARALLEL
!end if

!!$OMP PARALLEL DO PRIVATE         (inode,&
!tah,coh,t,h,it,iax,ix,ixm,iah,ih,izm,izd,&
!       indd,indm,minx,maxx,idx,rdx,fx,gx,&
!                 ax,tgx,xx,ah,tgh,hx,zm,zd) schedule(dynamic,1)  num_threads(nodes) 
                               
!    do ix=1, img%xm%n
!       inode = omp_get_thread_num()+1
!       minx=max(ix-nx,1)
!       maxx=min(ix+nx,sog%xm%n)
!       do iax=1,nax
!          ax = oax + (iax-1)*dax
!          tgx = -tand(ax)
!          do ixm=minx,maxx
!             xx=(ixm-ix+1)*sog%xm%d
!             do iah=1,nah
!                ah = oah + (iah-1)*dah
!                coh = -1/cosd(ah)
!                tah = tand(ah)
                do it=1,img%th%n
                   t=(it-1)*img%th%d+img%th%o
                   do ih=1,img%sh%n
                      h=(ih-1)*img%sh%d+img%sh%o
                      do izm=1, img%z%n
                         zm  = img%z%o + (izm-1) * img%z%d
!                         zd  = zm - coh*t*vv(izm,ix) - tah*h
                         if ((t*vv(izm,ix))**2 - h*h < 0. ) cycle
                         zd  = zm - sqrt((t*vv(izm,ix))**2 - h*h)
                         rdx = (zd - img%z%o) / img%z%d
                         idx = rdx
                         izd = idx + 1
                         if ( izd >= img%z%n .or. izd < 1) cycle
! if (mod(ix,20).eq.1) write(0,*) it,ih,izd,izm
!                         indd=(ixm-1)*img%z%n*img%th%n+(ih-1)*img%z%n+izd
!                         indm=(ix-1)*img%z%n*nah+(iax-1)*img%z%n*nah+(iah-1)*img%z%n+izm
!                         if (indm > size(modl) .or. indm < 1 .or. indd+1 > size(data) .or. indd+1 < 1) cycle
                         if (adj) then
                            modl(izm,ih,ix) = modl(izm,ih,ix) +  gx*data(izd,it,ix) +fx*data(izd+1,it,ix)
                         else
                            data(izd+1,it,ix)  = data(izd+1,it,ix)  + modl(izm,ih,ix)*fx
                            data(izd,  it,ix)  = data(izd,  it,ix)  + modl(izm,ih,ix)*gx
                         end if 
                      end do
                   end do
!                end do
!             end do
!          end do
       end do
       if (mod(ix,20).eq.1) write(0,*) "DONE CIG", ix, " FROM ", img%xm%n,data
!    end do
!!$OMP END PARALLEL DO
 end subroutine  
!===============================================================================
  function owslant_close() result(stat)
    integer                     :: stat 
    deallocate(vv)
    stat=0
  end function
end module 
