module OW_tim2ang
! 3D-linear Radon transform

  use igrad1
  use OW_type
  use OW_parms_areal

  implicit none

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

  contains

!===============================================================================
  function owslant_init(v_in) result(stat)
  real,optional          :: v_in(img%xm%n,img%z%n)
  integer                :: i,j,stat

write(0,*) size(v_in,1),size(v_in,2)
  allocate(vv(img%xm%n,img%z%n))
  do j=1,img%z%n
  do i=1,img%xm%n
     vv(i,j)=sum(v_in(i,1:j))/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(:,:,:) 
    real                        :: data(:,:,:)

    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,ixm,iah,ih,izm,izd,it,ix
    integer		:: minx,maxx,indd,indm,idx,nodes
    real			:: ax,tgx,xx,ah,tgh,hx,zm
    real                      :: sech,secx,zd,rdx,fx,gx,tah,coh,t,h
    real,allocatable	:: temp(:)

    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
else
    nodes=node
end if
if (not(adj)) then
   allocate(temp(img%th%n))
   do izm=1,img%z%n
      do ix=1,img%xm%n
         call igrad1_lop2(.false.,.false.,temp,data(ix,:,izm))
         data(ix,:,izm)=temp
      end do
   end do
   deallocate(temp)
end if

!$OMP PARALLEL DO PRIVATE         (inode,&
tah,coh,t,h,it,iax,ix,ixm,iah,ih,izm,izd,&
       sech,secx,                        &
       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-nxw,1)
       maxx=min(ix+nxw,img%xm%n)
       do ixm=minx,maxx
          do iax=1,nax
             ax = oax + (iax-1)*dax
             secx = -1/cosd(ax)
             do iah=1,nah
                ah = oah + (iah-1)*dah
                sech = -1/cosd(ah)
                do it=1,img%th%n
                   t=(it-1)*img%th%d+img%th%o
                   do izm=1, img%z%n
                      zm  = img%z%o + (izm-1) * img%z%d
                      zd  = zm - t/vv(ix,izm)*sech*secx
                      rdx = (zd - img%z%o) / img%z%d
                      idx = rdx
                      izd = idx + 1
                      if ( izd >= img%z%n .or. izd < 1) cycle
                      fx = rdx - idx
                      gx = 1.0 - fx
                      if (adj) then
                         modl(ix,iah,izm) = modl(ix,iah,izm) +  gx*data(ixm,it,izd) +fx*data(ixm,it,izd+1)
                      else
                         data(ixm,it,izd+1)  = data(ixm,it,izd+1)  + modl(ix,iah,izm)*fx
                         data(ixm,it,izd  )  = data(ixm,it,izd  )  + modl(ix,iah,izm)*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 ", img%xm%n
     end do
!$OMP END PARALLEL DO
     if (adj) then
        allocate(temp(nah))
        do izm=1,img%z%n
           do ix=1,img%xm%n
              call igrad1_lop2(.true.,.false.,temp,modl(ix,:,izm))
              modl(ix,:,izm)=temp
           end do
        end do
        deallocate(temp)
     end if
 end subroutine  
!===============================================================================
  function owslant_close() result(stat)
    integer                     :: stat 
    deallocate(vv)
    stat=0
  end function
end module 
