module OW_3diprad
! 3D-linear Radon transform

  use OW_type
  use OW_sltpar

  implicit none

  real, allocatable, private, save	:: dips(:)

  contains

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

!! internal variables
    integer		:: inode
    integer 		:: iax,ix,ixd,iah,ih,izm,izd
    integer		:: minx,maxx,indd,indm,idx
    real			:: ax,tgx,xx,ah,tgh,hx,zm,zd,rdx,fx,gx

    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,&
               iax,ix,ixd,izm,izd,&
              indd,indm,minx,maxx,&
                    idx,rdx,fx,gx,&
                    ax,tgx,xx,zm,zd) schedule(dynamic,1)  num_threads(nodes) 
                               
    do ix=1, sog%xm%n
       inode = omp_get_thread_num()+1
       minx=max(ix-nx,1)
       maxx=min(ix+nx,sog%xm%n)
       do iax=1,nax
          ax = oax + (iax-1)*dax
          tgx = -tand(ax)
          do ixd=minx,maxx
             xx=(ixd-ix+1)*sog%xm%d
             do izm=1, sog%z%n
                zm  = sog%z%o + (izm-1) * sog%z%d
                zd  = zm - tgx*xx
                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=(ixd-1)*sog%z%n+izd
                indm=(ix-1)*sog%z%n*nax+(iax-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
       if (mod(ix,20).eq.1) write(0,*) "DONE CIG", ix, " FROM ", sog%xm%n
    end do
!$OMP END PARALLEL DO
 end subroutine  

!
!===============================================================================
  function ow3curvrad_init(dips_in) result(stat)
    integer                     :: stat 
    real,dimension(:)           :: dips_in
   
    allocate(dips(size(dips_in))); dips=dips_in
    stat =0 

  end function
!===============================================================================
  function ow3curvrad_close() result(stat)
    integer                     :: stat 
   
    deallocate(dips)
    stat=0 
  end function
!===============================================================================
  function ow3curvrad(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 ow3curvrad_exe(adj,add,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine ow3curvrad_exe(adj,add,modl,data)
    logical,intent(in)        :: adj,add
    real		  	:: modl(:) 
    real		  	:: data(:)

!! internal variables
    integer		:: inode
    integer 		:: ixm,ir,it,izm,izd,ixd
    integer		:: indd,indm
    real			:: xm,radius,t,zm,a,b,xd,zd,ang1,ang2,pi,fac
    integer, external         :: omp_get_thread_num,omp_get_num_threads
    integer,allocatable	:: weight(:,:,:,:,:)

pi = atan(1.)*4.
if (nodes.eq.9999) then
!$OMP PARALLEL
    nodes = omp_get_num_threads()
!$OMP END PARALLEL
end if

if (applw) then
allocate(weight(sog%z%n,nrad,sog%xm%n,sog%z%n,sog%xm%n));weight=1;
!$OMP PARALLEL DO PRIVATE  (inode,&
              ixm,xm,ir,radius,it,&
    fac,ang1,ang2,t,izm,zm,a,b,xd,&
              ixd,zd,izd,indd,indm) schedule(dynamic,1)  num_threads(nodes) 
                               
    do ixm=1, sog%xm%n
       inode = omp_get_thread_num()+1
       xm = sog%xm%o + (ixm-1) * sog%xm%d
       do ir=1,nrad
          radius=(ir-1)*drad+orad
          do izm=1, sog%z%n
             ang1=-dips((ixm-1)*sog%z%n+izm)
             ang2=90
             if (ang1.gt.tiny(1.)) then
                ang2=-atand(1/tand(ang1))
             end if
             zm  = sog%z%o + (izm-1)        * sog%z%d
             a   = xm      - sind(ang2)    * radius
             b   = zm      - cosd(ang2)    * radius
             do it=1,nts
                t   = ots     + (it-1)*dts
                xd  = a       - sind(ang2+t)  * radius
                zd  = b       - cosd(ang2+t)  * radius
                ixd = 1.5     + (xd - sog%xm%o)/sog%xm%d
                izd = 1.5     + (zd - sog%z%o)/sog%z%d
                if ( izd > sog%z%n .or. izd < 1 .or. ixd > sog%xm%n .or. ixd < 1) cycle
                indd=(ixd-1)*sog%z%n+izd
                indm=(ixm-1)*sog%z%n*nrad+(ir-1)*sog%z%n+izm
                if (indm > size(modl) .or. indm < 1 .or. indd > size(data) .or. indd < 1) cycle
                weight(izm,ir,ixm,izd,ixd)=weight(izm,ir,ixm,izd,ixd)+1
             end do
          end do
       end do
       if (mod(ixm,20).eq.1) write(0,*) " Weights ==> DONE CIG", ixm, " FROM ", sog%xm%n
    end do
!$OMP END PARALLEL DO
end if

!$OMP PARALLEL DO PRIVATE  (inode,&
              ixm,xm,ir,radius,it,&
    fac,ang1,ang2,t,izm,zm,a,b,xd,&
              ixd,zd,izd,indd,indm) schedule(dynamic,1)  num_threads(nodes) 

    do ixm=1, sog%xm%n
       inode = omp_get_thread_num()+1
       xm = sog%xm%o + (ixm-1) * sog%xm%d
       do ir=1,nrad
          radius=(ir-1)*drad+orad
          do izm=1, sog%z%n
             ang1=dips((ixm-1)*sog%z%n+izm)
             zm  = sog%z%o + (izm-1)        * sog%z%d
             a   = xm      + cosd(ang1)    * radius
             b   = zm      + sind(ang1)    * radius
             do it=1,nts
                t   = ots     + (it-1)*dts
                xd  = a       + sind(ang1+t)  * radius
                zd  = b       + cosd(ang1+t)  * radius
                ixd = 1.5     + (xd - sog%xm%o)/sog%xm%d
                izd = 1.5     + (zd - sog%z%o)/sog%z%d
                if ( izd > sog%z%n .or. izd < 1 .or. ixd > sog%xm%n .or. ixd < 1) cycle
                indd=(ixd-1)*sog%z%n+izd
                indm=(ixm-1)*sog%z%n*nrad+(ir-1)*sog%z%n+izm
                if (indm > size(modl) .or. indm < 1 .or. indd > size(data) .or. indd < 1) cycle
                if (adj) then
                    if (applw) then
                       modl(indm) = modl(indm) + data(indd)/weight(izm,ir,ixm,izd,ixd)
                    else
                       modl(indm) = modl(indm) + data(indd)
                    end if
                else
                    if (applw) then
                       data(indd) = data(indd) + modl(indm)/weight(izm,ir,ixm,izd,ixd)
                    else
                       data(indd) = data(indd) + modl(indm)
                    end if
                end if 
             end do
          end do
       end do
       if (mod(ixm,20).eq.1) write(0,*) " Curvature ==> DONE CIG", ixm, " FROM ", sog%xm%n
    end do
!$OMP END PARALLEL DO

if (exist_file("weight")) then
     call to_history("n1",sog%z%n,"weight");call to_history("d1",sog%z%d,"weight");call to_history("o1",sog%z%o,"weight")
     call to_history("n2",nrad,"weight");call to_history("d2",drad,"weight");call to_history("o2",orad,"weight")  
     call to_history("n3",sog%xm%n,"weight");call to_history("d3",sog%xm%d,"weight");call to_history("o3",sog%xm%o,"weight")  
     call to_history("n4",sog%z%n,"weight");call to_history("d4",sog%z%d,"weight");call to_history("o4",sog%z%o,"weight")
     call to_history("n5",sog%xm%n,"weight");call to_history("d5",sog%xm%d,"weight");call to_history("o5",sog%xm%o,"weight")  
     call to_history("label1","z","weight");call to_history("label2","x","weight")
     call to_history("esize",4,"weight")
     do ixd=1, sog%xm%n
        do izd=1, sog%z%n
           do ixm=1, sog%xm%n
              do ir=1,nrad
                 call srite("weight",real(weight(:,ir,ixm,izd,ixd)),4*sog%z%n)
              end do
           end do
        end do
     end do
end if

 end subroutine  
!
!       minx=max(ix-nx,1)
!       maxx=min(ix+nx,sog%xm%n)
!                zd  = zm - tgx*xx
!                rdx = (zd - sog%z%o) / sog%z%d
!                idx = rdx
!                izd = idx + 1
!                fx = rdx - idx
!                gx = 1.0 - fx

!
!
!=================================================================================  
  integer function ident_lop( adj, add, model, data)
    real, dimension(:)            :: model, data
    logical                       :: adj,add

    if (adj) then
      model=0
      model = model + data
    else
      data=0
      data  = data + model
    end if 
    ident_lop = 0
  end function
!=================================================================================  

end module 
