module OW1_3dslt
! 3D-linear Radon transform

  use OW_type
  use OW_sltpar

  implicit none

  contains

!===============================================================================
  function ow3dslant(adj, add, modl, data) result(stat)
    integer                     :: stat 
    logical,intent(in)          :: adj,add 
    real,dimension(:)           :: modl
    real,dimension(:)           :: data

    if (adj) then
       modl=0.
    else
       data=0.
    end if
    call ow3dslant_exe(adj,add,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine ow3dslant_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
    integer		:: minx,maxx,indd,indm,idx
    real			:: ax,tgx,xx,ah,tgh,hx,zm,zd,rdx,fx,gx

    integer, external         :: omp_get_thread_num,omp_get_num_threads

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

!$OMP PARALLEL DO PRIVATE  (inode,&
        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, sog%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
                tgh = -tand(ah)
                if (rec) then
                   if (1-tgx*tgh.le.tiny(1.)) then 
                      cycle
                   else
                      tgs=(tgx+tgh)/(1-tgx*tgh)
                   end if
                else
                   if (1+tgx*tgh.le.tiny(1.)) then
                      cycle
                   else
                      tgs=(tgx-tgh)/(1+tgx*tgh)
                   end if
                end if
                do ih=1,sog%sh%n
                   hx=(ih-1)*sog%sh%d+sog%sh%o
                   do izm=1, sog%z%n
                      zm  = sog%z%o + (izm-1) * sog%z%d
                      if (rec) then
                         zd  = zm + tgs*hx
                      else
                         zd  = zm - tgs*hx
                      end if
                      rdx = (zd - sog%z%o) / sog%z%d
                      idx = rdx
                      izd = idx + 1
!                      izd = 1.5 + (zd-sog%z%o)/sog%z%d
                      if ( izd > sog%z%n .or. izd < 1) cycle
                      indd=(ixm-1)*sog%z%n*sog%sh%n+(ih-1)*sog%z%n+izd
!                      indm=(ix-1)*sog%z%n*nah*nax+(iax-1)*sog%z%n*nah+(iah-1)*sog%z%n+izm
                      indm=(ix-1)*sog%z%n*nah+(iax-1)*sog%z%n*nah+(iah-1)*sog%z%n+izm
                      if (indm > size(modl) .or. indm < 1 .or. indd+1 > size(data) .or. indd+1 < 1) cycle
!                         if (adj) then
!                            modl(indm) = modl(indm) +  data(indd)
!                         else
!                            data(indd)  = data(indd)  + modl(indm)
!                         end if 
                      if (adj) then
                         modl(indm) = modl(indm) +  gx*data(indd) +fx*data(indd+1)
                      else
                         data(indd+1)  = data(indd+1)  + modl(indm)*fx
                         data(indd  )  = data(indd  )  + modl(indm)*gx
                      end if 
                   end do
                end do
             end do
          end do
       end do
       if (mod(ix,20).eq.1) write(0,*) "DONE CIG", ix, " FROM ", sog%xm%n
    end do
!$OMP END PARALLEL DO
 end subroutine  

!                      indd=(ix-1)*sog%z%n*sog%sh%n+(ih-1)*sog%z%n+izd
!                      indm=(ix-1)*sog%z%n*nah*nax+(iax-1)*sog%z%n*nah+(iah-1)*sog%z%n+izm

!
!=================================================================================  
  integer function ident_lop( adj, add, model, data)
    real, dimension(:)            :: model, data
    logical                       :: adj,add

    if (adj) then
      model=0
      model = model + data
    else
      data=0
      data  = data + model
    end if 
    ident_lop = 0
  end function
!=================================================================================  

end module 
