module wei_off2ang2d_sd_mod
  ! simple linear Radon transform  z = z0 + tan(gamma)h
  ! input SODCIG(z,x,h) => ADCIG(z,x,gamma)
  use adj_mod
  implicit none
  integer                    , private :: nz, nh, nang
  real                       , private :: oz, oh, oang
  real                       , private :: dz, dh, dang
  real, dimension(:), pointer, private :: tg_ang
  integer                    , private :: style
  real, parameter            , private :: PI = 3.14159265359
contains
    
  integer function wei_off2ang2d_sd_init(az_in, ah_in, aang_in, style_in) result(stat)
    use wei_type_mod
    type(axis) :: az_in, ah_in, ax_in, aang_in
    integer    :: style_in
    integer    :: iang
    real       :: ang

    nz   = az_in%n  ; oz   = az_in%o  ; dz   = az_in%d
    nh   = ah_in%n  ; oh   = ah_in%o  ; dh   = ah_in%d
    nang = aang_in%n; oang = aang_in%o; dang = aang_in%d

    style = style_in

    allocate(tg_ang(nang))

    do iang=1, nang
       ang = oang + (iang-1)*dang
       ang = ang * PI / 180.
       tg_ang(iang) = tan(ang)
    end do

    stat = 0
  end function
!===============================================================================
  integer function wei_off2ang2d_sd_clos() result(stat)
    deallocate(tg_ang)
    stat = 0
  end function

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

    if (.not.add) then
       if (adj) then
          modl = 0.
       else
          data = 0.
       end if  
    end if
    call wei_off2ang2d_sd_compute(adj, modl, data)
    stat=0
  end function 
!=================================================================================  
  subroutine wei_off2ang2d_sd_compute(adj, modl, data)
    logical,intent(in)              :: adj
    real, dimension(:,:)  , pointer :: modl 
    real, dimension(:,:)  , pointer :: data 
    integer                         :: iang, ih, ix, izm, izd
    integer                         :: idx
    real                            :: h, x, zd, zm
    real                            :: rdx, fx, gx

    do iang=1, nang 
       !if (abs(tg(iax)) > dz/dhx) cycle 
!write(0,*)"idip, iang", idip, iang 
       do ih=1, nh
          h = oh + (ih-1)*dh               ! offset
          do izm=1, nz  
             zm  = oz + (izm-1) * dz             ! depth
             zd  = zm - tg_ang(iang) * h 
             !! .. nearest neighbour interpolation
             if (style==0) then
                izd = 1.5 + (zd-oz)/dz
                if ( izd >nz .or. izd<1) cycle   ! out of the range
                if (adj) then
                   modl(izm,iang) = modl(izm,iang) + data(izd,ih)
                else
                   data(izd,ih)   = data(izd,ih)   + modl(izm,iang)
                end if 
             !! .. linear interpolation
             else if(style==1) then
                rdx = (zd - oz) / 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,iang) = modl(izm,iang) &
                                  + data(izd, ih)*gx + data(izd+1, ih)*fx
                else
                   data(izd  , ih) = data(izd  , ih) + modl(izm,iang)*gx
                   data(izd+1, ih) = data(izd+1, ih) + modl(izm,iang)*fx
                end if
             end if
                
          end do
       end do
    end do
  end subroutine 
end module 






