module off2ang2d2_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 :: mode
  real   , private :: z0, dz 
  real   , private :: hx0, dhx
  real   , private :: axmin, dax
  real, dimension(:), pointer, private :: tg,co
  real, parameter, private :: PI = 3.14159265359
contains
    
  subroutine off2ang2d1_sd_init(z0_in, hx0_in, dz_in, dhx_in, nz_in, nhx_in, nax_in, dax_in, axmin_in, mode_in) 
    
    real    :: z0_in, dz_in 
    real    :: hx0_in, dhx_in 
    real    :: axmin_in, dax_in
    integer :: nz_in
    integer :: nhx_in, nax_in
    integer :: mode_in
    integer :: iax
    real    :: ax

    z0  = z0_in ; dz  = dz_in 
    hx0 = hx0_in; dhx = dhx_in 
    nz  = nz_in ; nhx = nhx_in
    nax = nax_in; dax = dax_in 
    mode = mode_in
    axmin = axmin_in

    allocate(tg(nax),co(nax))

    do iax=1, nax
       ax = axmin + (iax-1)*dax
       tg(iax) = -tand(ax)
       co(iax)=abs(cosd(ax))
    end do

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

!===============================================================================
  function off2ang2d1_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 off2ang2d1_sd_compute(adj,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine off2ang2d1_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                     	:: calfa
    real                     	:: rdx, fx, gx

    calfa=cos(30*3.14159/180)
    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)/calfa**2 * hx
             !! .. linear interpolation
             rdx = (zd - z0) / dz
             idx = rdx
             izd = idx + 1
             if (izd > nz .or. izd < 1) cycle				   ! out of the range
             fx = rdx - idx
             gx = 1.0 - fx
             if (mode==0) then						   ! uses the cosine weight
                if ( izd >nz .or. izd<1) cycle					   ! out of the range
                if (adj) then
                    modl(izm+nz*(iax-1)) = modl(izm+nz*(iax-1)) &
                                  + 1/co(iax)*(data(izd+nz*(ihx-1))*gx + data(izd+1+nz*(ihx-1))*fx)
                else
                   data(izd+nz*(ihx-1)) = data(izd+nz*(ihx-1)) + 1/co(iax)*modl(izm+nz*(iax-1))*gx
                   data(izd+1+nz*(ihx-1)) = data(izd+1+nz*(ihx-1)) + 1/co(iax)*modl(izm+nz*(iax-1))*fx
                end if 
             else if(mode==1) then					   ! does not use the cosine weight
                if (adj) then
!----------------------------
! Passes dotprod test
                    modl(izm+nz*(iax-1)) = modl(izm+nz*(iax-1)) &
                                  + (data(izd+nz*(ihx-1))*gx + data(izd+1+nz*(ihx-1))*fx)
!----------------------------
! dot prod bad
!                   modl(iax+nax*(izm-1)) = modl(iax+nax*(izm-1)) &
!                                  + 1/co(iax)*(data(ihx+nhx*(izd-1))*gx + data(ihx+nhx*(izd))*fx)
!                    modl(iax+nax*(izm-1)) = modl(iax+nax*(izm-1)) &
!                                  + (data(ihx+nhx*(izd-1))*gx + data(ihx+nhx*(izd))*fx)
               else
!----------------------------
! Passes dotprod test
                   data(izd+nz*(ihx-1)) = data(izd+nz*(ihx-1)) + modl(izm+nz*(iax-1))*gx
                   data(izd+1+nz*(ihx-1)) = data(izd+1+nz*(ihx-1)) + modl(izm+nz*(iax-1))*fx
!----------------------------
! dot prod bad
!                   data(ihx+nhx*(izd-1)) = data(ihx+nhx*(izd-1)) + 1/co(iax)*modl(iax+nax*(izm-1))*gx
!                   data(ihx+nhx*(izd)) = data(ihx+nhx*(izd)) + 1/co(iax)*modl(iax+nax*(izm-1))*fx
!                   data(ihx+nhx*(izd-1)) = data(ihx+nhx*(izd-1)) + modl(iax+nax*(izm-1))*gx
!                   data(ihx+nhx*(izd)) = data(ihx+nhx*(izd)) + modl(iax+nax*(izm-1))*fx
                end if
             end if
          end do
       end do
    end do
  end subroutine 
  
end module