module OW_3dsltrt
! 3D-linear Radon transform

  use OW_type
  use OW_sltpar

  implicit none
  real, allocatable, private  :: dips(:,:)

  contains

!===============================================================================
  function ow3dslantrt_init(dips_in) result(stat)
    integer                     :: stat
    real,optional               :: dips_in(:,:)

    if (present(dips_in)) then
       allocate(dips(sog%z%n,sog%xm%n))
       dips=dips_in
    end if
    stat=0
  end function 
!===============================================================================
  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
!         w=1-abs(float(ix-ixm)/float(nxw));w=w**wexp
         if (ix-ixm >= 0) then
            w=1-(float(ix-ixm)/float(nxw))**wexp
         else
            w=1-(float(ixm-ix)/float(nxw))**wexp
         end if
         do iah=1,nah
            ah = oah + (iah-1)*dah
            tgh = -tand(ah)
            do ih=1,sog%sh%n
               hx=(ih-1)*sog%sh%d+sog%sh%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
               hw=1 - abs(hx)/2/abs(sog%sh%o)
               hw=hw**4
               dhx=abs(hx)*sind(ax)*cosd(ax)*abs(tgh)
               idhx = 1.5 + dhx/sog%sh%d
               dhz=abs(hx)*sind(ax)*sind(ax)*abs(tgh)
               do izm=1,sog%z%n
                  if (allocated(dips)) then
!                     y = (1 - abs(dips(izm,ix)-ax)/2/abs(oax))
                     if (dips(izm,ix) - ax >=0) then
                        y = 1 - ((dips(izm,ix)-ax)/abs(oax))**yexp
                     else
                        y = 1 - ((ax-dips(izm,ix))/abs(oax))**yexp
                     end if
                     if (y < 0) y=0.
                  else
                     y = 1           
                  end if
                  zm  = zmin + (izm-1) * sog%z%d
                  if (rec) then
                     if (upshft) then
                        zd  = zm - tgh*hx - tgx*xx + tgs*hx - dhz
                     else
                        zd  = zm - tgh*hx - tgx*xx + tgs*hx
                     end if
                  else
                     if (upshft) then
                        zd  = zm - tgh*hx - tgx*xx - tgs*hx - dhz
                     else
                        zd  = zm - tgh*hx - tgx*xx - tgs*hx
                     end if
                  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%sh%n+(ih-1)*sog%z%n+izd
                  indm=(ix -1)*sog%z%n*sog%sh%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))*y*w!/hw
                  else
                     wdata(indm) = wdata(indm) + wmodl(indd)
                     data(indd+1)  = data(indd+1)  + modl(indm)*fx*y*w!/hw
                     data(indd  )  = data(indd  )  + modl(indm)*gx*y*w!/hw
                  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))
      if (exist_file("wmodl")) then
         call to_history("n1",nz,"wmodl");call to_history("d1",sog%z%d,"wmodl");call to_history("o1",zmin,"wmodl")
         call to_history("n2",sog%sh%n,"wmodl");call to_history("d2",sog%sh%d,"wmodl");call to_history("o2",sog%sh%o,"wmodl")
         call to_history("n3",sog%xm%n,"wmodl");call to_history("d3",sog%xm%d,"wmodl");call to_history("o3",sog%xm%o,"wmodl")  
         call to_history("label1","z","wmodl");call to_history("label2","ah","wmodl")
         call to_history("label3","x","wmodl")
         call to_history("esize",4,"wmodl")
         call srite("wmodl",wmodl,4*nz*sog%sh%n*sog%xm%n)
      end if
   else
      data=data/(wdata/maxval(wdata))
      if (exist_file("wdata")) then
         call to_history("n1",nz,"wdata");call to_history("d1",sog%z%d,"wdata");call to_history("o1",zmin,"wdata")
         call to_history("n2",sog%sh%n,"wdata");call to_history("d2",sog%sh%d,"wdata");call to_history("o2",sog%sh%o,"wdata")
         call to_history("n3",sog%xm%n,"wdata");call to_history("d3",sog%xm%d,"wdata");call to_history("o3",sog%xm%o,"wdata")  
         call to_history("label1","z","wdata");call to_history("label2","ah","wdata")
         call to_history("label3","x","wdata")
         call to_history("esize",4,"wdata")
         call srite("wdata",wdata,4*nz*sog%sh%n*sog%xm%n)
      end if
   end if
   deallocate(wmodl,wdata)
 end subroutine  
!===============================================================================
  function ow3dslantrt_close() result(stat)
    integer             :: stat
    deallocate(dips)
    stat=0
  end function 
!=================================================================================  
end module 
