module off2ang2d_sd_mod  
  ! 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
    
  integer function off2ang2d_sd_init(z0_in, hx0_in, dz_in, dhx_in, &
                                     nz_in, nhx_in, axmin_in, axmax_in, style_in) result(stat)
    
    
    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.

!write(0,*)"z0, dz, hx0, dhx, nz, nhx", z0, dz, hx0, dhx, nz, nhx


    allocate(tg(nax))

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


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

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

!===============================================================================
  integer function off2ang2d_sd_oper(adj, add, modl, data) result(stat)
    logical,intent(in)          :: adj,add 
    real,dimension(:,:), pointer:: modl,data 

    if (.not.add) then
       if (adj) then
          modl = 0.
       else
          data = 0.
       end if  
    end if
    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 






