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

  use sep
  use util_

  implicit none
  character(len=128) :: name
  logical            :: adj

  real,    pointer :: mapCC(:,:)
  real,    pointer :: mapRC(:,:)
  complex, pointer ::  rays(:,:)

  type(myaxis) :: ax,az,at,ag
  integer      :: ix,iz,it,ig, jx,jz
  real         :: x,z

  integer      :: nsz,nsx
  real         :: zz,xx
  real         :: dz,dx
  real         :: r,sincr,dr

  call sep_init()
  call from_param("adj",adj,.false.)
  call from_param("nsz",nsz,11)
  call from_param("nsx",nsx,11)

  name="rays"
  call ainaxis(at,1,name)
  call ainaxis(ag,2,name)

  if(adj) then
     az%l="az"
     call printaxis(az)
     call pinaxis(az)
     call printaxis(az)

     ax%l="ax"
     call printaxis(ax)
     call pinaxis(ax)
     call printaxis(ax)

     call houaxis(az,1)
     call houaxis(ax,2)
  else
     call hinaxis(az,1)
     call hinaxis(ax,2)

     call houaxis(at,1)
     call houaxis(ag,2)
  end if

  call sep_close()
  dr = sqrt(az%d**2 + ax%d**2)

  allocate(mapCC(az%n,ax%n)); mapCC=0.
  allocate(mapRC(at%n,ag%n)); mapRC=0.
  allocate( rays(at%n,ag%n))

  call sreed("rays",rays,size(rays)*8)
  if(adj) then
     call sreed("in" ,mapRC,size(mapRC)*4)
  else
     call sreed("in" ,mapCC,size(mapCC)*4)
  end if

  do ig=1,ag%n
     if(mod(ig,100)==0) write(0,*) ig,ag%n
     do it=1,at%n   
        z=aimag(rays(it,ig))
        x= real(rays(it,ig))

        iz=floor((z-az%o)/az%d)+1
        ix=floor((x-ax%o)/ax%d)+1

        if(iz<   1+nsz) cycle
        if(iz>az%n-nsz) cycle
        if(ix<   1+nsx) cycle
        if(ix>ax%n-nsx) cycle

        do jz=iz-nsz,iz+nsz
           zz=az%o+(jz-1)*az%d
           dz=zz-z
           
           do jx=ix-nsx,ix+nsx
              xx=ax%o+(jx-1)*ax%d
              dx=xx-x
              
              r=sqrt(dz**2 + dx**2)
              if(abs(r)<epsilon(r)) then
                 sincr = 1.0
              else
                 r = r / dr
                 sincr = sin(r)/r
              end if

              if(adj) then
                 mapCC(jz,jx) = mapCC(jz,jx) + mapRC(it,ig)*sincr
              else
                 mapRC(it,ig) = mapRC(it,ig) + mapCC(jz,jx)*sincr
              end if
           end do
        end do

     end do
  end do

  if(adj) then
     call srite("out",mapCC,size(mapCC)*4)
  else
     call srite("out",mapRC,size(mapRC)*4)
  end if

  call exit(0)
end program C2RSINC_






