module OW_3dsltrt1
! 3D-linear Radon transform

  use OW_type
  use OW_sltpar

  implicit none
  integer, private                     :: iy,miny, maxy
  contains

!===============================================================================
  function ow3dslantrt(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 ow3dslantrt_exe(adj,add,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine ow3dslantrt_exe(adj,add,modl,data)
    logical,intent(in)        :: adj,add
    real		  	:: modl(:) 
    real		  	:: data(:)
    real,allocatable          :: wmodl(:),wdata(:)

!! internal variables
    integer		:: inode
    integer 		:: iax,ix,ixm,iah,ih,izm,izd,idhx
    integer		:: minx,maxx,indd,indm,idx
    real			:: ax,tgx,xx,ah,tgh,hx,dhx,dhz
    real                      :: hw,zm,zd,rdx,fx,gx,w,y,tgs

    integer, external         :: omp_get_thread_num,omp_get_num_threads

   allocate(wmodl(size(modl)))
   wmodl=1.
   allocate(wdata(size(data)))
   wdata=1.

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

!$OMP PARALLEL DO PRIVATE  (inode,dhx,&
idhx,dhz,hw,iax,ix,ixm,iah,ih,izm,izd,&
    indd,indm,minx,maxx,idx,rdx,fx,gx,&
     tgs,y,w,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-nxw,1)
   maxx=min(ix+nxw,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,nahx
            ah = oahx + (iah-1)*dahx
            tgh = -tand(ah)
            do ih=1,sog%shx%n
               hx=(ih-1)*sog%shx%d+sog%shx%o
               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 izm=1,sog%z%n
                  zm  = zmin + (izm-1) * sog%z%d
                  if (rec) then
                     zd  = zm - tgh*hx - tgx*xx + tgs*hx
                  else
                     zd  = zm - tgh*hx - tgx*xx - tgs*hx
                  end if
                  rdx = (zd - sog%z%o) / sog%z%d
                  idx = rdx
                  izd = idx + 1
                  fx = rdx - idx
                  gx = 1.0 - fx
                  if ( izd > sog%z%n .or. izd < 1) cycle
                  indd=(ixm-1)*sog%z%n*sog%shx%n+(ih-1)*sog%z%n+izd
                  indm=(ix -1)*sog%z%n*sog%shx%n+(ih-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
                     wmodl(indm) = wmodl(indm) + wdata(indd)
                     modl(indm) = modl(indm) +  (gx*data(indd) +fx*data(indd+1))
                  else
                     wdata(indm) = wdata(indm) + wmodl(indd)
                     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
   if (adj) then
      modl=modl/(wmodl/maxval(wmodl))
   else
      data=data/(wdata/maxval(wdata))
   end if
   deallocate(wmodl,wdata)
 end subroutine  
!!===============================================================================
!  function ow3dslantrt25_init(iy_in,miny_in,maxy_in) result(stat)
!    integer                     :: iy_in,miny_in,maxy_in
!    iy=iy_in
!    maxy=maxy_in
!    miny=miny_in
!    stat=0
!  end function
!!===============================================================================
!  function ow3dslantrt25(adj, add, modl, data) result(stat)
!    integer                     :: stat 
!    logical,intent(in)          :: adj,add 
!    real,dimension(:)           :: modl
!    real,dimension(:)           :: data
!    integer, external           :: omp_get_thread_num,omp_get_num_threads
!
!    if (nodes.eq.9999) then
!    !$OMP PARALLEL
!        node = omp_get_num_threads()
!    !$OMP END PARALLEL
!    else
!        node=nodes
!    end if
!
!    if (adj) then
!       modl=0.
!    else
!       data=0.
!    end if
!    call ow3dslantrt2_exe(adj,add,modl,data )
!    stat=0
!  end function 
!=================================================================================  
!  subroutine ow3dslantrt2_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,idhx
!    integer		:: minx,maxx,indd,indm,idx
!    integer                   :: iay,iym,iahy,ihy
!    real			:: ax,tgx,xx,ah,tgh,hx,dhx,dhz
!    real                      :: hw,zm,zd,rdx,fx,gx,w,y,tgs
!    real                      :: ay,tgy,hy,tgsy,ahy,tghy,yy
!    real                      :: starty,finishy,startx,finishx
!
!    integer, external         :: omp_get_thread_num,omp_get_num_threads
!
!!$OMP PARALLEL DO PRIVATE  (inode,dhx,&
!idhx,dhz,hw,iax,ix,ixm,iah,ih,izm,izd,&
!    indd,indm,minx,maxx,idx,rdx,fx,gx,&
!    tgs,y,w,ax,tgx,xx,ah,tgh,hx,zm,zd,&
!    iay,ay,tgy,iym,yy,hy,&
!               iahy,tghy,ihy,ahy,tgsy) schedule(dynamic,1) ordered num_threads(node) 
!                               
!   do iay=1,nay
!      inode = omp_get_thread_num()+1
!      ay = oay + (iay-1)*day
!      tgy = -tand(ay)
!      do iym=miny,maxy
!         yy=(iym-iy+1)*sog%ym%d
!         do iahy=1,nahy
!            ahy = oahy + (iahy-1)*dahy
!            tghy = -tand(ahy)
!            do ihy=1,sog%shy%n
!               hy=(ihy-1)*sog%shy%d+sog%shy%o
!               if (rec) then
!                  if (1-tgy*tghy.le.tiny(1.)) then 
!                     cycle
!                  else
!                     tgsy=(tgy+tghy)/(1-tgy*tghy)
!                  end if
!               else
!                  if (1+tgy*tghy.le.tiny(1.)) then
!                     cycle
!                  else
!                     tgsy=(tgy-tghy)/(1+tgy*tghy)
!                  end if
!               end if
!               do ix=1, sog%xm%n
!                  call owexec_time (startx)
!                  minx=max(ix-nxw,1)
!                  maxx=min(ix+nxw,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,nahx
!                           ah = oahx + (iah-1)*dahx
!                           tgh = -tand(ah)
!                           do ih=1,sog%shx%n
!                              hx=(ih-1)*sog%shx%d+sog%shx%o
!                              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 izm=1,sog%z%n
!                                 zm  = zmin + (izm-1) * sog%z%d
!                                 if (rec) then
!                                    zd  = zm - tgh*hx - tgx*xx + tgs*hx - tghy*hy - tgy*yy + tgsy*hy
!                                 else
!                                    zd  = zm - tgh*hx - tgx*xx - tgs*hx - tghy*hy - tgy*yy - tgsy*hy
!                                 end if
!                                 rdx = (zd - sog%z%o) / sog%z%d
!                                 idx = rdx
!                                 izd = idx + 1
!                                 fx = rdx - idx
!                                 gx = 1.0 - fx
!                                 if ( izd > sog%z%n .or. izd < 1) cycle
!                                 indd=(iym-1)*sog%z%n*sog%shx%n*sog%xm%n*sog%shy%n+(ihy-1)*sog%z%n*sog%shx%n*sog%xm%n+&
!                                      (ixm-1)*sog%z%n*sog%shx%n+(ih-1)*sog%z%n+izd
!                                 indm=(iy -1)*sog%z%n*sog%shx%n*sog%xm%n*sog%shy%n+(ihy-1)*sog%z%n*sog%shx%n*sog%xm%n+&
!                                      (ix -1)*sog%z%n*sog%shx%n+(ih-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) +  (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
!               end do
!            end do
!         end do
!      end do
!   end do
!!$OMP END PARALLEL DO
!
! end subroutine  
!!=================================================================================  
end module 
