module C2R
  !! interpolates from cartesian coordinates to ray coordinates
  !! using sinc interpolation
  !! Paul Sava (paul@sep.stanford.edu)
  !! 08/15/03
!  use sep
!  use utils_

  implicit none

  integer,private :: nsx,nsy,nsz,norm,cnx,cny,cnz,rnx,rny,rnz
  real,private  :: xmin,xmax,ymin,ymax,zmin,zmax,dxx,dyy,dzz
  real, allocatable, dimension(:,:,:),private :: fold

contains

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

  subroutine C2R_init(insx,insy,insz,inorm,idxx,idyy,idzz,&
                      icnx,icny,icnz,irnx,irnz,irny,ixmin,iymin,izmin)
    integer :: insx,insy,insz,inorm,icnx,icny,icnz,irnx,irny,irnz
    real    :: idxx,idyy,idzz,ixmin,iymin,izmin
    nsx = insx;  nsy = insy;  nsz = insz
    norm= inorm
    xmin=ixmin; ymin=iymin; zmin=izmin
    cnx = icnx;  cny = icny;  cnz=icnz
    rnx = irnx;  rny = irny;  rnz=irnz
    dxx = idxx;  dyy = idyy;  dzz = idzz

    allocate(fold(cnz,cnx,cny))

  end subroutine C2R_init

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

  subroutine C2R_run(RC,rays,CC)
    real    :: CC(:,:,:),RC(:,:,:)
    real    :: rays(:,:,:,:)
    integer :: ix,iz,iy,it,ig,im,jx,jy,jz
    real    :: x,z,y, zz,xx,yy
    real    :: ddz,ddx,ddy,r,sincr,dr

    write(0,*) 'MIN/MAX Image to interp',minval(RC),maxval(RC)
    write(0,*) 'RAY DIM',rnz,rnx,rny
    write(0,*) 'CART DIM',cnz,cnx,cny

    if (rny .eq. 1) then
       dr = sqrt( dxx**2 + dzz**2 )
    else
       dr = sqrt( dxx**2 + dyy**2 + dzz**2)
    end if

    do im=1,rny
       do ig=1,rnx
          do it=1,rnz   
             x=rays(it,ig,im,1)
             y=rays(it,ig,im,2)
             z=rays(it,ig,im,3)

             iz=floor((z-zmin)/dzz)+1
             ix=floor((x-xmin)/dxx)+1
             iy=max(1,floor((y-ymin)/dyy)+1)

             if(iz<   1+nsz) cycle
             if(iz> cnz-nsz) cycle
             if(ix<   1+nsx) cycle
             if(ix> cnx-nsx) cycle
             if(iy<   1+nsy) cycle
             if(iy> cny-nsy) 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

                   do jy=iy-nsy,iy+nsy
                      yy=ymin+(jy-1)*dyy
                      ddy=yy-y
                      r=sqrt(ddz**2 + ddx**2 + ddy**2)

                      if(abs(r)<epsilon(r)) then
                         sincr = 1.0
                      else
                         r = r / dr
                         sincr = sin(r)/r
                      end if

                      CC(jz,jx,jy) = CC(jz,jx,jy) + RC(it,ig,im)*sincr
                      fold(jz,jx,jy) =  fold(jz,jx,jy) + 1.! sincr

                   end do
                end do
             end do

          end do !! End Tau
       end do !! End Mu
    end do !! End Gamma

    !   if (norm .eq. 1) then
    !      where (fold .gt. 0.)
    CC=CC/(fold+0.0001)
    !      end where
    !   end if

  end subroutine C2R_run

end module C2R