module OW_3drot_det
! 3D-linear Radon transform

  use OW_type
  use OW_sltpar

  integer, private		:: oxchunk,xblk
  real, allocatable, private  :: dips(:,:)
  contains

!===============================================================================
  function ow3drotdet_init(oxchunk_in,xblk_in,dips_in) result(stat)
    real                        :: dips_in(:,:)
    integer                     :: oxchunk_in,xblk_in,stat
    oxchunk=oxchunk_in
    xblk=xblk_in
    allocate(dips(nz,xblk))
    dips=dips_in
    stat=0
  end function 
!===============================================================================
  function ow3drotdet(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 ow3drotdet_exe(adj,add,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine ow3drotdet_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,idx
    integer		:: minx,maxx,indd,indm,idhx,minax,maxax
    real			:: ax,axx,tgx,xx,ah,tgh,hx,zm,zd,tgs,rdx,fx,gx,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,&
minax,maxax,iax,ix,iah,ih,izm,izd,&
  axx,idx,indd,indm,tgs,rdx,fx,gx,&
dhx,idhx,dhz,ax,tgx,ah,tgh,hx,zm,zd) schedule(dynamic,1)  num_threads(nodes) 
                               
do ix=1,xblk
   inode = omp_get_thread_num()+1
   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
         do izm=1,nz !nzmin, nzmax
            zm  = zmin + (izm-1) * sog%z%d
            axx=float(int(dips(izm,ix)))-float(mod(int(dips(izm,ix)),2))
            minax=int((axx-20-oax)/dax)+1.5
            maxax=int((axx+20-oax)/dax)+1.5
            do iax=minax,maxax
!            iax=int((ax-oax)/dax) + 1.5
               ax=(iax-minax)*dax+(axx-20)
               ax=axx
               tgx=-tand(ax)
               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)
               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
               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
!!.. linear interpolation
!               rdx = (zd - sog%z%o) / sog%z%d
!               idx = rdx
!               izd = idx + 1
!               if (izd > sog%z%n.or. izd < 1) cycle
!               fx = rdx - idx
!               gx = 1.0 - fx
!               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 
!....................
!               izd = 1.5 + (zd-sog%z%o)/sog%z%d
!               if ( izd > sog%z%n .or. izd < 1) cycle
               izd = 1.5 + (zd-zmin)/sog%z%d
               if ( izd > nz .or. izd < 1) cycle
               indm=(ix-1)*nz*sog%sh%n+(ih-1)*nz+izm
               if (upshft) then
                  indd=(ix+idhx-1)*nz*sog%sh%n*nax*nah+(iax-1)*nz*sog%sh%n*nah+(iah-1)*nz*sog%sh%n+(ih-1)*nz+izd
               else
                  indd=(ix-1)*nz*nax*nah+(iax-1)*nz*nah+(iah-1)*nz+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
end do

!$OMP END PARALLEL DO
 write(0,*) "DONE ROTATION UP TO CIG", oxchunk-1+xblk, " FROM ", sog%xm%n
 end subroutine
end module 
