program C2Rsinc3D
  !! 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(:,:,:)
  real,    pointer ::  rays(:,:,:,:,:)
  real,    pointer ::  fold(:,:,:)

  type(myaxis) :: ax,az,ay,at,ag,am
  integer      :: ix,iz,iy,it,ig,im,jx,jy,jz
  real         :: x,z,y

  integer      :: nsz,nsx,nsy,norm
  real         :: zz,xx,yy
  real         :: dz,dx,dy
  real         :: r,sincr,dr

  call sep_init()
  call from_param("adj ",adj,.false.)
  call from_param("nsz ",nsz ,5)
  call from_param("nsx ",nsx ,5)
  call from_param("nsy ",nsy ,5)
  call from_param("norm",norm,1)

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

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

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

     ay%l="ay"
     call pinaxis(ay)

     call houaxis(az,1)
     call houaxis(ax,2)
     call houaxis(ay,3)

  else
     call hinaxis(az,1)
     call hinaxis(ax,2)
     call hinaxis(ay,3)

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

  
  if (ay%n .eq. 1) then
     dr = sqrt(az%d**2 + ax%d**2 + ay%d**2)
  else
     dr = sqrt(az%d**2 + ax%d**2 )
  end if

  allocate(mapCC(az%n,ax%n,ay%n)); mapCC=0.
  allocate(mapRC(at%n,ag%n,am%n)); mapRC=0.
  allocate( rays(at%n,ag%n,am%n,1,3))
  allocate( fold(az%n,ax%n,ay%n)); fold =0.

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

write(0,*) 'X max/min', maxval(rays(:,:,:,1,1)),minval(rays(:,:,:,1,1))
write(0,*) 'Y max/min', maxval(rays(:,:,:,1,2)),minval(rays(:,:,:,1,2))
write(0,*) 'Z max/min', maxval(rays(:,:,:,1,3)),minval(rays(:,:,:,1,3))

  do ig=1,ag%n
     do im=1,am%n
     if(mod(ag%n*(ig-1)+im,100)==0) write(0,*) ag%n*(ig-1)+im,ag%n*am%n
        do it=1,at%n   
           x=rays(it,ig,im,1,1)
           y=rays(it,ig,im,1,2)
           z=rays(it,ig,im,1,3)
           iz=floor((z-az%o)/az%d)+1
           ix=floor((x-ax%o)/ax%d)+1
           iy=floor((y-ay%o)/ay%d)+1

           if(iz<   1+nsz) cycle
           if(iz>az%n-nsz) cycle
           if(ix<   1+nsx) cycle
           if(ix>ax%n-nsx) cycle
           if(iy<   1+nsy) cycle
           if(iy>ay%n-nsy) 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

                 do jy=iy-nsy,iy+nsy
                    yy=ay%o+(jy-1)*ay%d
                    dy=yy-y

                    r=sqrt(dz**2 + dx**2 + dy**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,jy) = mapCC(jz,jx,jy) + mapRC(it,ig,im)*sincr
                        fold(jz,jx,jy) =  fold(jz,jx,jy) + sincr
                    else
                       mapRC(it,ig,im) = mapRC(it,ig,im) + mapCC(jz,jx,jy)*sincr
                    end if
                 end do
              end do
           end do

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

  if(adj) then
   if (norm .eq. 1) then
        where (fold .gt. 0.)
           mapCC=mapCC/fold
        end where
        write(0,*) 'NORMALIZED Data by fold'
        call to_history("n1",az%n,"fold.H")
        call to_history("n2",ax%n,"fold.H")
        call to_history("n3",ay%n,"fold.H")
        call to_history("d1",az%d,"fold.H")
        call to_history("d2",ax%d,"fold.H")
        call to_history("d3",ay%d,"fold.H")
        call to_history("o1",az%o,"fold.H")
        call to_history("o2",ax%o,"fold.H")
        call to_history("o3",ay%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 sep_close()

  call exit(0)
end program C2Rsinc3D






