# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
program CC2RC_
  !! interpolates from cartesian coordinates to ray coordinates
  !! 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

  real :: xmin,xmax,zmin,zmax
  real :: fx,fz
  real :: w00,w01,w10,w11, sumw

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

  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()
  xmin = ax%o
  xmax = ax%o + (ax%n-1)*ax%d
  zmin = az%o
  zmax = az%o + (az%n-1)*az%d

  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
     do it=1,at%n

        z=aimag(rays(it,ig))
        x= real(rays(it,ig))

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

        iz=floor((z-az%o)/az%d)+1
        ix=floor((x-ax%o)/ax%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)

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

        !! 00 01 --> x
        !! . |
        !! 10 11 v
        !! z
        w00=(1.-fx)*(1.-fz)
        w10=( fx)*(1.-fz)
        w01=(1.-fx)*( fz)
        w11=( fx)*( fz)

        !! normalize the sum of the weights
        sumw = w00+w10+w01+w11
        w00 = w00 / sumw
        w10 = w10 / sumw
        w01 = w01 / sumw
        w11 = w11 / sumw

        if(adj) then
           mapCC(iz,ix) = mapCC(iz,ix) + w00 * mapRC(it,ig)
           mapCC(jz,ix) = mapCC(jz,ix) + w10 * mapRC(it,ig)
           mapCC(iz,jx) = mapCC(iz,jx) + w01 * mapRC(it,ig)
           mapCC(jz,jx) = mapCC(jz,jx) + w11 * mapRC(it,ig)
        else
           mapRC(it,ig) = mapRC(it,ig) + w00 * mapCC(iz,ix)
           mapRC(it,ig) = mapRC(it,ig) + w10 * mapCC(jz,ix)
           mapRC(it,ig) = mapRC(it,ig) + w01 * mapCC(iz,jx)
           mapRC(it,ig) = mapRC(it,ig) + w11 * mapCC(jz,jx)
        end if

     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 CC2RC_
