module off2ang2d_sd_trans_pmod  
  ! simple linear Radon transform  z = z0 + p*h
  use adj_mod
  implicit none
  integer, private :: nz 
  integer, private :: nhx
  integer, private :: nax
  integer, private :: style
  real   , private :: z0, dz 
  real   , private :: hx0, dhx
  real   , private :: axmin, axmax
  real, dimension(:), pointer, private :: tg 
  real, parameter, private :: PI = 3.14159265359
contains
    
  subroutine off2ang2d_sd_init(z0_in, hx0_in, dz_in, dhx_in, &
                                     nz_in, nhx_in, axmin_in, axmax_in, style_in) 
    
    real    :: z0_in, dz_in 
    real    :: hx0_in, dhx_in 
    real    :: axmin_in, axmax_in
    integer :: nz_in 
    integer :: nhx_in 
    integer :: style_in
    integer :: iax
    real    :: dax, ax

    z0  = z0_in ; dz  = dz_in 
    hx0 = hx0_in; dhx = dhx_in 
    nz  = nz_in ; nhx = nhx_in
    nax = nhx_in
    style = style_in
    axmin = axmin_in * PI / 180.
    axmax = axmax_in * PI / 180.

    allocate(tg(nax))

    if (nax == 1) then
        dax = 0
    else
        dax = (axmax - axmin) / (nax - 1)
    end if


    do iax=1, nax
       ax = axmax - (iax-1)*dax
       tg(iax) = tan(ax)
    end do
 

  end subroutine
!===============================================================================
!  integer function off2ang2d_sd_clos() result(stat)
!  integer                     :: stat
!    deallocate(tg)
!    stat = 0
!  end function

!===============================================================================
  function off2ang2d_sd_oper(adj, add, modl, data) result(stat)
    integer                     :: stat 
    logical,intent(in)          :: adj,add 
    real,dimension(:)           :: modl,data 
    call adjnull (adj,add,modl,data )
    call off2ang2d_sd_compute(adj,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine off2ang2d_sd_compute(adj,modl,data)
    logical,intent(in)       :: adj 
    real, dimension(nz,nax)  :: modl 
    real, dimension(nz,nhx)  :: data 
    integer                  :: izm,izd
    integer                  :: ihx,iax
    integer                  :: stat, idx
    real                     :: zd, zm
    real                     :: hx
    real                     :: rdx, fx, gx
          
    do iax=1, nax 
!       if (abs(tg(iax)) > dz/dhx) cycle 
       do ihx=1, nhx
          hx = hx0 + (ihx-1) * dhx               ! offset
          do izm=1, nz  
             zm  = z0 + (izm-1) * dz             ! depth
             zd  = zm - tg(iax) * hx
             !! .. nearest neighbour interpolation
             if (style==0) then
                izd = 1.5 + (zd-z0)/dz
                if ( izd >nz .or. izd<1) cycle   ! out of the range
                if (adj) then
                   modl(izm,iax) = modl(izm,iax) + data(izd,ihx)
                else
                   data(izd,ihx) = data(izd,ihx) + modl(izm,iax)
                end if 
             !! .. linear interpolation
             else if(style==1) then
                rdx = (zd - z0) / dz
                idx = rdx
                izd = idx + 1
                if (izd >= nz .or. izd < 1) cycle
                fx = rdx - idx
                gx = 1.0 - fx
                if (adj) then
                   modl(izm, iax) = modl(izm, iax) &
                                  + data(izd, ihx)*gx + data(izd+1, ihx)*fx
                else
                   data(izd, ihx)   = data(izd, ihx)   + modl(izm, iax)*gx
                   data(izd+1, ihx) = data(izd+1, ihx) + modl(izm, iax)*fx
                end if
             end if
                
          end do
       end do
    end do
    
  end subroutine 
  
  
end module