program C2Rsinc_norm_
  !! 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
  character(len=128) :: name
  logical            :: adj
  real,    pointer   :: mapCC(:,:),mapRC(:,:)
  real,    pointer   :: fold(:,:)
  complex, pointer   ::  rays(:,:)
  type(myaxis) :: ax,az,at,ag
  integer      :: ix,iz,it,ig,jx,jz,nsz,nsx
  real         :: x,z,zz,xx,dz,dx,r,sincr,dr
  integer      :: norm
  call sep_init()
  call from_param("adj",adj,.false.)
  call from_param("nsz",nsz,11)
  call from_param("nsx",nsx,11)
  call from_param("norm",norm,0)
  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()
  allocate(mapCC(az%n,ax%n)); mapCC=0.
  allocate(mapRC(at%n,ag%n)); mapRC=0.
  allocate( fold(az%n,ax%n)); fold =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
  dr = sqrt(az%d**2 + ax%d**2)

  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 .or. iz>az%n-nsz .or. ix<1+nsx .or. 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
                 fold(jz,jx)= fold(jz,jx)+1
              else
                 mapRC(it,ig) = mapRC(it,ig) + mapCC(jz,jx)*sincr
              end if
           end do
        end do

     end do
  end do

  if(adj) then
     if (norm .eq. 1) then
        where (fold .gt. 0)
           mapCC=mapCC/fold
        end where
        call to_history("n1",az%n,"fold.H")
        call to_history("n2",ax%n,"fold.H")
        call to_history("d1",az%d,"fold.H")
        call to_history("d2",ax%d,"fold.H")
        call to_history("o1",az%o,"fold.H")
        call to_history("o2",ax%o,"fold.H")
        call srite("fold.H",fold,size(fold)*4)
     end if
     call srite("out",mapCC,size(mapCC)*4)
  else
     call srite("out",mapRC,size(mapRC)*4)
  end if
  call exit(0)
end program C2Rsinc_norm_






