module off2ang2d_sd_mod2
  ! simple linear Radon transform  z = z0 + tan(alpha)x + tan(gamma)h
  ! input SODCIG(z,x,h) => ADCIG(z,x,gamma,alpha)
  use adj_mod
  implicit none
  integer                    , private :: nz, nh, nx, nang, ndip
  real                       , private :: oz, oh, ox, oang, odip
  real                       , private :: dz, dh, dx, dang, ddip
  real, dimension(:), pointer, private :: tg_ang, tg_dip
  integer                    , private :: style
  real, parameter            , private :: PI = 3.14159265359
contains
    
  integer function off2ang2d_sd_init(az_in, ah_in, ax_in, aang_in, adip_in, style_in) result(stat)
    use wei_type_mod
    type(axis) :: az_in, ah_in, ax_in, aang_in, adip_in
    integer    :: style_in
    integer    :: iang, idip
    real       :: ang, dip

    nz   = az_in%n  ; oz   = az_in%o  ; dz   = az_in%d
    nh   = ah_in%n  ; oh   = ah_in%o  ; dh   = ah_in%d
    nx   = ax_in%n  ; ox   = ax_in%o  ; dx   = ax_in%d
    nang = aang_in%n; oang = aang_in%o; dang = aang_in%d
    ndip = adip_in%n; odip = adip_in%o; ddip = adip_in%d

    style = style_in

    allocate(tg_ang(nang), tg_dip(ndip))

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

    do idip=1, ndip
       dip = odip + (idip-1)*ddip
       dip = dip * PI / 180.
       tg_dip(idip) = tan(dip)
    end do

write(0,*)"az", az_in   
write(0,*)"ah", ah_in
write(0,*)"aang", aang_in
write(0,*)"adip", adip_in
    stat = 0
  end function
!===============================================================================
  integer function off2ang2d_sd_clos() result(stat)
    deallocate(tg_ang)
    deallocate(tg_dip)
    stat = 0
  end function

!===============================================================================
  integer function off2ang2d_sd_oper(adj, add, cmpx, modl, data) result(stat)
    logical,intent(in)              :: adj,add
    real                            :: cmpx
    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 off2ang2d_sd_compute(adj, cmpx, modl, data)
    stat=0
  end function 
!=================================================================================  
  subroutine off2ang2d_sd_compute(adj, cmpx, modl, data)
    logical,intent(in)              :: adj
    real                            :: cmpx
    real, dimension(:,:,:), pointer :: modl 
    real, dimension(:,:)  , pointer :: data 
    integer                         :: idip, iang, ih, ix, izm, izd
    integer                         :: idx
    real                            :: h, x, zd, zm
    real                            :: rdx, fx, gx

    do idip=1, ndip
    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 - tg_dip(idip) * x
             !! .. 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,idip) = modl(izm,iang,idip) + data(izd,ih)
                else
                   data(izd,ih) = data(izd,ih) + modl(izm,iang,idip)
                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,idip) = modl(izm,iang,idip) &
                                         + data(izd, ih)*gx + data(izd+1, ih)*fx
                else
                   data(izd  , ih) = data(izd  , ih) + modl(izm,iang,idip)*gx
                   data(izd+1, ih) = data(izd+1, ih) + modl(izm,iang,idip)*fx
                end if
             end if
                
          end do
       end do
    end do
    end do
  end subroutine 
  !! local dip estimation
  subroutine off2ang2d_sd_dip()
    
  end subroutine
end module 






