program CC2RC3D_
  !! interpolates from cartesian coordinates 
  !! to ray coordinates in 3D
  !! Jeff Shragge (jeff@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(:,:,:,:)

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

  real         :: xmin,xmax,zmin,zmax,ymin,ymax
  real         :: fx,fz,fy
  real         :: w000,w001,w010,w011,w100,w101,w110,w111,sumw

  call sep_init()
  call from_param("adj",adj,.false.)

  name="rays"
  call ainaxis(at,1,name)
  call ainaxis(ag,2,name)
  call ainaxis(am,3,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)

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

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

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

  call sep_close()
  xmin = ax%o
  xmax = ax%o + (ax%n-1)*ax%d
  zmin = az%o
  zmax = az%o + (az%n-1)*az%d
  ymin = ay%o
  ymax = ay%o + (ay%n-1)*ay%d

  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,3))

  call sreed("rays",rays,size(rays))

  if(adj) then
     call sreed("in" ,mapRC,size(mapRC)*4)
     write(0,*) 'READ IN RAYS'
  else
     call sreed("in" ,mapCC,size(mapCC)*4)
     write(0,*) 'READ IN CARTESIAN'
  end if

  do im=1,am%n
     do ig=1,ag%n
        do it=1,at%n

           x= rays(it,ig,im,1)
           y= rays(it,ig,im,2)
           z= rays(it,ig,im,3)

           z = min(z,zmax)
           z = max(z,zmin)
           x = min(x,xmax)
           x = max(x,xmin)
           y = min(y,ymax)
           y = max(y,ymin)

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

           iz = min(iz  ,az%n)
           iz = max(iz  , 1  )
           jz = min(iz+1,az%n)

           ix = min(ix  ,ax%n)
           ix = max(ix  , 1  )
           jx = min(ix+1,ax%n)

           iy = min(iy  ,ay%n)
           iy = max(iy  , 1  )
           jy = min(iy+1,ay%n)

           fz= z - az%o-(iz-1)*az%d
           fx= x - ax%o-(ix-1)*ax%d
           fy= y - ay%o-(iy-1)*ay%d
      
           fz = max(0., fz)
           fz = min(1., fz)

           fx = max(0., fx)
           fx = min(1., fx)

           fy = max(0., fy)
           fy = min(1., fy)

           !! zxy   zxy
           !! 000   010     --> x
           !!    .           |
           !! 100   110      v
           !!                z
           w000 = (1.-fz)*(1.-fx)*(1.-fy)
           w001 = (1.-fz)*(1.-fx)*(   fy)
           w010 = (1.-fz)*(   fx)*(1.-fy)
           w100 = (   fz)*(1.-fx)*(1.-fy)
           w101 = (   fz)*(1.-fx)*(   fy)
           w110 = (   fz)*(   fx)*(1.-fy)
           w011 = (1.-fz)*(   fx)*(   fy)
           w111 = (   fz)*(   fx)*(   fy)

           !! normalize the sum of the weights
           sumw = w000+w001+w010+w011+w100+w101+w110+w111
           w000 = w000 / sumw
           w010 = w001 / sumw
           w011 = w010 / sumw
           w011 = w011 / sumw
           w100 = w100 / sumw
           w110 = w101 / sumw
           w111 = w110 / sumw
           w111 = w111 / sumw

           if(adj) then
              mapCC(iz,ix,iy) = mapCC(iz,ix,iy) + w000 * mapRC(it,ig,im)
              mapCC(jz,ix,iy) = mapCC(jz,ix,iy) + w100 * mapRC(it,ig,im)
              mapCC(iz,jx,iy) = mapCC(iz,jx,iy) + w010 * mapRC(it,ig,im)
              mapCC(jz,jx,iy) = mapCC(jz,jx,iy) + w110 * mapRC(it,ig,im)
              mapCC(iz,ix,jy) = mapCC(iz,ix,jy) + w001 * mapRC(it,ig,im)
              mapCC(jz,ix,jy) = mapCC(jz,ix,jy) + w101 * mapRC(it,ig,im)
              mapCC(iz,jx,jy) = mapCC(iz,jx,jy) + w111 * mapRC(it,ig,im)
              mapCC(jz,jx,jy) = mapCC(jz,jx,jy) + w111 * mapRC(it,ig,im)
           else
              mapRC(it,ig,im) = mapRC(it,ig,im) + w000 * mapCC(iz,ix,iy)
              mapRC(it,ig,im) = mapRC(it,ig,im) + w100 * mapCC(jz,ix,iy)
              mapRC(it,ig,im) = mapRC(it,ig,im) + w010 * mapCC(iz,jx,iy)
              mapRC(it,ig,im) = mapRC(it,ig,im) + w110 * mapCC(jz,jx,iy)
              mapRC(it,ig,im) = mapRC(it,ig,im) + w001 * mapCC(iz,ix,jy)
              mapRC(it,ig,im) = mapRC(it,ig,im) + w101 * mapCC(jz,ix,jy)
              mapRC(it,ig,im) = mapRC(it,ig,im) + w011 * mapCC(iz,jx,jy)
              mapRC(it,ig,im) = mapRC(it,ig,im) + w111 * mapCC(jz,jx,jy)
         end if

        end do
     end do
  end do
  if(adj) then
     call srite("out",mapCC,size(mapCC)*4)
     write(0,*) 'WROTE OUT CARTESIAN'
  else
     call srite("out",mapRC,size(mapRC)*4)
     write(0,*) 'WROTE OUT RAYS'
  end if

  call exit(0)
end program CC2RC3D_






