module OW_3drot2
! 3D-linear Radon transform

  use OW_type
  use OW_sltpar

  implicit none

  contains

!===============================================================================
  function ow3drot(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 ow3drot_exe(adj,add,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine ow3drot_exe(adj,add,modl,data)
    logical,intent(in)        :: adj,add
    real		  	:: modl(:) 
    real		  	:: data(:)

!! internal variables
    integer		:: inode
    integer                   :: ix,ixm,ih,iz,izm
    integer                   :: minx,maxx,minz,maxz,izd,indd,indm,idhx
    real			:: xx,hx,zm,zz,zd,tgx,tgh,tgs,dhx,dhz

    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,&
                       ix,ixm,ih,iz,izm,&
 minx,maxx,minz,maxz,izd,indd,indm,idhx,&
     xx,hx,zm,zz,zd,tgx,tgh,tgs,dhx,dhz) 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 ixm=minx,maxx
      xx=(ixm-ix+1)*sog%xm%d
      do ih=1,sog%sh%n
         hx=(ih-1)*sog%sh%d+sog%sh%o
         do iz=1,nz
            zm = zmin + (iz-1) * sog%z%d
            minz=max(iz-nzw,1)
            maxz=min(iz+nzw,sog%z%n)
            do izm=minz,maxz
               zz=(izm-iz+1)*sog%z%d
               if (xx.eq.0.0 .and. zz.ne.0.) cycle
               if (hx.eq.0.0 .and. zz.ne.0.) cycle
               if (xx.ne.0.0) then
                  tgx=-zz/xx
               else if (xx.eq.0.0 .and. zz.eq.0.0) then
                  tgx=-sog%xm%d/sog%z%d
               end if
               if (hx.ne.0.0) then
                  tgh=-zz/hx
               else if (hx.eq.0.0 .and. zz.eq.0.0) then
                  tgh=-sog%sh%d/sog%z%d
               end if
               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
               dhx=hx*sin(atan(tgx))*cos(atan(tgx))*tgh
               idhx = 1.5 + dhx/sog%sh%d
               dhz=hx*sin(atan(tgx))*sin(atan(tgx))*tgh
               if (rec) then
                  if (upshft) then
                     zd = zm + tgs*hx + dhz
                  else
                     zd = zm + tgs*hx
                  end if
               else
                  if (upshft) then
                     zd = zm - tgs*hx + dhz
                  else
                     zd = zm - tgs*hx
                  end if
               end if
               izd = 1.5 + (zd-zmin)/sog%z%d
               if ( izd > sog%z%n .or. izd < 1) cycle
!               if ( izd > maxz .or. izd < minz) cycle
               indm=(ix-1)*nz*sog%sh%n+(ih-1)*nz+iz
               if (upshft) then
                  indd=(ixm+idhx-1)*sog%z%n*sog%sh%n+(ih-1)*sog%z%n+izd
               else
                  indd=(ixm-1)*sog%z%n*sog%sh%n+(ih-1)*sog%z%n+izd
               end if
               if (indm > size(modl) .or. indm < 1 .or. indd > size(data) .or. indd < 1) cycle
               if (adj) then
                  modl(indm) = modl(indm) +  data(indd)!/abs(cos(atan(tgs)))
               else
                  data(indd) = data(indd) + modl(indm)!/abs(cos(atan(tgs)))
               end if 
            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
end module 
