module C2R2D

  implicit none

  integer, private :: nsx,nsz,norm,cnx,cnz,rnx,rnz
  real   , private  :: xmin,xmax,zmin,zmax,dxx,dzz

contains

  !----------------------------------------------------------------  

 subroutine C2R_init(insx,insz,inorm,idxx,idzz,&
                     icnx,icnz,irnx,irnz,ixmin,izmin)
    integer :: insx,insz,inorm,icnx,icnz,irnx,irnz
    real    :: idxx,idzz,ixmin,izmin
    nsx = insx;  nsz = insz;   norm= inorm
    xmin=ixmin;  zmin= izmin
    cnx = icnx;  cnz = icnz
    rnx = irnx;  rnz = irnz
    dxx = idxx;  dzz = idzz
    
    xmax=(cnx-1)*dxx+xmin
    zmax=(cnz-1)*dzz+zmin

  end subroutine C2R_init

  !---------------------------------------------------------------- 

  subroutine C2R_run(RC,rays,CC,gnorm,adj)
    real    :: CC(:,:),RC(:,:),rays(:,:,:),gnorm(:,:)
    real,allocatable:: fold(:,:)
    integer :: ix,iz,it,ig,im,jx,jz
    real    :: x,z,zz,xx
    real    :: ddz,ddx,r,sincr,dr
    logical :: adj

!    real         :: xmin,xmax,zmin,zmax
    real         :: fx,fz
    real         :: w00,w01,w10,w11, sumw


    if (adj) then
       allocate(fold(cnz,cnx))
    else
       allocate(fold(rnz,rnx))
    end if

!    if (adj) write(0,*) 'MIN/MAX Image to interp',minval(RC),maxval(RC)
!    write(0,*) ' RAY DIM',rnz,rnx
!    write(0,*) 'CART DIM',cnz,cnx
!    if (adj) write(0,*) 'MIN/MAX Gnorm',minval(gnorm),maxval(gnorm)
!    if (adj) write(0,*) 'MIN/MAX CC',minval(CC),maxval(CC)
    dr = sqrt( dxx**2 + dzz**2 )

    if (adj) then
       do ig=1,rnx
          do it=1,rnz   
             x=rays(it,ig,1)
             z=rays(it,ig,2)
             
             iz=floor((z-zmin)/dzz)+1
             ix=floor((x-xmin)/dxx)+1
             
             if(iz<   1+nsz) cycle
             if(iz> cnz-nsz) cycle
             if(ix<   1+nsx) cycle
             if(ix> cnx-nsx) cycle
             
             do jz=iz-nsz,iz+nsz
                zz=zmin+(jz-1)*dzz
                ddz=zz-z
                
                do jx=ix-nsx,ix+nsx
                   xx=xmin+(jx-1)*dxx
                   ddx=xx-x
                   
                   r=sqrt( ddz**2 + ddx**2 )
                   
                   if(abs(r)<epsilon(r)) then
                      sincr = 1.0
                   else
                      r = r / dr
                      sincr = sin(r)/r
                   end if
                   
                   CC(jz,jx)=CC(jz,jx)+RC(it,ig)*sincr!*sqrt(abs(gnorm(it,ig)))
                   fold(jz,jx) = fold(jz,jx) + 1.
                   
                end do
             end do
             
          end do !! End Tau
       end do !! End Mu

       where (fold .gt. 0.)
          CC=CC /fold
       end where

    else
       do ig=1,rnx
          do it=1,rnz

             x=rays(it,ig,1)
             z=rays(it,ig,2)

             z = min(z,zmax)
             z = max(z,zmin)
             x = min(x,xmax)
             x = max(x,xmin)
             
             iz=floor((z-zmin)/dzz)+1
             ix=floor((x-xmin)/dxx)+1
             
             iz = min(iz  ,cnz)
             iz = max(iz  , 1  )
             jz = min(iz+1,cnz)
             
             ix = min(ix  ,cnx)
             ix = max(ix  , 1  )
             jx = min(ix+1,cnx)
             
             fz= (z-zmin-(iz-1)*dzz)
             fx= (x-xmin-(ix-1)*dxx)
             fz = max(0., fz)
             fz = min(1., fz)
             fx = max(0., fx)
             fx = min(1., fx)
             
             !! 00   01       --> x
             !!    .         |
             !! 10   11      v
             !!               z
             w00=(1.-fx)*(1.-fz)
             w10=(   fx)*(1.-fz)
             w01=(1.-fx)*(   fz)
             w11=(   fx)*(   fz)
             
             !! normalize the sum of the weights
             sumw = w00+w10+w01+w11
             w00 = w00 / sumw
             w10 = w10 / sumw
             w01 = w01 / sumw
             w11 = w11 / sumw

             RC(it,ig) = RC(it,ig) + w00 * CC(iz,ix)
             RC(it,ig) = RC(it,ig) + w10 * CC(jz,ix)
             RC(it,ig) = RC(it,ig) + w01 * CC(iz,jx)
             RC(it,ig) = RC(it,ig) + w11 * CC(jz,jx)

             w00=0.;w01=0.;w10=0.;w11=0.;sumw=0.
          end do
       end do


    end if

    
    deallocate ( fold )
  end subroutine C2R_run

end module C2R2D






