module rotsog_mod1
  ! simple linear Radon transform  z = z0 + p*h

  use adj_mod

  implicit none

  integer, private :: nz 
  integer, private :: nhx,nmx,nax
  real   , private :: z0, dz ,axmin
  real   , private :: mx0, dmx 
  real   , private :: hx0, dhx,dax
  real, dimension(:), pointer, private :: tgama,talfa
  real, parameter, private :: PI = 3.14159265359
contains
    
  subroutine rotsog_init(z0_in, hx0_in, mx0_in, dz_in, dhx_in, dmx_in, nz_in, nhx_in, nmx_in, axmin_in, dax_in, nax_in) 
    
    real    :: z0_in, dz_in ,axmin_in
    real    :: hx0_in, dhx_in , dax_in, ax
    real    :: mx0_in, dmx_in 
    integer :: nz_in, nhx_in, nmx_in, nax_in, iax

    z0  = z0_in ; dz  = dz_in;  mx0  = mx0_in 
    hx0 = hx0_in; dhx = dhx_in; dmx = dmx_in 
    nz  = nz_in ; nhx = nhx_in; nmx = nmx_in
    axmin=axmin_in; dax=dax_in; nax= nax_in

    allocate(tgama(nax),talfa(nax))

!    do iax=1, nax
!       ax = axmin + (iax-1)*dax
!       tgama(iax) = -tan(ax*pi/180)
!       co(iax)=abs(cosd(ax))
!    end do

    do iax=1, 21
       ax = -50. + (iax-1)*2.5
       tgama(iax) = -tan(ax*pi/180)
    end do

    do iax=1, 21
       ax = -20. + (iax-1)*2.
       talfa(iax) = -tan(ax*pi/180)
!       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 rotsog_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 rotsog_exe(adj,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine rotsog_exe(adj,modl,data)
    logical,intent(in)       	:: adj 
    real		  	:: modl(:,:,:)
    real		  	:: data(:,:,:)
    integer                  	:: izm,izd
    integer                  	:: ihx,igx,iax,ix,iz,ixb,ib
    integer                  	:: stat, idx,ihnew
    real                     	:: zd, zm,tteta_s
    real			:: deltax,deltaz
    real                     	:: hx,hnew,cgama,cteta_S
    real                     	:: rdx, fx, gx

     do ihx=1, nhx
	write(0,*) ihx , " from  ", nhx
        hx = hx0 + (ihx-1) * dhx               ! offset
        do iax=1, 21! nax
           do igx=1, 41! nax
              cgama=abs(cos((-50. + (igx-1)*2.5)*pi/180))
              if (1-talfa(iax)*tgama(igx) > tiny(0.)) then
                 tteta_s=-(talfa(iax)+tgama(igx))/(1-talfa(iax)*tgama(igx))
              else
	         tteta_s=sign(huge(1.),-(talfa(iax)+tgama(igx)))
	      end if
              cteta_s=abs(cos(atan(tteta_s)))
              do ix=1,nmx
		 do ib=-4,4
		    ixb=ix+ib
                 do iz=1,nz
                    zm  = z0 + (iz-1) * dz; zd  = zm - 0.5*(tgama(igx) * hx + ib*dmx * talfa(iax))
		    hnew = -(zd - zm)/tteta_s
!		    hnew = -(zd)/tteta_s
!		    hnew = - tgama(igx) * hx/tteta_s
                    ihnew = int((hnew - hx0)/dhx)+1
		    if (ihnew.ge.1.and.ihnew.le.nhx) then
		       if (ixb .ge. 1 .and. ixb .le. nmx) modl(iz,ihnew,ix) = modl(iz,ihnew,ix) + cgama/cteta_s*data(iz,ihx,ixb)
                    end if
                 end do
              end do
              end do
           end do
        end do
     end do

!              !! .. 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(iax+nax*(izm-1)) = modl(iax+nax*(izm-1)) + data(ihx+nhx*(izd-1))
!                 else
!	 	   data(ihx+nhx*(izd-1)) = data(ihx+nhx*(izd-1)) + modl(iax+nax*(izm-1))
!                 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				   ! out of the range
!                 fx = rdx - idx
!                 gx = 1.0 - fx
!                 if (adj) then
!                     modl(izm+nz*(iax-1)) = modl(izm+nz*(iax-1)) &
!                                   + (data(izd+nz*(ihx-1))*gx + data(izd+1+nz*(ihx-1))*fx)
!                else
!                    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
!                 end if
!              end if



!	if (style.eq.0) then
!	  if (mode .eq. 0) then
!             call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 0)
!             stat = off2ang2d1_sd_oper(adj, F, modl, data)
!	  else if (mode .eq. 1) then
!             call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 1)
!             call solver_tiny(m=modl, d=data, Fop=off2ang2d1_sd_oper, stepper=cgstep, niter=niter)
!	  else
!             call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 1)
!	     call solver_prc(m=modl, d=data,Fop=off2ang2d1_sd_oper,Sop=ident_lop,stepper=cgstep,nSop=nAop,niter=niter,eps=eps,&
!	         verb=.true.)
!	  end if
!	else
!           call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 1)
!           stat = off2ang2d1_sd_oper(adj, F, modl, data)
!	end if	
!        call srite("out", modl, 4*nz*nax)
!     end do

!    do iax=1, nax 
!       if (abs(tg(iax)) > dz/dhx) cycle 
!    end do
  end subroutine 
  
end module 
